OSDN Git Service

* config.gcc (i[34567]86-*-mingw32*): Enable threads by default.
[pf3gnuchains/gcc-fork.git] / gcc / fortran / simplify.c
1 /* Simplify intrinsic functions at compile-time.
2    Copyright (C) 2000, 2001, 2002, 2003, 2004 Free Software Foundation,
3    Inc.
4    Contributed by Andy Vaught & Katherine Holcomb
5
6 This file is part of GCC.
7
8 GCC is free software; you can redistribute it and/or modify it under
9 the terms of the GNU General Public License as published by the Free
10 Software Foundation; either version 2, or (at your option) any later
11 version.
12
13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with GCC; see the file COPYING.  If not, write to the Free
20 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
21 02111-1307, USA.  */
22
23 #include "config.h"
24 #include "system.h"
25 #include "flags.h"
26
27 #include <string.h>
28
29 #include "gfortran.h"
30 #include "arith.h"
31 #include "intrinsic.h"
32
33 static mpf_t mpf_zero, mpf_half, mpf_one;
34 static mpz_t mpz_zero;
35
36 gfc_expr gfc_bad_expr;
37
38
39 /* Note that 'simplification' is not just transforming expressions.
40    For functions that are not simplified at compile time, range
41    checking is done if possible.
42
43    The return convention is that each simplification function returns:
44
45      A new expression node corresponding to the simplified arguments.
46      The original arguments are destroyed by the caller, and must not
47      be a part of the new expression.
48
49      NULL pointer indicating that no simplification was possible and
50      the original expression should remain intact.  If the
51      simplification function sets the type and/or the function name
52      via the pointer gfc_simple_expression, then this type is
53      retained.
54
55      An expression pointer to gfc_bad_expr (a static placeholder)
56      indicating that some error has prevented simplification.  For
57      example, sqrt(-1.0).  The error is generated within the function
58      and should be propagated upwards
59
60    By the time a simplification function gets control, it has been
61    decided that the function call is really supposed to be the
62    intrinsic.  No type checking is strictly necessary, since only
63    valid types will be passed on.  On the other hand, a simplification
64    subroutine may have to look at the type of an argument as part of
65    its processing.
66
67    Array arguments are never passed to these subroutines.
68
69    The functions in this file don't have much comment with them, but
70    everything is reasonably straight-forward.  The Standard, chapter 13
71    is the best comment you'll find for this file anyway.  */
72
73 /* Static table for converting non-ascii character sets to ascii.
74    The xascii_table[] is the inverse table.  */
75
76 static int ascii_table[256] = {
77   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',
78   '\b', '\t', '\n', '\v', '\0', '\r', '\0', '\0',
79   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',
80   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',
81   ' ', '!', '\'', '#', '$', '%', '&', '\'',
82   '(', ')', '*', '+', ',', '-', '.', '/',
83   '0', '1', '2', '3', '4', '5', '6', '7',
84   '8', '9', ':', ';', '<', '=', '>', '?',
85   '@', 'A', 'B', 'C', 'D', 'E', 'F', 'G',
86   'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O',
87   'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W',
88   'X', 'Y', 'Z', '[', '\\', ']', '^', '_',
89   '`', 'a', 'b', 'c', 'd', 'e', 'f', 'g',
90   'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
91   'p', 'q', 'r', 's', 't', 'u', 'v', 'w',
92   'x', 'y', 'z', '{', '|', '}', '~', '\?'
93 };
94
95 static int xascii_table[256];
96
97
98 /* Range checks an expression node.  If all goes well, returns the
99    node, otherwise returns &gfc_bad_expr and frees the node.  */
100
101 static gfc_expr *
102 range_check (gfc_expr * result, const char *name)
103 {
104
105   if (gfc_range_check (result) == ARITH_OK)
106     return result;
107
108   gfc_error ("Result of %s overflows its kind at %L", name, &result->where);
109   gfc_free_expr (result);
110   return &gfc_bad_expr;
111 }
112
113
114 /* A helper function that gets an optional and possibly missing
115    kind parameter.  Returns the kind, -1 if something went wrong.  */
116
117 static int
118 get_kind (bt type, gfc_expr * k, const char *name, int default_kind)
119 {
120   int kind;
121
122   if (k == NULL)
123     return default_kind;
124
125   if (k->expr_type != EXPR_CONSTANT)
126     {
127       gfc_error ("KIND parameter of %s at %L must be an initialization "
128                  "expression", name, &k->where);
129
130       return -1;
131     }
132
133   if (gfc_extract_int (k, &kind) != NULL
134       || gfc_validate_kind (type, kind) == -1)
135     {
136
137       gfc_error ("Invalid KIND parameter of %s at %L", name, &k->where);
138       return -1;
139     }
140
141   return kind;
142 }
143
144
145 /********************** Simplification functions *****************************/
146
147 gfc_expr *
148 gfc_simplify_abs (gfc_expr * e)
149 {
150   gfc_expr *result;
151   mpf_t a, b;
152
153   if (e->expr_type != EXPR_CONSTANT)
154     return NULL;
155
156   switch (e->ts.type)
157     {
158     case BT_INTEGER:
159       result = gfc_constant_result (BT_INTEGER, e->ts.kind, &e->where);
160
161       mpz_abs (result->value.integer, e->value.integer);
162
163       result = range_check (result, "IABS");
164       break;
165
166     case BT_REAL:
167       result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
168
169       mpf_abs (result->value.real, e->value.real);
170
171       result = range_check (result, "ABS");
172       break;
173
174     case BT_COMPLEX:
175       result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
176
177       mpf_init (a);
178       mpf_mul (a, e->value.complex.r, e->value.complex.r);
179
180       mpf_init (b);
181       mpf_mul (b, e->value.complex.i, e->value.complex.i);
182
183       mpf_add (a, a, b);
184       mpf_sqrt (result->value.real, a);
185
186       mpf_clear (a);
187       mpf_clear (b);
188
189       result = range_check (result, "CABS");
190       break;
191
192     default:
193       gfc_internal_error ("gfc_simplify_abs(): Bad type");
194     }
195
196   return result;
197 }
198
199
200 gfc_expr *
201 gfc_simplify_achar (gfc_expr * e)
202 {
203   gfc_expr *result;
204   int index;
205
206   if (e->expr_type != EXPR_CONSTANT)
207     return NULL;
208
209   /* We cannot assume that the native character set is ASCII in this
210      function.  */
211   if (gfc_extract_int (e, &index) != NULL || index < 0 || index > 127)
212     {
213       gfc_error ("Extended ASCII not implemented: argument of ACHAR at %L "
214                  "must be between 0 and 127", &e->where);
215       return &gfc_bad_expr;
216     }
217
218   result = gfc_constant_result (BT_CHARACTER, gfc_default_character_kind (),
219                                 &e->where);
220
221   result->value.character.string = gfc_getmem (2);
222
223   result->value.character.length = 1;
224   result->value.character.string[0] = ascii_table[index];
225   result->value.character.string[1] = '\0';     /* For debugger */
226   return result;
227 }
228
229
230 gfc_expr *
231 gfc_simplify_acos (gfc_expr * x)
232 {
233   gfc_expr *result;
234   mpf_t negative, square, term;
235
236   if (x->expr_type != EXPR_CONSTANT)
237     return NULL;
238
239   if (mpf_cmp_si (x->value.real, 1) > 0 || mpf_cmp_si (x->value.real, -1) < 0)
240     {
241       gfc_error ("Argument of ACOS at %L must be between -1 and 1",
242                  &x->where);
243       return &gfc_bad_expr;
244     }
245
246   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
247
248   if (mpf_cmp_si (x->value.real, 1) == 0)
249     {
250       mpf_set_ui (result->value.real, 0);
251       return range_check (result, "ACOS");
252     }
253
254   if (mpf_cmp_si (x->value.real, -1) == 0)
255     {
256       mpf_set (result->value.real, pi);
257       return range_check (result, "ACOS");
258     }
259
260   mpf_init (negative);
261   mpf_init (square);
262   mpf_init (term);
263
264   mpf_pow_ui (square, x->value.real, 2);
265   mpf_ui_sub (term, 1, square);
266   mpf_sqrt (term, term);
267   mpf_div (term, x->value.real, term);
268   mpf_neg (term, term);
269   arctangent (&term, &negative);
270   mpf_add (result->value.real, half_pi, negative);
271
272   mpf_clear (negative);
273   mpf_clear (square);
274   mpf_clear (term);
275
276   return range_check (result, "ACOS");
277 }
278
279
280 gfc_expr *
281 gfc_simplify_adjustl (gfc_expr * e)
282 {
283   gfc_expr *result;
284   int count, i, len;
285   char ch;
286
287   if (e->expr_type != EXPR_CONSTANT)
288     return NULL;
289
290   len = e->value.character.length;
291
292   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
293
294   result->value.character.length = len;
295   result->value.character.string = gfc_getmem (len + 1);
296
297   for (count = 0, i = 0; i < len; ++i)
298     {
299       ch = e->value.character.string[i];
300       if (ch != ' ')
301         break;
302       ++count;
303     }
304
305   for (i = 0; i < len - count; ++i)
306     {
307       result->value.character.string[i] =
308         e->value.character.string[count + i];
309     }
310
311   for (i = len - count; i < len; ++i)
312     {
313       result->value.character.string[i] = ' ';
314     }
315
316   result->value.character.string[len] = '\0';   /* For debugger */
317
318   return result;
319 }
320
321
322 gfc_expr *
323 gfc_simplify_adjustr (gfc_expr * e)
324 {
325   gfc_expr *result;
326   int count, i, len;
327   char ch;
328
329   if (e->expr_type != EXPR_CONSTANT)
330     return NULL;
331
332   len = e->value.character.length;
333
334   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
335
336   result->value.character.length = len;
337   result->value.character.string = gfc_getmem (len + 1);
338
339   for (count = 0, i = len - 1; i >= 0; --i)
340     {
341       ch = e->value.character.string[i];
342       if (ch != ' ')
343         break;
344       ++count;
345     }
346
347   for (i = 0; i < count; ++i)
348     {
349       result->value.character.string[i] = ' ';
350     }
351
352   for (i = count; i < len; ++i)
353     {
354       result->value.character.string[i] =
355         e->value.character.string[i - count];
356     }
357
358   result->value.character.string[len] = '\0';   /* For debugger */
359
360   return result;
361 }
362
363
364 gfc_expr *
365 gfc_simplify_aimag (gfc_expr * e)
366 {
367   gfc_expr *result;
368
369   if (e->expr_type != EXPR_CONSTANT)
370     return NULL;
371
372   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
373   mpf_set (result->value.real, e->value.complex.i);
374
375   return range_check (result, "AIMAG");
376 }
377
378
379 gfc_expr *
380 gfc_simplify_aint (gfc_expr * e, gfc_expr * k)
381 {
382   gfc_expr *rtrunc, *result;
383   int kind;
384
385   kind = get_kind (BT_REAL, k, "AINT", e->ts.kind);
386   if (kind == -1)
387     return &gfc_bad_expr;
388
389   if (e->expr_type != EXPR_CONSTANT)
390     return NULL;
391
392   rtrunc = gfc_copy_expr (e);
393
394   mpf_trunc (rtrunc->value.real, e->value.real);
395
396   result = gfc_real2real (rtrunc, kind);
397   gfc_free_expr (rtrunc);
398
399   return range_check (result, "AINT");
400 }
401
402
403 gfc_expr *
404 gfc_simplify_dint (gfc_expr * e)
405 {
406   gfc_expr *rtrunc, *result;
407
408   if (e->expr_type != EXPR_CONSTANT)
409     return NULL;
410
411   rtrunc = gfc_copy_expr (e);
412
413   mpf_trunc (rtrunc->value.real, e->value.real);
414
415   result = gfc_real2real (rtrunc, gfc_default_double_kind ());
416   gfc_free_expr (rtrunc);
417
418   return range_check (result, "DINT");
419
420 }
421
422
423 gfc_expr *
424 gfc_simplify_anint (gfc_expr * e, gfc_expr * k)
425 {
426   gfc_expr *rtrunc, *result;
427   int kind, cmp;
428
429   kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
430   if (kind == -1)
431     return &gfc_bad_expr;
432
433   if (e->expr_type != EXPR_CONSTANT)
434     return NULL;
435
436   result = gfc_constant_result (e->ts.type, kind, &e->where);
437
438   rtrunc = gfc_copy_expr (e);
439
440   cmp = mpf_cmp_ui (e->value.real, 0);
441
442   if (cmp > 0)
443     {
444       mpf_add (rtrunc->value.real, e->value.real, mpf_half);
445       mpf_trunc (result->value.real, rtrunc->value.real);
446     }
447   else if (cmp < 0)
448     {
449       mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
450       mpf_trunc (result->value.real, rtrunc->value.real);
451     }
452   else
453     mpf_set_ui (result->value.real, 0);
454
455   gfc_free_expr (rtrunc);
456
457   return range_check (result, "ANINT");
458 }
459
460
461 gfc_expr *
462 gfc_simplify_dnint (gfc_expr * e)
463 {
464   gfc_expr *rtrunc, *result;
465   int cmp;
466
467   if (e->expr_type != EXPR_CONSTANT)
468     return NULL;
469
470   result =
471     gfc_constant_result (BT_REAL, gfc_default_double_kind (), &e->where);
472
473   rtrunc = gfc_copy_expr (e);
474
475   cmp = mpf_cmp_ui (e->value.real, 0);
476
477   if (cmp > 0)
478     {
479       mpf_add (rtrunc->value.real, e->value.real, mpf_half);
480       mpf_trunc (result->value.real, rtrunc->value.real);
481     }
482   else if (cmp < 0)
483     {
484       mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
485       mpf_trunc (result->value.real, rtrunc->value.real);
486     }
487   else
488     mpf_set_ui (result->value.real, 0);
489
490   gfc_free_expr (rtrunc);
491
492   return range_check (result, "DNINT");
493 }
494
495
496 gfc_expr *
497 gfc_simplify_asin (gfc_expr * x)
498 {
499   gfc_expr *result;
500   mpf_t negative, square, term;
501
502   if (x->expr_type != EXPR_CONSTANT)
503     return NULL;
504
505   if (mpf_cmp_si (x->value.real, 1) > 0 || mpf_cmp_si (x->value.real, -1) < 0)
506     {
507       gfc_error ("Argument of ASIN at %L must be between -1 and 1",
508                  &x->where);
509       return &gfc_bad_expr;
510     }
511
512   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
513
514   if (mpf_cmp_si (x->value.real, 1) == 0)
515     {
516       mpf_set (result->value.real, half_pi);
517       return range_check (result, "ASIN");
518     }
519
520   if (mpf_cmp_si (x->value.real, -1) == 0)
521     {
522       mpf_init (negative);
523       mpf_neg (negative, half_pi);
524       mpf_set (result->value.real, negative);
525       mpf_clear (negative);
526       return range_check (result, "ASIN");
527     }
528
529   mpf_init (square);
530   mpf_init (term);
531
532   mpf_pow_ui (square, x->value.real, 2);
533   mpf_ui_sub (term, 1, square);
534   mpf_sqrt (term, term);
535   mpf_div (term, x->value.real, term);
536   arctangent (&term, &result->value.real);
537
538   mpf_clear (square);
539   mpf_clear (term);
540
541   return range_check (result, "ASIN");
542 }
543
544
545 gfc_expr *
546 gfc_simplify_atan (gfc_expr * x)
547 {
548   gfc_expr *result;
549
550   if (x->expr_type != EXPR_CONSTANT)
551     return NULL;
552
553   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
554
555   arctangent (&x->value.real, &result->value.real);
556
557   return range_check (result, "ATAN");
558
559 }
560
561
562 gfc_expr *
563 gfc_simplify_atan2 (gfc_expr * y, gfc_expr * x)
564 {
565   gfc_expr *result;
566
567   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
568     return NULL;
569
570   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
571
572
573   if (mpf_sgn (y->value.real) == 0 && mpf_sgn (x->value.real) == 0)
574     {
575       gfc_error
576         ("If first argument of ATAN2 %L is zero, the second argument "
577           "must not be zero", &x->where);
578       gfc_free_expr (result);
579       return &gfc_bad_expr;
580     }
581
582   arctangent2 (&y->value.real, &x->value.real, &result->value.real);
583
584   return range_check (result, "ATAN2");
585
586 }
587
588
589 gfc_expr *
590 gfc_simplify_bit_size (gfc_expr * e)
591 {
592   gfc_expr *result;
593   int i;
594
595   i = gfc_validate_kind (e->ts.type, e->ts.kind);
596   if (i == -1)
597     gfc_internal_error ("In gfc_simplify_bit_size(): Bad kind");
598
599   result = gfc_constant_result (BT_INTEGER, e->ts.kind, &e->where);
600   mpz_set_ui (result->value.integer, gfc_integer_kinds[i].bit_size);
601
602   return result;
603 }
604
605
606 gfc_expr *
607 gfc_simplify_btest (gfc_expr * e, gfc_expr * bit)
608 {
609   int b;
610
611   if (e->expr_type != EXPR_CONSTANT || bit->expr_type != EXPR_CONSTANT)
612     return NULL;
613
614   if (gfc_extract_int (bit, &b) != NULL || b < 0)
615     return gfc_logical_expr (0, &e->where);
616
617   return gfc_logical_expr (mpz_tstbit (e->value.integer, b), &e->where);
618 }
619
620
621 gfc_expr *
622 gfc_simplify_ceiling (gfc_expr * e, gfc_expr * k)
623 {
624   gfc_expr *ceil, *result;
625   int kind;
626
627   kind = get_kind (BT_REAL, k, "CEILING", gfc_default_real_kind ());
628   if (kind == -1)
629     return &gfc_bad_expr;
630
631   if (e->expr_type != EXPR_CONSTANT)
632     return NULL;
633
634   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
635
636   ceil = gfc_copy_expr (e);
637
638   mpf_ceil (ceil->value.real, e->value.real);
639   mpz_set_f (result->value.integer, ceil->value.real);
640
641   gfc_free_expr (ceil);
642
643   return range_check (result, "CEILING");
644 }
645
646
647 gfc_expr *
648 gfc_simplify_char (gfc_expr * e, gfc_expr * k)
649 {
650   gfc_expr *result;
651   int c, kind;
652
653   kind = get_kind (BT_CHARACTER, k, "CHAR", gfc_default_character_kind ());
654   if (kind == -1)
655     return &gfc_bad_expr;
656
657   if (e->expr_type != EXPR_CONSTANT)
658     return NULL;
659
660   if (gfc_extract_int (e, &c) != NULL || c < 0 || c > 255)
661     {
662       gfc_error ("Bad character in CHAR function at %L", &e->where);
663       return &gfc_bad_expr;
664     }
665
666   result = gfc_constant_result (BT_CHARACTER, kind, &e->where);
667
668   result->value.character.length = 1;
669   result->value.character.string = gfc_getmem (2);
670
671   result->value.character.string[0] = c;
672   result->value.character.string[1] = '\0';     /* For debugger */
673
674   return result;
675 }
676
677
678 /* Common subroutine for simplifying CMPLX and DCMPLX.  */
679
680 static gfc_expr *
681 simplify_cmplx (const char *name, gfc_expr * x, gfc_expr * y, int kind)
682 {
683   gfc_expr *result;
684
685   result = gfc_constant_result (BT_COMPLEX, kind, &x->where);
686
687   mpf_set_ui (result->value.complex.i, 0);
688
689   switch (x->ts.type)
690     {
691     case BT_INTEGER:
692       mpf_set_z (result->value.complex.r, x->value.integer);
693       break;
694
695     case BT_REAL:
696       mpf_set (result->value.complex.r, x->value.real);
697       break;
698
699     case BT_COMPLEX:
700       mpf_set (result->value.complex.r, x->value.complex.r);
701       mpf_set (result->value.complex.i, x->value.complex.i);
702       break;
703
704     default:
705       gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (x)");
706     }
707
708   if (y != NULL)
709     {
710       switch (y->ts.type)
711         {
712         case BT_INTEGER:
713           mpf_set_z (result->value.complex.i, y->value.integer);
714           break;
715
716         case BT_REAL:
717           mpf_set (result->value.complex.i, y->value.real);
718           break;
719
720         default:
721           gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (y)");
722         }
723     }
724
725   return range_check (result, name);
726 }
727
728
729 gfc_expr *
730 gfc_simplify_cmplx (gfc_expr * x, gfc_expr * y, gfc_expr * k)
731 {
732   int kind;
733
734   if (x->expr_type != EXPR_CONSTANT
735       || (y != NULL && y->expr_type != EXPR_CONSTANT))
736     return NULL;
737
738   kind = get_kind (BT_REAL, k, "CMPLX", gfc_default_real_kind ());
739   if (kind == -1)
740     return &gfc_bad_expr;
741
742   return simplify_cmplx ("CMPLX", x, y, kind);
743 }
744
745
746 gfc_expr *
747 gfc_simplify_conjg (gfc_expr * e)
748 {
749   gfc_expr *result;
750
751   if (e->expr_type != EXPR_CONSTANT)
752     return NULL;
753
754   result = gfc_copy_expr (e);
755   mpf_neg (result->value.complex.i, result->value.complex.i);
756
757   return range_check (result, "CONJG");
758 }
759
760
761 gfc_expr *
762 gfc_simplify_cos (gfc_expr * x)
763 {
764   gfc_expr *result;
765   mpf_t xp, xq;
766
767   if (x->expr_type != EXPR_CONSTANT)
768     return NULL;
769
770   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
771
772   switch (x->ts.type)
773     {
774     case BT_REAL:
775       cosine (&x->value.real, &result->value.real);
776       break;
777     case BT_COMPLEX:
778       mpf_init (xp);
779       mpf_init (xq);
780
781       cosine (&x->value.complex.r, &xp);
782       hypercos (&x->value.complex.i, &xq);
783       mpf_mul (result->value.complex.r, xp, xq);
784
785       sine (&x->value.complex.r, &xp);
786       hypersine (&x->value.complex.i, &xq);
787       mpf_mul (xp, xp, xq);
788       mpf_neg (result->value.complex.i, xp);
789
790       mpf_clear (xp);
791       mpf_clear (xq);
792       break;
793     default:
794       gfc_internal_error ("in gfc_simplify_cos(): Bad type");
795     }
796
797   return range_check (result, "COS");
798
799 }
800
801
802 gfc_expr *
803 gfc_simplify_cosh (gfc_expr * x)
804 {
805   gfc_expr *result;
806
807   if (x->expr_type != EXPR_CONSTANT)
808     return NULL;
809
810   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
811
812   hypercos (&x->value.real, &result->value.real);
813
814   return range_check (result, "COSH");
815 }
816
817
818 gfc_expr *
819 gfc_simplify_dcmplx (gfc_expr * x, gfc_expr * y)
820 {
821
822   if (x->expr_type != EXPR_CONSTANT
823       || (y != NULL && y->expr_type != EXPR_CONSTANT))
824     return NULL;
825
826   return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind ());
827 }
828
829
830 gfc_expr *
831 gfc_simplify_dble (gfc_expr * e)
832 {
833   gfc_expr *result;
834
835   if (e->expr_type != EXPR_CONSTANT)
836     return NULL;
837
838   switch (e->ts.type)
839     {
840     case BT_INTEGER:
841       result = gfc_int2real (e, gfc_default_double_kind ());
842       break;
843
844     case BT_REAL:
845       result = gfc_real2real (e, gfc_default_double_kind ());
846       break;
847
848     case BT_COMPLEX:
849       result = gfc_complex2real (e, gfc_default_double_kind ());
850       break;
851
852     default:
853       gfc_internal_error ("gfc_simplify_dble(): bad type at %L", &e->where);
854     }
855
856   return range_check (result, "DBLE");
857 }
858
859
860 gfc_expr *
861 gfc_simplify_digits (gfc_expr * x)
862 {
863   int i, digits;
864
865   i = gfc_validate_kind (x->ts.type, x->ts.kind);
866   if (i == -1)
867     goto bad;
868
869   switch (x->ts.type)
870     {
871     case BT_INTEGER:
872       digits = gfc_integer_kinds[i].digits;
873       break;
874
875     case BT_REAL:
876     case BT_COMPLEX:
877       digits = gfc_real_kinds[i].digits;
878       break;
879
880     default:
881     bad:
882       gfc_internal_error ("gfc_simplify_digits(): Bad type");
883     }
884
885   return gfc_int_expr (digits);
886 }
887
888
889 gfc_expr *
890 gfc_simplify_dim (gfc_expr * x, gfc_expr * y)
891 {
892   gfc_expr *result;
893
894   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
895     return NULL;
896
897   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
898
899   switch (x->ts.type)
900     {
901     case BT_INTEGER:
902       if (mpz_cmp (x->value.integer, y->value.integer) > 0)
903         mpz_sub (result->value.integer, x->value.integer, y->value.integer);
904       else
905         mpz_set (result->value.integer, mpz_zero);
906
907       break;
908
909     case BT_REAL:
910       if (mpf_cmp (x->value.real, y->value.real) > 0)
911         mpf_sub (result->value.real, x->value.real, y->value.real);
912       else
913         mpf_set (result->value.real, mpf_zero);
914
915       break;
916
917     default:
918       gfc_internal_error ("gfc_simplify_dim(): Bad type");
919     }
920
921   return range_check (result, "DIM");
922 }
923
924
925 gfc_expr *
926 gfc_simplify_dprod (gfc_expr * x, gfc_expr * y)
927 {
928   gfc_expr *mult1, *mult2, *result;
929
930   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
931     return NULL;
932
933   result =
934     gfc_constant_result (BT_REAL, gfc_default_double_kind (), &x->where);
935
936   mult1 = gfc_real2real (x, gfc_default_double_kind ());
937   mult2 = gfc_real2real (y, gfc_default_double_kind ());
938
939   mpf_mul (result->value.real, mult1->value.real, mult2->value.real);
940
941   gfc_free_expr (mult1);
942   gfc_free_expr (mult2);
943
944   return range_check (result, "DPROD");
945 }
946
947
948 gfc_expr *
949 gfc_simplify_epsilon (gfc_expr * e)
950 {
951   gfc_expr *result;
952   int i;
953
954   i = gfc_validate_kind (e->ts.type, e->ts.kind);
955   if (i == -1)
956     gfc_internal_error ("gfc_simplify_epsilon(): Bad kind");
957
958   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
959
960   mpf_set (result->value.real, gfc_real_kinds[i].epsilon);
961
962   return range_check (result, "EPSILON");
963 }
964
965
966 gfc_expr *
967 gfc_simplify_exp (gfc_expr * x)
968 {
969   gfc_expr *result;
970   mpf_t xp, xq;
971   double ln2, absval, rhuge;
972
973   if (x->expr_type != EXPR_CONSTANT)
974     return NULL;
975
976   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
977
978   /* Exactitude doesn't matter here */
979   ln2 = .6931472;
980   rhuge = ln2 * mpz_get_d (gfc_integer_kinds[0].huge);
981
982   switch (x->ts.type)
983     {
984     case BT_REAL:
985       absval = mpf_get_d (x->value.real);
986       if (absval < 0)
987         absval = -absval;
988       if (absval > rhuge)
989         {
990           /* Underflow (set arg to zero) if x is negative and its
991              magnitude is greater than the maximum C long int times
992              ln2, because the exponential method in arith.c will fail
993              for such values.  */
994           if (mpf_cmp_ui (x->value.real, 0) < 0)
995             {
996               if (pedantic == 1)
997                 gfc_warning_now
998                   ("Argument of EXP at %L is negative and too large, "
999                    "setting result to zero", &x->where);
1000               mpf_set_ui (result->value.real, 0);
1001               return range_check (result, "EXP");
1002             }
1003           /* Overflow if magnitude of x is greater than C long int
1004              huge times ln2.  */
1005           else
1006             {
1007               gfc_error ("Argument of EXP at %L too large", &x->where);
1008               gfc_free_expr (result);
1009               return &gfc_bad_expr;
1010             }
1011         }
1012       exponential (&x->value.real, &result->value.real);
1013       break;
1014
1015     case BT_COMPLEX:
1016       /* Using Euler's formula.  */
1017       absval = mpf_get_d (x->value.complex.r);
1018       if (absval < 0)
1019         absval = -absval;
1020       if (absval > rhuge)
1021         {
1022           if (mpf_cmp_ui (x->value.complex.r, 0) < 0)
1023             {
1024               if (pedantic == 1)
1025                 gfc_warning_now
1026                   ("Real part of argument of EXP at %L is negative "
1027                    "and too large, setting result to zero", &x->where);
1028
1029               mpf_set_ui (result->value.complex.r, 0);
1030               mpf_set_ui (result->value.complex.i, 0);
1031               return range_check (result, "EXP");
1032             }
1033           else
1034             {
1035               gfc_error ("Real part of argument of EXP at %L too large",
1036                          &x->where);
1037               gfc_free_expr (result);
1038               return &gfc_bad_expr;
1039             }
1040         }
1041       mpf_init (xp);
1042       mpf_init (xq);
1043       exponential (&x->value.complex.r, &xq);
1044       cosine (&x->value.complex.i, &xp);
1045       mpf_mul (result->value.complex.r, xq, xp);
1046       sine (&x->value.complex.i, &xp);
1047       mpf_mul (result->value.complex.i, xq, xp);
1048       mpf_clear (xp);
1049       mpf_clear (xq);
1050       break;
1051
1052     default:
1053       gfc_internal_error ("in gfc_simplify_exp(): Bad type");
1054     }
1055
1056   return range_check (result, "EXP");
1057 }
1058
1059
1060 gfc_expr *
1061 gfc_simplify_exponent (gfc_expr * x)
1062 {
1063   mpf_t i2, absv, ln2, lnx;
1064   gfc_expr *result;
1065
1066   if (x->expr_type != EXPR_CONSTANT)
1067     return NULL;
1068
1069   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
1070                                 &x->where);
1071
1072   if (mpf_cmp (x->value.real, mpf_zero) == 0)
1073     {
1074       mpz_set_ui (result->value.integer, 0);
1075       return result;
1076     }
1077
1078   mpf_init_set_ui (i2, 2);
1079   mpf_init (absv);
1080   mpf_init (ln2);
1081   mpf_init (lnx);
1082
1083   natural_logarithm (&i2, &ln2);
1084
1085   mpf_abs (absv, x->value.real);
1086   natural_logarithm (&absv, &lnx);
1087
1088   mpf_div (lnx, lnx, ln2);
1089   mpf_trunc (lnx, lnx);
1090   mpf_add_ui (lnx, lnx, 1);
1091   mpz_set_f (result->value.integer, lnx);
1092
1093   mpf_clear (i2);
1094   mpf_clear (ln2);
1095   mpf_clear (lnx);
1096   mpf_clear (absv);
1097
1098   return range_check (result, "EXPONENT");
1099 }
1100
1101
1102 gfc_expr *
1103 gfc_simplify_float (gfc_expr * a)
1104 {
1105   gfc_expr *result;
1106
1107   if (a->expr_type != EXPR_CONSTANT)
1108     return NULL;
1109
1110   result = gfc_int2real (a, gfc_default_real_kind ());
1111   return range_check (result, "FLOAT");
1112 }
1113
1114
1115 gfc_expr *
1116 gfc_simplify_floor (gfc_expr * e, gfc_expr * k)
1117 {
1118   gfc_expr *result;
1119   mpf_t floor;
1120   int kind;
1121
1122   kind = get_kind (BT_REAL, k, "FLOOR", gfc_default_real_kind ());
1123   if (kind == -1)
1124     gfc_internal_error ("gfc_simplify_floor(): Bad kind");
1125
1126   if (e->expr_type != EXPR_CONSTANT)
1127     return NULL;
1128
1129   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
1130
1131   mpf_init (floor);
1132   mpf_floor (floor, e->value.real);
1133   mpz_set_f (result->value.integer, floor);
1134   mpf_clear (floor);
1135
1136   return range_check (result, "FLOOR");
1137 }
1138
1139
1140 gfc_expr *
1141 gfc_simplify_fraction (gfc_expr * x)
1142 {
1143   gfc_expr *result;
1144   mpf_t i2, absv, ln2, lnx, pow2;
1145   unsigned long exp2;
1146
1147   if (x->expr_type != EXPR_CONSTANT)
1148     return NULL;
1149
1150   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
1151
1152   if (mpf_cmp (x->value.real, mpf_zero) == 0)
1153     {
1154       mpf_set (result->value.real, mpf_zero);
1155       return result;
1156     }
1157
1158   mpf_init_set_ui (i2, 2);
1159   mpf_init (absv);
1160   mpf_init (ln2);
1161   mpf_init (lnx);
1162   mpf_init (pow2);
1163
1164   natural_logarithm (&i2, &ln2);
1165
1166   mpf_abs (absv, x->value.real);
1167   natural_logarithm (&absv, &lnx);
1168
1169   mpf_div (lnx, lnx, ln2);
1170   mpf_trunc (lnx, lnx);
1171   mpf_add_ui (lnx, lnx, 1);
1172
1173   exp2 = (unsigned long) mpf_get_d (lnx);
1174   mpf_pow_ui (pow2, i2, exp2);
1175
1176   mpf_div (result->value.real, absv, pow2);
1177
1178   mpf_clear (i2);
1179   mpf_clear (ln2);
1180   mpf_clear (absv);
1181   mpf_clear (lnx);
1182   mpf_clear (pow2);
1183
1184   return range_check (result, "FRACTION");
1185 }
1186
1187
1188 gfc_expr *
1189 gfc_simplify_huge (gfc_expr * e)
1190 {
1191   gfc_expr *result;
1192   int i;
1193
1194   i = gfc_validate_kind (e->ts.type, e->ts.kind);
1195   if (i == -1)
1196     goto bad_type;
1197
1198   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
1199
1200   switch (e->ts.type)
1201     {
1202     case BT_INTEGER:
1203       mpz_set (result->value.integer, gfc_integer_kinds[i].huge);
1204       break;
1205
1206     case BT_REAL:
1207       mpf_set (result->value.real, gfc_real_kinds[i].huge);
1208       break;
1209
1210     bad_type:
1211     default:
1212       gfc_internal_error ("gfc_simplify_huge(): Bad type");
1213     }
1214
1215   return result;
1216 }
1217
1218
1219 gfc_expr *
1220 gfc_simplify_iachar (gfc_expr * e)
1221 {
1222   gfc_expr *result;
1223   int index;
1224
1225   if (e->expr_type != EXPR_CONSTANT)
1226     return NULL;
1227
1228   if (e->value.character.length != 1)
1229     {
1230       gfc_error ("Argument of IACHAR at %L must be of length one", &e->where);
1231       return &gfc_bad_expr;
1232     }
1233
1234   index = xascii_table[(int) e->value.character.string[0] & 0xFF];
1235
1236   result = gfc_int_expr (index);
1237   result->where = e->where;
1238
1239   return range_check (result, "IACHAR");
1240 }
1241
1242
1243 gfc_expr *
1244 gfc_simplify_iand (gfc_expr * x, gfc_expr * y)
1245 {
1246   gfc_expr *result;
1247
1248   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1249     return NULL;
1250
1251   result = gfc_constant_result (BT_INTEGER, x->ts.kind, &x->where);
1252
1253   mpz_and (result->value.integer, x->value.integer, y->value.integer);
1254
1255   return range_check (result, "IAND");
1256 }
1257
1258
1259 gfc_expr *
1260 gfc_simplify_ibclr (gfc_expr * x, gfc_expr * y)
1261 {
1262   gfc_expr *result;
1263   int k, pos;
1264
1265   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1266     return NULL;
1267
1268   if (gfc_extract_int (y, &pos) != NULL || pos < 0)
1269     {
1270       gfc_error ("Invalid second argument of IBCLR at %L", &y->where);
1271       return &gfc_bad_expr;
1272     }
1273
1274   k = gfc_validate_kind (x->ts.type, x->ts.kind);
1275   if (k == -1)
1276     gfc_internal_error ("gfc_simplify_ibclr(): Bad kind");
1277
1278   if (pos > gfc_integer_kinds[k].bit_size)
1279     {
1280       gfc_error ("Second argument of IBCLR exceeds bit size at %L",
1281                  &y->where);
1282       return &gfc_bad_expr;
1283     }
1284
1285   result = gfc_copy_expr (x);
1286
1287   mpz_clrbit (result->value.integer, pos);
1288   return range_check (result, "IBCLR");
1289 }
1290
1291
1292 gfc_expr *
1293 gfc_simplify_ibits (gfc_expr * x, gfc_expr * y, gfc_expr * z)
1294 {
1295   gfc_expr *result;
1296   int pos, len;
1297   int i, k, bitsize;
1298   int *bits;
1299
1300   if (x->expr_type != EXPR_CONSTANT
1301       || y->expr_type != EXPR_CONSTANT
1302       || z->expr_type != EXPR_CONSTANT)
1303     return NULL;
1304
1305   if (gfc_extract_int (y, &pos) != NULL || pos < 0)
1306     {
1307       gfc_error ("Invalid second argument of IBITS at %L", &y->where);
1308       return &gfc_bad_expr;
1309     }
1310
1311   if (gfc_extract_int (z, &len) != NULL || len < 0)
1312     {
1313       gfc_error ("Invalid third argument of IBITS at %L", &z->where);
1314       return &gfc_bad_expr;
1315     }
1316
1317   k = gfc_validate_kind (BT_INTEGER, x->ts.kind);
1318   if (k == -1)
1319     gfc_internal_error ("gfc_simplify_ibits(): Bad kind");
1320
1321   bitsize = gfc_integer_kinds[k].bit_size;
1322
1323   if (pos + len > bitsize)
1324     {
1325       gfc_error
1326         ("Sum of second and third arguments of IBITS exceeds bit size "
1327          "at %L", &y->where);
1328       return &gfc_bad_expr;
1329     }
1330
1331   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
1332
1333   bits = gfc_getmem (bitsize * sizeof (int));
1334
1335   for (i = 0; i < bitsize; i++)
1336     bits[i] = 0;
1337
1338   for (i = 0; i < len; i++)
1339     bits[i] = mpz_tstbit (x->value.integer, i + pos);
1340
1341   for (i = 0; i < bitsize; i++)
1342     {
1343       if (bits[i] == 0)
1344         {
1345           mpz_clrbit (result->value.integer, i);
1346         }
1347       else if (bits[i] == 1)
1348         {
1349           mpz_setbit (result->value.integer, i);
1350         }
1351       else
1352         {
1353           gfc_internal_error ("IBITS: Bad bit");
1354         }
1355     }
1356
1357   gfc_free (bits);
1358
1359   return range_check (result, "IBITS");
1360 }
1361
1362
1363 gfc_expr *
1364 gfc_simplify_ibset (gfc_expr * x, gfc_expr * y)
1365 {
1366   gfc_expr *result;
1367   int k, pos;
1368
1369   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1370     return NULL;
1371
1372   if (gfc_extract_int (y, &pos) != NULL || pos < 0)
1373     {
1374       gfc_error ("Invalid second argument of IBSET at %L", &y->where);
1375       return &gfc_bad_expr;
1376     }
1377
1378   k = gfc_validate_kind (x->ts.type, x->ts.kind);
1379   if (k == -1)
1380     gfc_internal_error ("gfc_simplify_ibset(): Bad kind");
1381
1382   if (pos > gfc_integer_kinds[k].bit_size)
1383     {
1384       gfc_error ("Second argument of IBSET exceeds bit size at %L",
1385                  &y->where);
1386       return &gfc_bad_expr;
1387     }
1388
1389   result = gfc_copy_expr (x);
1390
1391   mpz_setbit (result->value.integer, pos);
1392   return range_check (result, "IBSET");
1393 }
1394
1395
1396 gfc_expr *
1397 gfc_simplify_ichar (gfc_expr * e)
1398 {
1399   gfc_expr *result;
1400   int index;
1401
1402   if (e->expr_type != EXPR_CONSTANT)
1403     return NULL;
1404
1405   if (e->value.character.length != 1)
1406     {
1407       gfc_error ("Argument of ICHAR at %L must be of length one", &e->where);
1408       return &gfc_bad_expr;
1409     }
1410
1411   index = (int) e->value.character.string[0];
1412
1413   if (index < CHAR_MIN || index > CHAR_MAX)
1414     {
1415       gfc_error ("Argument of ICHAR at %L out of range of this processor",
1416                  &e->where);
1417       return &gfc_bad_expr;
1418     }
1419
1420   result = gfc_int_expr (index);
1421   result->where = e->where;
1422   return range_check (result, "ICHAR");
1423 }
1424
1425
1426 gfc_expr *
1427 gfc_simplify_ieor (gfc_expr * x, gfc_expr * y)
1428 {
1429   gfc_expr *result;
1430
1431   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1432     return NULL;
1433
1434   result = gfc_constant_result (BT_INTEGER, x->ts.kind, &x->where);
1435
1436   mpz_xor (result->value.integer, x->value.integer, y->value.integer);
1437
1438   return range_check (result, "IEOR");
1439 }
1440
1441
1442 gfc_expr *
1443 gfc_simplify_index (gfc_expr * x, gfc_expr * y, gfc_expr * b)
1444 {
1445   gfc_expr *result;
1446   int back, len, lensub;
1447   int i, j, k, count, index = 0, start;
1448
1449   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1450     return NULL;
1451
1452   if (b != NULL && b->value.logical != 0)
1453     back = 1;
1454   else
1455     back = 0;
1456
1457   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
1458                                 &x->where);
1459
1460   len = x->value.character.length;
1461   lensub = y->value.character.length;
1462
1463   if (len < lensub)
1464     {
1465       mpz_set_si (result->value.integer, 0);
1466       return result;
1467     }
1468
1469   if (back == 0)
1470     {
1471
1472       if (lensub == 0)
1473         {
1474           mpz_set_si (result->value.integer, 1);
1475           return result;
1476         }
1477       else if (lensub == 1)
1478         {
1479           for (i = 0; i < len; i++)
1480             {
1481               for (j = 0; j < lensub; j++)
1482                 {
1483                   if (y->value.character.string[j] ==
1484                       x->value.character.string[i])
1485                     {
1486                       index = i + 1;
1487                       goto done;
1488                     }
1489                 }
1490             }
1491         }
1492       else
1493         {
1494           for (i = 0; i < len; i++)
1495             {
1496               for (j = 0; j < lensub; j++)
1497                 {
1498                   if (y->value.character.string[j] ==
1499                       x->value.character.string[i])
1500                     {
1501                       start = i;
1502                       count = 0;
1503
1504                       for (k = 0; k < lensub; k++)
1505                         {
1506                           if (y->value.character.string[k] ==
1507                               x->value.character.string[k + start])
1508                             count++;
1509                         }
1510
1511                       if (count == lensub)
1512                         {
1513                           index = start + 1;
1514                           goto done;
1515                         }
1516                     }
1517                 }
1518             }
1519         }
1520
1521     }
1522   else
1523     {
1524
1525       if (lensub == 0)
1526         {
1527           mpz_set_si (result->value.integer, len + 1);
1528           return result;
1529         }
1530       else if (lensub == 1)
1531         {
1532           for (i = 0; i < len; i++)
1533             {
1534               for (j = 0; j < lensub; j++)
1535                 {
1536                   if (y->value.character.string[j] ==
1537                       x->value.character.string[len - i])
1538                     {
1539                       index = len - i + 1;
1540                       goto done;
1541                     }
1542                 }
1543             }
1544         }
1545       else
1546         {
1547           for (i = 0; i < len; i++)
1548             {
1549               for (j = 0; j < lensub; j++)
1550                 {
1551                   if (y->value.character.string[j] ==
1552                       x->value.character.string[len - i])
1553                     {
1554                       start = len - i;
1555                       if (start <= len - lensub)
1556                         {
1557                           count = 0;
1558                           for (k = 0; k < lensub; k++)
1559                             if (y->value.character.string[k] ==
1560                                 x->value.character.string[k + start])
1561                               count++;
1562
1563                           if (count == lensub)
1564                             {
1565                               index = start + 1;
1566                               goto done;
1567                             }
1568                         }
1569                       else
1570                         {
1571                           continue;
1572                         }
1573                     }
1574                 }
1575             }
1576         }
1577     }
1578
1579 done:
1580   mpz_set_si (result->value.integer, index);
1581   return range_check (result, "INDEX");
1582 }
1583
1584
1585 gfc_expr *
1586 gfc_simplify_int (gfc_expr * e, gfc_expr * k)
1587 {
1588   gfc_expr *rpart, *rtrunc, *result;
1589   int kind;
1590
1591   kind = get_kind (BT_REAL, k, "INT", gfc_default_real_kind ());
1592   if (kind == -1)
1593     return &gfc_bad_expr;
1594
1595   if (e->expr_type != EXPR_CONSTANT)
1596     return NULL;
1597
1598   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
1599
1600   switch (e->ts.type)
1601     {
1602     case BT_INTEGER:
1603       mpz_set (result->value.integer, e->value.integer);
1604       break;
1605
1606     case BT_REAL:
1607       rtrunc = gfc_copy_expr (e);
1608       mpf_trunc (rtrunc->value.real, e->value.real);
1609       mpz_set_f (result->value.integer, rtrunc->value.real);
1610       gfc_free_expr (rtrunc);
1611       break;
1612
1613     case BT_COMPLEX:
1614       rpart = gfc_complex2real (e, kind);
1615       rtrunc = gfc_copy_expr (rpart);
1616       mpf_trunc (rtrunc->value.real, rpart->value.real);
1617       mpz_set_f (result->value.integer, rtrunc->value.real);
1618       gfc_free_expr (rpart);
1619       gfc_free_expr (rtrunc);
1620       break;
1621
1622     default:
1623       gfc_error ("Argument of INT at %L is not a valid type", &e->where);
1624       gfc_free_expr (result);
1625       return &gfc_bad_expr;
1626     }
1627
1628   return range_check (result, "INT");
1629 }
1630
1631
1632 gfc_expr *
1633 gfc_simplify_ifix (gfc_expr * e)
1634 {
1635   gfc_expr *rtrunc, *result;
1636
1637   if (e->expr_type != EXPR_CONSTANT)
1638     return NULL;
1639
1640   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
1641                                 &e->where);
1642
1643   rtrunc = gfc_copy_expr (e);
1644
1645   mpf_trunc (rtrunc->value.real, e->value.real);
1646   mpz_set_f (result->value.integer, rtrunc->value.real);
1647
1648   gfc_free_expr (rtrunc);
1649   return range_check (result, "IFIX");
1650 }
1651
1652
1653 gfc_expr *
1654 gfc_simplify_idint (gfc_expr * e)
1655 {
1656   gfc_expr *rtrunc, *result;
1657
1658   if (e->expr_type != EXPR_CONSTANT)
1659     return NULL;
1660
1661   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
1662                                 &e->where);
1663
1664   rtrunc = gfc_copy_expr (e);
1665
1666   mpf_trunc (rtrunc->value.real, e->value.real);
1667   mpz_set_f (result->value.integer, rtrunc->value.real);
1668
1669   gfc_free_expr (rtrunc);
1670   return range_check (result, "IDINT");
1671 }
1672
1673
1674 gfc_expr *
1675 gfc_simplify_ior (gfc_expr * x, gfc_expr * y)
1676 {
1677   gfc_expr *result;
1678
1679   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1680     return NULL;
1681
1682   result = gfc_constant_result (BT_INTEGER, x->ts.kind, &x->where);
1683
1684   mpz_ior (result->value.integer, x->value.integer, y->value.integer);
1685   return range_check (result, "IOR");
1686 }
1687
1688
1689 gfc_expr *
1690 gfc_simplify_ishft (gfc_expr * e, gfc_expr * s)
1691 {
1692   gfc_expr *result;
1693   int shift, ashift, isize, k;
1694   long e_int;
1695
1696   if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
1697     return NULL;
1698
1699   if (gfc_extract_int (s, &shift) != NULL)
1700     {
1701       gfc_error ("Invalid second argument of ISHFT at %L", &s->where);
1702       return &gfc_bad_expr;
1703     }
1704
1705   k = gfc_validate_kind (BT_INTEGER, e->ts.kind);
1706   if (k == -1)
1707     gfc_internal_error ("gfc_simplify_ishft(): Bad kind");
1708
1709   isize = gfc_integer_kinds[k].bit_size;
1710
1711   if (shift >= 0)
1712     ashift = shift;
1713   else
1714     ashift = -shift;
1715
1716   if (ashift > isize)
1717     {
1718       gfc_error
1719         ("Magnitude of second argument of ISHFT exceeds bit size at %L",
1720          &s->where);
1721       return &gfc_bad_expr;
1722     }
1723
1724   e_int = mpz_get_si (e->value.integer);
1725   if (e_int > INT_MAX || e_int < INT_MIN)
1726     gfc_internal_error ("ISHFT: unable to extract integer");
1727
1728   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
1729
1730   if (shift == 0)
1731     {
1732       mpz_set (result->value.integer, e->value.integer);
1733       return range_check (result, "ISHFT");
1734     }
1735
1736   if (shift > 0)
1737     mpz_set_si (result->value.integer, e_int << shift);
1738   else
1739     mpz_set_si (result->value.integer, e_int >> ashift);
1740
1741   return range_check (result, "ISHFT");
1742 }
1743
1744
1745 gfc_expr *
1746 gfc_simplify_ishftc (gfc_expr * e, gfc_expr * s, gfc_expr * sz)
1747 {
1748   gfc_expr *result;
1749   int shift, ashift, isize, delta, k;
1750   int i, *bits;
1751
1752   if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
1753     return NULL;
1754
1755   if (gfc_extract_int (s, &shift) != NULL)
1756     {
1757       gfc_error ("Invalid second argument of ISHFTC at %L", &s->where);
1758       return &gfc_bad_expr;
1759     }
1760
1761   k = gfc_validate_kind (e->ts.type, e->ts.kind);
1762   if (k == -1)
1763     gfc_internal_error ("gfc_simplify_ishftc(): Bad kind");
1764
1765   if (sz != NULL)
1766     {
1767       if (gfc_extract_int (sz, &isize) != NULL || isize < 0)
1768         {
1769           gfc_error ("Invalid third argument of ISHFTC at %L", &sz->where);
1770           return &gfc_bad_expr;
1771         }
1772     }
1773   else
1774     isize = gfc_integer_kinds[k].bit_size;
1775
1776   if (shift >= 0)
1777     ashift = shift;
1778   else
1779     ashift = -shift;
1780
1781   if (ashift > isize)
1782     {
1783       gfc_error
1784         ("Magnitude of second argument of ISHFTC exceeds third argument "
1785          "at %L", &s->where);
1786       return &gfc_bad_expr;
1787     }
1788
1789   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
1790
1791   bits = gfc_getmem (isize * sizeof (int));
1792
1793   for (i = 0; i < isize; i++)
1794     bits[i] = mpz_tstbit (e->value.integer, i);
1795
1796   delta = isize - ashift;
1797
1798   if (shift == 0)
1799     {
1800       mpz_set (result->value.integer, e->value.integer);
1801       gfc_free (bits);
1802       return range_check (result, "ISHFTC");
1803     }
1804
1805   else if (shift > 0)
1806     {
1807       for (i = 0; i < delta; i++)
1808         {
1809           if (bits[i] == 0)
1810             mpz_clrbit (result->value.integer, i + shift);
1811           if (bits[i] == 1)
1812             mpz_setbit (result->value.integer, i + shift);
1813         }
1814
1815       for (i = delta; i < isize; i++)
1816         {
1817           if (bits[i] == 0)
1818             mpz_clrbit (result->value.integer, i - delta);
1819           if (bits[i] == 1)
1820             mpz_setbit (result->value.integer, i - delta);
1821         }
1822
1823       gfc_free (bits);
1824       return range_check (result, "ISHFTC");
1825     }
1826   else
1827     {
1828       for (i = 0; i < ashift; i++)
1829         {
1830           if (bits[i] == 0)
1831             mpz_clrbit (result->value.integer, i + delta);
1832           if (bits[i] == 1)
1833             mpz_setbit (result->value.integer, i + delta);
1834         }
1835
1836       for (i = ashift; i < isize; i++)
1837         {
1838           if (bits[i] == 0)
1839             mpz_clrbit (result->value.integer, i + shift);
1840           if (bits[i] == 1)
1841             mpz_setbit (result->value.integer, i + shift);
1842         }
1843
1844       gfc_free (bits);
1845       return range_check (result, "ISHFTC");
1846     }
1847 }
1848
1849
1850 gfc_expr *
1851 gfc_simplify_kind (gfc_expr * e)
1852 {
1853
1854   if (e->ts.type == BT_DERIVED)
1855     {
1856       gfc_error ("Argument of KIND at %L is a DERIVED type", &e->where);
1857       return &gfc_bad_expr;
1858     }
1859
1860   return gfc_int_expr (e->ts.kind);
1861 }
1862
1863
1864 static gfc_expr *
1865 gfc_simplify_bound (gfc_expr * array, gfc_expr * dim, int upper)
1866 {
1867   gfc_ref *ref;
1868   gfc_array_spec *as;
1869   int i;
1870
1871   if (array->expr_type != EXPR_VARIABLE)
1872       return NULL;
1873
1874   if (dim == NULL)
1875     return NULL;
1876
1877   if (dim->expr_type != EXPR_CONSTANT)
1878     return NULL;
1879
1880   /* Follow any component references.  */
1881   as = array->symtree->n.sym->as;
1882   ref = array->ref;
1883   while (ref->next != NULL)
1884     {
1885       if (ref->type == REF_COMPONENT)
1886         as = ref->u.c.sym->as;
1887       ref = ref->next;
1888     }
1889
1890   if (ref->type != REF_ARRAY || ref->u.ar.type != AR_FULL)
1891     return NULL;
1892   
1893   i = mpz_get_si (dim->value.integer);
1894   if (upper) 
1895     return as->upper[i-1];
1896   else
1897     return as->lower[i-1];
1898 }
1899
1900
1901 gfc_expr *
1902 gfc_simplify_lbound (gfc_expr * array, gfc_expr * dim)
1903 {
1904   return gfc_simplify_bound (array, dim, 0);
1905 }
1906
1907
1908 gfc_expr *
1909 gfc_simplify_len (gfc_expr * e)
1910 {
1911   gfc_expr *result;
1912
1913   if (e->expr_type != EXPR_CONSTANT)
1914     return NULL;
1915
1916   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
1917                                 &e->where);
1918
1919   mpz_set_si (result->value.integer, e->value.character.length);
1920   return range_check (result, "LEN");
1921 }
1922
1923
1924 gfc_expr *
1925 gfc_simplify_len_trim (gfc_expr * e)
1926 {
1927   gfc_expr *result;
1928   int count, len, lentrim, i;
1929
1930   if (e->expr_type != EXPR_CONSTANT)
1931     return NULL;
1932
1933   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
1934                                 &e->where);
1935
1936   len = e->value.character.length;
1937
1938   for (count = 0, i = 1; i <= len; i++)
1939     if (e->value.character.string[len - i] == ' ')
1940       count++;
1941     else
1942       break;
1943
1944   lentrim = len - count;
1945
1946   mpz_set_si (result->value.integer, lentrim);
1947   return range_check (result, "LEN_TRIM");
1948 }
1949
1950
1951 gfc_expr *
1952 gfc_simplify_lge (gfc_expr * a, gfc_expr * b)
1953 {
1954
1955   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
1956     return NULL;
1957
1958   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) >= 0,
1959                            &a->where);
1960 }
1961
1962
1963 gfc_expr *
1964 gfc_simplify_lgt (gfc_expr * a, gfc_expr * b)
1965 {
1966
1967   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
1968     return NULL;
1969
1970   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) > 0,
1971                            &a->where);
1972 }
1973
1974
1975 gfc_expr *
1976 gfc_simplify_lle (gfc_expr * a, gfc_expr * b)
1977 {
1978
1979   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
1980     return NULL;
1981
1982   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) <= 0,
1983                            &a->where);
1984 }
1985
1986
1987 gfc_expr *
1988 gfc_simplify_llt (gfc_expr * a, gfc_expr * b)
1989 {
1990
1991   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
1992     return NULL;
1993
1994   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) < 0,
1995                            &a->where);
1996 }
1997
1998
1999 gfc_expr *
2000 gfc_simplify_log (gfc_expr * x)
2001 {
2002   gfc_expr *result;
2003   mpf_t xr, xi;
2004
2005   if (x->expr_type != EXPR_CONSTANT)
2006     return NULL;
2007
2008   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
2009
2010   switch (x->ts.type)
2011     {
2012     case BT_REAL:
2013       if (mpf_cmp (x->value.real, mpf_zero) <= 0)
2014         {
2015           gfc_error
2016             ("Argument of LOG at %L cannot be less than or equal to zero",
2017              &x->where);
2018           gfc_free_expr (result);
2019           return &gfc_bad_expr;
2020         }
2021
2022       natural_logarithm (&x->value.real, &result->value.real);
2023       break;
2024
2025     case BT_COMPLEX:
2026       if ((mpf_cmp (x->value.complex.r, mpf_zero) == 0)
2027           && (mpf_cmp (x->value.complex.i, mpf_zero) == 0))
2028         {
2029           gfc_error ("Complex argument of LOG at %L cannot be zero",
2030                      &x->where);
2031           gfc_free_expr (result);
2032           return &gfc_bad_expr;
2033         }
2034
2035       mpf_init (xr);
2036       mpf_init (xi);
2037
2038       arctangent2 (&x->value.complex.i, &x->value.complex.r,
2039         &result->value.complex.i);
2040
2041       mpf_mul (xr, x->value.complex.r, x->value.complex.r);
2042       mpf_mul (xi, x->value.complex.i, x->value.complex.i);
2043       mpf_add (xr, xr, xi);
2044       mpf_sqrt (xr, xr);
2045       natural_logarithm (&xr, &result->value.complex.r);
2046
2047       mpf_clear (xr);
2048       mpf_clear (xi);
2049
2050       break;
2051
2052     default:
2053       gfc_internal_error ("gfc_simplify_log: bad type");
2054     }
2055
2056   return range_check (result, "LOG");
2057 }
2058
2059
2060 gfc_expr *
2061 gfc_simplify_log10 (gfc_expr * x)
2062 {
2063   gfc_expr *result;
2064
2065   if (x->expr_type != EXPR_CONSTANT)
2066     return NULL;
2067
2068   if (mpf_cmp (x->value.real, mpf_zero) <= 0)
2069     {
2070       gfc_error
2071         ("Argument of LOG10 at %L cannot be less than or equal to zero",
2072          &x->where);
2073       return &gfc_bad_expr;
2074     }
2075
2076   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
2077
2078   common_logarithm (&x->value.real, &result->value.real);
2079
2080   return range_check (result, "LOG10");
2081 }
2082
2083
2084 gfc_expr *
2085 gfc_simplify_logical (gfc_expr * e, gfc_expr * k)
2086 {
2087   gfc_expr *result;
2088   int kind;
2089
2090   kind = get_kind (BT_LOGICAL, k, "LOGICAL", gfc_default_logical_kind ());
2091   if (kind < 0)
2092     return &gfc_bad_expr;
2093
2094   if (e->expr_type != EXPR_CONSTANT)
2095     return NULL;
2096
2097   result = gfc_constant_result (BT_LOGICAL, kind, &e->where);
2098
2099   result->value.logical = e->value.logical;
2100
2101   return result;
2102 }
2103
2104
2105 /* This function is special since MAX() can take any number of
2106    arguments.  The simplified expression is a rewritten version of the
2107    argument list containing at most one constant element.  Other
2108    constant elements are deleted.  Because the argument list has
2109    already been checked, this function always succeeds.  sign is 1 for
2110    MAX(), -1 for MIN().  */
2111
2112 static gfc_expr *
2113 simplify_min_max (gfc_expr * expr, int sign)
2114 {
2115   gfc_actual_arglist *arg, *last, *extremum;
2116   gfc_intrinsic_sym * specific;
2117
2118   last = NULL;
2119   extremum = NULL;
2120   specific = expr->value.function.isym;
2121
2122   arg = expr->value.function.actual;
2123
2124   for (; arg; last = arg, arg = arg->next)
2125     {
2126       if (arg->expr->expr_type != EXPR_CONSTANT)
2127         continue;
2128
2129       if (extremum == NULL)
2130         {
2131           extremum = arg;
2132           continue;
2133         }
2134
2135       switch (arg->expr->ts.type)
2136         {
2137         case BT_INTEGER:
2138           if (mpz_cmp (arg->expr->value.integer,
2139                        extremum->expr->value.integer) * sign > 0)
2140             mpz_set (extremum->expr->value.integer, arg->expr->value.integer);
2141
2142           break;
2143
2144         case BT_REAL:
2145           if (mpf_cmp (arg->expr->value.real, extremum->expr->value.real) *
2146               sign > 0)
2147             mpf_set (extremum->expr->value.real, arg->expr->value.real);
2148
2149           break;
2150
2151         default:
2152           gfc_internal_error ("gfc_simplify_max(): Bad type in arglist");
2153         }
2154
2155       /* Delete the extra constant argument.  */
2156       if (last == NULL)
2157         expr->value.function.actual = arg->next;
2158       else
2159         last->next = arg->next;
2160
2161       arg->next = NULL;
2162       gfc_free_actual_arglist (arg);
2163       arg = last;
2164     }
2165
2166   /* If there is one value left, replace the function call with the
2167      expression.  */
2168   if (expr->value.function.actual->next != NULL)
2169     return NULL;
2170
2171   /* Convert to the correct type and kind.  */
2172   if (expr->ts.type != BT_UNKNOWN) 
2173     return gfc_convert_constant (expr->value.function.actual->expr,
2174         expr->ts.type, expr->ts.kind);
2175
2176   if (specific->ts.type != BT_UNKNOWN) 
2177     return gfc_convert_constant (expr->value.function.actual->expr,
2178         specific->ts.type, specific->ts.kind); 
2179  
2180   return gfc_copy_expr (expr->value.function.actual->expr);
2181 }
2182
2183
2184 gfc_expr *
2185 gfc_simplify_min (gfc_expr * e)
2186 {
2187
2188   return simplify_min_max (e, -1);
2189 }
2190
2191
2192 gfc_expr *
2193 gfc_simplify_max (gfc_expr * e)
2194 {
2195
2196   return simplify_min_max (e, 1);
2197 }
2198
2199
2200 gfc_expr *
2201 gfc_simplify_maxexponent (gfc_expr * x)
2202 {
2203   gfc_expr *result;
2204   int i;
2205
2206   i = gfc_validate_kind (BT_REAL, x->ts.kind);
2207   if (i == -1)
2208     gfc_internal_error ("gfc_simplify_maxexponent(): Bad kind");
2209
2210   result = gfc_int_expr (gfc_real_kinds[i].max_exponent);
2211   result->where = x->where;
2212
2213   return result;
2214 }
2215
2216
2217 gfc_expr *
2218 gfc_simplify_minexponent (gfc_expr * x)
2219 {
2220   gfc_expr *result;
2221   int i;
2222
2223   i = gfc_validate_kind (BT_REAL, x->ts.kind);
2224   if (i == -1)
2225     gfc_internal_error ("gfc_simplify_minexponent(): Bad kind");
2226
2227   result = gfc_int_expr (gfc_real_kinds[i].min_exponent);
2228   result->where = x->where;
2229
2230   return result;
2231 }
2232
2233
2234 gfc_expr *
2235 gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
2236 {
2237   gfc_expr *result;
2238   mpf_t quot, iquot, term;
2239
2240   if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
2241     return NULL;
2242
2243   result = gfc_constant_result (a->ts.type, a->ts.kind, &a->where);
2244
2245   switch (a->ts.type)
2246     {
2247     case BT_INTEGER:
2248       if (mpz_cmp_ui (p->value.integer, 0) == 0)
2249         {
2250           /* Result is processor-dependent.  */
2251           gfc_error ("Second argument MOD at %L is zero", &a->where);
2252           gfc_free_expr (result);
2253           return &gfc_bad_expr;
2254         }
2255       mpz_tdiv_r (result->value.integer, a->value.integer, p->value.integer);
2256       break;
2257
2258     case BT_REAL:
2259       if (mpf_cmp_ui (p->value.real, 0) == 0)
2260         {
2261           /* Result is processor-dependent.  */
2262           gfc_error ("Second argument of MOD at %L is zero", &p->where);
2263           gfc_free_expr (result);
2264           return &gfc_bad_expr;
2265         }
2266
2267       mpf_init (quot);
2268       mpf_init (iquot);
2269       mpf_init (term);
2270
2271       mpf_div (quot, a->value.real, p->value.real);
2272       mpf_trunc (iquot, quot);
2273       mpf_mul (term, iquot, p->value.real);
2274       mpf_sub (result->value.real, a->value.real, term);
2275
2276       mpf_clear (quot);
2277       mpf_clear (iquot);
2278       mpf_clear (term);
2279       break;
2280
2281     default:
2282       gfc_internal_error ("gfc_simplify_mod(): Bad arguments");
2283     }
2284
2285   return range_check (result, "MOD");
2286 }
2287
2288
2289 gfc_expr *
2290 gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
2291 {
2292   gfc_expr *result;
2293   mpf_t quot, iquot, term;
2294
2295   if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
2296     return NULL;
2297
2298   result = gfc_constant_result (a->ts.type, a->ts.kind, &a->where);
2299
2300   switch (a->ts.type)
2301     {
2302     case BT_INTEGER:
2303       if (mpz_cmp_ui (p->value.integer, 0) == 0)
2304         {
2305           /* Result is processor-dependent. This processor just opts
2306              to not handle it at all.  */
2307           gfc_error ("Second argument of MODULO at %L is zero", &a->where);
2308           gfc_free_expr (result);
2309           return &gfc_bad_expr;
2310         }
2311       mpz_fdiv_r (result->value.integer, a->value.integer, p->value.integer);
2312
2313       break;
2314
2315     case BT_REAL:
2316       if (mpf_cmp_ui (p->value.real, 0) == 0)
2317         {
2318           /* Result is processor-dependent.  */
2319           gfc_error ("Second argument of MODULO at %L is zero", &p->where);
2320           gfc_free_expr (result);
2321           return &gfc_bad_expr;
2322         }
2323
2324       mpf_init (quot);
2325       mpf_init (iquot);
2326       mpf_init (term);
2327
2328       mpf_div (quot, a->value.real, p->value.real);
2329       mpf_floor (iquot, quot);
2330       mpf_mul (term, iquot, p->value.real);
2331
2332       mpf_clear (quot);
2333       mpf_clear (iquot);
2334       mpf_clear (term);
2335
2336       mpf_sub (result->value.real, a->value.real, term);
2337       break;
2338
2339     default:
2340       gfc_internal_error ("gfc_simplify_modulo(): Bad arguments");
2341     }
2342
2343   return range_check (result, "MODULO");
2344 }
2345
2346
2347 /* Exists for the sole purpose of consistency with other intrinsics.  */
2348 gfc_expr *
2349 gfc_simplify_mvbits (gfc_expr * f  ATTRIBUTE_UNUSED,
2350                      gfc_expr * fp ATTRIBUTE_UNUSED,
2351                      gfc_expr * l  ATTRIBUTE_UNUSED,
2352                      gfc_expr * to ATTRIBUTE_UNUSED,
2353                      gfc_expr * tp ATTRIBUTE_UNUSED)
2354 {
2355   return NULL;
2356 }
2357
2358
2359 gfc_expr *
2360 gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
2361 {
2362   gfc_expr *result;
2363   float rval;
2364   double val, eps;
2365   int p, i, k, match_float;
2366
2367   /* FIXME: This implementation is dopey and probably not quite right,
2368      but it's a start.  */
2369
2370   if (x->expr_type != EXPR_CONSTANT)
2371     return NULL;
2372
2373   k = gfc_validate_kind (x->ts.type, x->ts.kind);
2374   if (k == -1)
2375     gfc_internal_error ("gfc_simplify_precision(): Bad kind");
2376
2377   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
2378
2379   val = mpf_get_d (x->value.real);
2380   p = gfc_real_kinds[k].digits;
2381
2382   eps = 1.;
2383   for (i = 1; i < p; ++i)
2384     {
2385       eps = eps / 2.;
2386     }
2387
2388   /* TODO we should make sure that 'float' matches kind 4 */
2389   match_float = gfc_real_kinds[k].kind == 4;
2390   if (mpf_cmp_ui (s->value.real, 0) > 0)
2391     {
2392       if (match_float)
2393         {
2394           rval = (float) val;
2395           rval = rval + eps;
2396           mpf_set_d (result->value.real, rval);
2397         }
2398       else
2399         {
2400           val = val + eps;
2401           mpf_set_d (result->value.real, val);
2402         }
2403     }
2404   else if (mpf_cmp_ui (s->value.real, 0) < 0)
2405     {
2406       if (match_float)
2407         {
2408           rval = (float) val;
2409           rval = rval - eps;
2410           mpf_set_d (result->value.real, rval);
2411         }
2412       else
2413         {
2414           val = val - eps;
2415           mpf_set_d (result->value.real, val);
2416         }
2417     }
2418   else
2419     {
2420       gfc_error ("Invalid second argument of NEAREST at %L", &s->where);
2421       gfc_free (result);
2422       return &gfc_bad_expr;
2423     }
2424
2425   return range_check (result, "NEAREST");
2426
2427 }
2428
2429
2430 static gfc_expr *
2431 simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
2432 {
2433   gfc_expr *rtrunc, *itrunc, *result;
2434   int kind, cmp;
2435
2436   kind = get_kind (BT_INTEGER, k, name, gfc_default_integer_kind ());
2437   if (kind == -1)
2438     return &gfc_bad_expr;
2439
2440   if (e->expr_type != EXPR_CONSTANT)
2441     return NULL;
2442
2443   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
2444
2445   rtrunc = gfc_copy_expr (e);
2446   itrunc = gfc_copy_expr (e);
2447
2448   cmp = mpf_cmp_ui (e->value.real, 0);
2449
2450   if (cmp > 0)
2451     {
2452       mpf_add (rtrunc->value.real, e->value.real, mpf_half);
2453       mpf_trunc (itrunc->value.real, rtrunc->value.real);
2454     }
2455   else if (cmp < 0)
2456     {
2457       mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
2458       mpf_trunc (itrunc->value.real, rtrunc->value.real);
2459     }
2460   else
2461     mpf_set_ui (itrunc->value.real, 0);
2462
2463   mpz_set_f (result->value.integer, itrunc->value.real);
2464
2465   gfc_free_expr (itrunc);
2466   gfc_free_expr (rtrunc);
2467
2468   return range_check (result, name);
2469 }
2470
2471
2472 gfc_expr *
2473 gfc_simplify_nint (gfc_expr * e, gfc_expr * k)
2474 {
2475
2476   return simplify_nint ("NINT", e, k);
2477 }
2478
2479
2480 gfc_expr *
2481 gfc_simplify_idnint (gfc_expr * e)
2482 {
2483
2484   return simplify_nint ("IDNINT", e, NULL);
2485 }
2486
2487
2488 gfc_expr *
2489 gfc_simplify_not (gfc_expr * e)
2490 {
2491   gfc_expr *result;
2492   int i;
2493
2494   if (e->expr_type != EXPR_CONSTANT)
2495     return NULL;
2496
2497   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
2498
2499   mpz_com (result->value.integer, e->value.integer);
2500
2501   /* Because of how GMP handles numbers, the result must be ANDed with
2502      the max_int mask.  For radices <> 2, this will require change.  */
2503
2504   i = gfc_validate_kind (BT_INTEGER, e->ts.kind);
2505   if (i == -1)
2506     gfc_internal_error ("gfc_simplify_not(): Bad kind");
2507
2508   mpz_and (result->value.integer, result->value.integer,
2509            gfc_integer_kinds[i].max_int);
2510
2511   return range_check (result, "NOT");
2512 }
2513
2514
2515 gfc_expr *
2516 gfc_simplify_null (gfc_expr * mold)
2517 {
2518   gfc_expr *result;
2519
2520   result = gfc_get_expr ();
2521   result->expr_type = EXPR_NULL;
2522
2523   if (mold == NULL)
2524     result->ts.type = BT_UNKNOWN;
2525   else
2526     {
2527       result->ts = mold->ts;
2528       result->where = mold->where;
2529     }
2530
2531   return result;
2532 }
2533
2534
2535 gfc_expr *
2536 gfc_simplify_precision (gfc_expr * e)
2537 {
2538   gfc_expr *result;
2539   int i;
2540
2541   i = gfc_validate_kind (e->ts.type, e->ts.kind);
2542   if (i == -1)
2543     gfc_internal_error ("gfc_simplify_precision(): Bad kind");
2544
2545   result = gfc_int_expr (gfc_real_kinds[i].precision);
2546   result->where = e->where;
2547
2548   return result;
2549 }
2550
2551
2552 gfc_expr *
2553 gfc_simplify_radix (gfc_expr * e)
2554 {
2555   gfc_expr *result;
2556   int i;
2557
2558   i = gfc_validate_kind (e->ts.type, e->ts.kind);
2559   if (i == -1)
2560     goto bad;
2561
2562   switch (e->ts.type)
2563     {
2564     case BT_INTEGER:
2565       i = gfc_integer_kinds[i].radix;
2566       break;
2567
2568     case BT_REAL:
2569       i = gfc_real_kinds[i].radix;
2570       break;
2571
2572     default:
2573     bad:
2574       gfc_internal_error ("gfc_simplify_radix(): Bad type");
2575     }
2576
2577   result = gfc_int_expr (i);
2578   result->where = e->where;
2579
2580   return result;
2581 }
2582
2583
2584 gfc_expr *
2585 gfc_simplify_range (gfc_expr * e)
2586 {
2587   gfc_expr *result;
2588   int i;
2589   long j;
2590
2591   i = gfc_validate_kind (e->ts.type, e->ts.kind);
2592   if (i == -1)
2593     goto bad_type;
2594
2595   switch (e->ts.type)
2596     {
2597     case BT_INTEGER:
2598       j = gfc_integer_kinds[i].range;
2599       break;
2600
2601     case BT_REAL:
2602     case BT_COMPLEX:
2603       j = gfc_real_kinds[i].range;
2604       break;
2605
2606     bad_type:
2607     default:
2608       gfc_internal_error ("gfc_simplify_range(): Bad kind");
2609     }
2610
2611   result = gfc_int_expr (j);
2612   result->where = e->where;
2613
2614   return result;
2615 }
2616
2617
2618 gfc_expr *
2619 gfc_simplify_real (gfc_expr * e, gfc_expr * k)
2620 {
2621   gfc_expr *result;
2622   int kind;
2623
2624   if (e->ts.type == BT_COMPLEX)
2625     kind = get_kind (BT_REAL, k, "REAL", e->ts.kind);
2626   else
2627     kind = get_kind (BT_REAL, k, "REAL", gfc_default_real_kind ());
2628
2629   if (kind == -1)
2630     return &gfc_bad_expr;
2631
2632   if (e->expr_type != EXPR_CONSTANT)
2633     return NULL;
2634
2635   switch (e->ts.type)
2636     {
2637     case BT_INTEGER:
2638       result = gfc_int2real (e, kind);
2639       break;
2640
2641     case BT_REAL:
2642       result = gfc_real2real (e, kind);
2643       break;
2644
2645     case BT_COMPLEX:
2646       result = gfc_complex2real (e, kind);
2647       break;
2648
2649     default:
2650       gfc_internal_error ("bad type in REAL");
2651       /* Not reached */
2652     }
2653
2654   return range_check (result, "REAL");
2655 }
2656
2657 gfc_expr *
2658 gfc_simplify_repeat (gfc_expr * e, gfc_expr * n)
2659 {
2660   gfc_expr *result;
2661   int i, j, len, ncopies, nlen;
2662
2663   if (e->expr_type != EXPR_CONSTANT || n->expr_type != EXPR_CONSTANT)
2664     return NULL;
2665
2666   if (n != NULL && (gfc_extract_int (n, &ncopies) != NULL || ncopies < 0))
2667     {
2668       gfc_error ("Invalid second argument of REPEAT at %L", &n->where);
2669       return &gfc_bad_expr;
2670     }
2671
2672   len = e->value.character.length;
2673   nlen = ncopies * len;
2674
2675   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
2676
2677   if (ncopies == 0)
2678     {
2679       result->value.character.string = gfc_getmem (1);
2680       result->value.character.length = 0;
2681       result->value.character.string[0] = '\0';
2682       return result;
2683     }
2684
2685   result->value.character.length = nlen;
2686   result->value.character.string = gfc_getmem (nlen + 1);
2687
2688   for (i = 0; i < ncopies; i++)
2689     for (j = 0; j < len; j++)
2690       result->value.character.string[j + i * len] =
2691         e->value.character.string[j];
2692
2693   result->value.character.string[nlen] = '\0';  /* For debugger */
2694   return result;
2695 }
2696
2697
2698 /* This one is a bear, but mainly has to do with shuffling elements.  */
2699
2700 gfc_expr *
2701 gfc_simplify_reshape (gfc_expr * source, gfc_expr * shape_exp,
2702                       gfc_expr * pad, gfc_expr * order_exp)
2703 {
2704
2705   int order[GFC_MAX_DIMENSIONS], shape[GFC_MAX_DIMENSIONS];
2706   int i, rank, npad, x[GFC_MAX_DIMENSIONS];
2707   gfc_constructor *head, *tail;
2708   mpz_t index, size;
2709   unsigned long j;
2710   size_t nsource;
2711   gfc_expr *e;
2712
2713   /* Unpack the shape array.  */
2714   if (source->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (source))
2715     return NULL;
2716
2717   if (shape_exp->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (shape_exp))
2718     return NULL;
2719
2720   if (pad != NULL
2721       && (pad->expr_type != EXPR_ARRAY
2722           || !gfc_is_constant_expr (pad)))
2723     return NULL;
2724
2725   if (order_exp != NULL
2726       && (order_exp->expr_type != EXPR_ARRAY
2727           || !gfc_is_constant_expr (order_exp)))
2728     return NULL;
2729
2730   mpz_init (index);
2731   rank = 0;
2732   head = tail = NULL;
2733
2734   for (;;)
2735     {
2736       e = gfc_get_array_element (shape_exp, rank);
2737       if (e == NULL)
2738         break;
2739
2740       if (gfc_extract_int (e, &shape[rank]) != NULL)
2741         {
2742           gfc_error ("Integer too large in shape specification at %L",
2743                      &e->where);
2744           gfc_free_expr (e);
2745           goto bad_reshape;
2746         }
2747
2748       gfc_free_expr (e);
2749
2750       if (rank >= GFC_MAX_DIMENSIONS)
2751         {
2752           gfc_error ("Too many dimensions in shape specification for RESHAPE "
2753                      "at %L", &e->where);
2754
2755           goto bad_reshape;
2756         }
2757
2758       if (shape[rank] < 0)
2759         {
2760           gfc_error ("Shape specification at %L cannot be negative",
2761                      &e->where);
2762           goto bad_reshape;
2763         }
2764
2765       rank++;
2766     }
2767
2768   if (rank == 0)
2769     {
2770       gfc_error ("Shape specification at %L cannot be the null array",
2771                  &shape_exp->where);
2772       goto bad_reshape;
2773     }
2774
2775   /* Now unpack the order array if present.  */
2776   if (order_exp == NULL)
2777     {
2778       for (i = 0; i < rank; i++)
2779         order[i] = i;
2780
2781     }
2782   else
2783     {
2784
2785       for (i = 0; i < rank; i++)
2786         x[i] = 0;
2787
2788       for (i = 0; i < rank; i++)
2789         {
2790           e = gfc_get_array_element (order_exp, i);
2791           if (e == NULL)
2792             {
2793               gfc_error
2794                 ("ORDER parameter of RESHAPE at %L is not the same size "
2795                  "as SHAPE parameter", &order_exp->where);
2796               goto bad_reshape;
2797             }
2798
2799           if (gfc_extract_int (e, &order[i]) != NULL)
2800             {
2801               gfc_error ("Error in ORDER parameter of RESHAPE at %L",
2802                          &e->where);
2803               gfc_free_expr (e);
2804               goto bad_reshape;
2805             }
2806
2807           gfc_free_expr (e);
2808
2809           if (order[i] < 1 || order[i] > rank)
2810             {
2811               gfc_error ("ORDER parameter of RESHAPE at %L is out of range",
2812                          &e->where);
2813               goto bad_reshape;
2814             }
2815
2816           order[i]--;
2817
2818           if (x[order[i]])
2819             {
2820               gfc_error ("Invalid permutation in ORDER parameter at %L",
2821                          &e->where);
2822               goto bad_reshape;
2823             }
2824
2825           x[order[i]] = 1;
2826         }
2827     }
2828
2829   /* Count the elements in the source and padding arrays.  */
2830
2831   npad = 0;
2832   if (pad != NULL)
2833     {
2834       gfc_array_size (pad, &size);
2835       npad = mpz_get_ui (size);
2836       mpz_clear (size);
2837     }
2838
2839   gfc_array_size (source, &size);
2840   nsource = mpz_get_ui (size);
2841   mpz_clear (size);
2842
2843   /* If it weren't for that pesky permutation we could just loop
2844      through the source and round out any shortage with pad elements.
2845      But no, someone just had to have the compiler do something the
2846      user should be doing.  */
2847
2848   for (i = 0; i < rank; i++)
2849     x[i] = 0;
2850
2851   for (;;)
2852     {
2853       /* Figure out which element to extract.  */
2854       mpz_set_ui (index, 0);
2855
2856       for (i = rank - 1; i >= 0; i--)
2857         {
2858           mpz_add_ui (index, index, x[order[i]]);
2859           if (i != 0)
2860             mpz_mul_ui (index, index, shape[order[i - 1]]);
2861         }
2862
2863       if (mpz_cmp_ui (index, INT_MAX) > 0)
2864         gfc_internal_error ("Reshaped array too large at %L", &e->where);
2865
2866       j = mpz_get_ui (index);
2867
2868       if (j < nsource)
2869         e = gfc_get_array_element (source, j);
2870       else
2871         {
2872           j = j - nsource;
2873
2874           if (npad == 0)
2875             {
2876               gfc_error
2877                 ("PAD parameter required for short SOURCE parameter at %L",
2878                  &source->where);
2879               goto bad_reshape;
2880             }
2881
2882           j = j % npad;
2883           e = gfc_get_array_element (pad, j);
2884         }
2885
2886       if (head == NULL)
2887         head = tail = gfc_get_constructor ();
2888       else
2889         {
2890           tail->next = gfc_get_constructor ();
2891           tail = tail->next;
2892         }
2893
2894       if (e == NULL)
2895         goto bad_reshape;
2896
2897       tail->where = e->where;
2898       tail->expr = e;
2899
2900       /* Calculate the next element.  */
2901       i = 0;
2902
2903 inc:
2904       if (++x[i] < shape[i])
2905         continue;
2906       x[i++] = 0;
2907       if (i < rank)
2908         goto inc;
2909
2910       break;
2911     }
2912
2913   mpz_clear (index);
2914
2915   e = gfc_get_expr ();
2916   e->where = source->where;
2917   e->expr_type = EXPR_ARRAY;
2918   e->value.constructor = head;
2919   e->shape = gfc_get_shape (rank);
2920
2921   for (i = 0; i < rank; i++)
2922     mpz_init_set_ui (e->shape[i], shape[order[i]]);
2923
2924   e->ts = head->expr->ts;
2925   e->rank = rank;
2926
2927   return e;
2928
2929 bad_reshape:
2930   gfc_free_constructor (head);
2931   mpz_clear (index);
2932   return &gfc_bad_expr;
2933 }
2934
2935
2936 gfc_expr *
2937 gfc_simplify_rrspacing (gfc_expr * x)
2938 {
2939   gfc_expr *result;
2940   mpf_t i2, absv, ln2, lnx, frac, pow2;
2941   unsigned long exp2;
2942   int i, p;
2943
2944   if (x->expr_type != EXPR_CONSTANT)
2945     return NULL;
2946
2947   i = gfc_validate_kind (x->ts.type, x->ts.kind);
2948   if (i == -1)
2949     gfc_internal_error ("gfc_simplify_rrspacing(): Bad kind");
2950
2951   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
2952
2953   p = gfc_real_kinds[i].digits;
2954
2955   if (mpf_cmp (x->value.real, mpf_zero) == 0)
2956     {
2957       mpf_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny);
2958       return result;
2959     }
2960
2961   mpf_init_set_ui (i2, 2);
2962   mpf_init (ln2);
2963   mpf_init (absv);
2964   mpf_init (lnx);
2965   mpf_init (frac);
2966   mpf_init (pow2);
2967
2968   natural_logarithm (&i2, &ln2);
2969
2970   mpf_abs (absv, x->value.real);
2971   natural_logarithm (&absv, &lnx);
2972
2973   mpf_div (lnx, lnx, ln2);
2974   mpf_trunc (lnx, lnx);
2975   mpf_add_ui (lnx, lnx, 1);
2976
2977   exp2 = (unsigned long) mpf_get_d (lnx);
2978   mpf_pow_ui (pow2, i2, exp2);
2979   mpf_div (frac, absv, pow2);
2980
2981   exp2 = (unsigned long) p;
2982   mpf_mul_2exp (result->value.real, frac, exp2);
2983
2984   mpf_clear (i2);
2985   mpf_clear (ln2);
2986   mpf_clear (absv);
2987   mpf_clear (lnx);
2988   mpf_clear (frac);
2989   mpf_clear (pow2);
2990
2991   return range_check (result, "RRSPACING");
2992 }
2993
2994
2995 gfc_expr *
2996 gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
2997 {
2998   int k, neg_flag, power, exp_range;
2999   mpf_t scale, radix;
3000   gfc_expr *result;
3001
3002   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
3003     return NULL;
3004
3005   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3006
3007   if (mpf_sgn (x->value.real) == 0)
3008     {
3009       mpf_set_ui (result->value.real, 0);
3010       return result;
3011     }
3012
3013   k = gfc_validate_kind (BT_REAL, x->ts.kind);
3014   if (k == -1)
3015     gfc_internal_error ("gfc_simplify_scale(): Bad kind");
3016
3017   exp_range = gfc_real_kinds[k].max_exponent - gfc_real_kinds[k].min_exponent;
3018
3019   /* This check filters out values of i that would overflow an int.  */
3020   if (mpz_cmp_si (i->value.integer, exp_range + 2) > 0
3021       || mpz_cmp_si (i->value.integer, -exp_range - 2) < 0)
3022     {
3023       gfc_error ("Result of SCALE overflows its kind at %L", &result->where);
3024       return &gfc_bad_expr;
3025     }
3026
3027   /* Compute scale = radix ** power.  */
3028   power = mpz_get_si (i->value.integer);
3029
3030   if (power >= 0)
3031     neg_flag = 0;
3032   else
3033     {
3034       neg_flag = 1;
3035       power = -power;
3036     }
3037
3038   mpf_init_set_ui (radix, gfc_real_kinds[k].radix);
3039   mpf_init (scale);
3040   mpf_pow_ui (scale, radix, power);
3041
3042   if (neg_flag)
3043     mpf_div (result->value.real, x->value.real, scale);
3044   else
3045     mpf_mul (result->value.real, x->value.real, scale);
3046
3047   mpf_clear (scale);
3048   mpf_clear (radix);
3049
3050   return range_check (result, "SCALE");
3051 }
3052
3053
3054 gfc_expr *
3055 gfc_simplify_scan (gfc_expr * e, gfc_expr * c, gfc_expr * b)
3056 {
3057   gfc_expr *result;
3058   int back;
3059   size_t i;
3060   size_t indx, len, lenc;
3061
3062   if (e->expr_type != EXPR_CONSTANT || c->expr_type != EXPR_CONSTANT)
3063     return NULL;
3064
3065   if (b != NULL && b->value.logical != 0)
3066     back = 1;
3067   else
3068     back = 0;
3069
3070   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
3071                                 &e->where);
3072
3073   len = e->value.character.length;
3074   lenc = c->value.character.length;
3075
3076   if (len == 0 || lenc == 0)
3077     {
3078       indx = 0;
3079     }
3080   else
3081     {
3082       if (back == 0)
3083         {
3084           indx =
3085             strcspn (e->value.character.string, c->value.character.string) + 1;
3086           if (indx > len)
3087             indx = 0;
3088         }
3089       else
3090         {
3091           i = 0;
3092           for (indx = len; indx > 0; indx--)
3093             {
3094               for (i = 0; i < lenc; i++)
3095                 {
3096                   if (c->value.character.string[i]
3097                         == e->value.character.string[indx - 1])
3098                     break;
3099                 }
3100               if (i < lenc)
3101                 break;
3102             }
3103         }
3104     }
3105   mpz_set_ui (result->value.integer, indx);
3106   return range_check (result, "SCAN");
3107 }
3108
3109
3110 gfc_expr *
3111 gfc_simplify_selected_int_kind (gfc_expr * e)
3112 {
3113   int i, kind, range;
3114   gfc_expr *result;
3115
3116   if (e->expr_type != EXPR_CONSTANT || gfc_extract_int (e, &range) != NULL)
3117     return NULL;
3118
3119   kind = INT_MAX;
3120
3121   for (i = 0; gfc_integer_kinds[i].kind != 0; i++)
3122     if (gfc_integer_kinds[i].range >= range
3123         && gfc_integer_kinds[i].kind < kind)
3124       kind = gfc_integer_kinds[i].kind;
3125
3126   if (kind == INT_MAX)
3127     kind = -1;
3128
3129   result = gfc_int_expr (kind);
3130   result->where = e->where;
3131
3132   return result;
3133 }
3134
3135
3136 gfc_expr *
3137 gfc_simplify_selected_real_kind (gfc_expr * p, gfc_expr * q)
3138 {
3139   int range, precision, i, kind, found_precision, found_range;
3140   gfc_expr *result;
3141
3142   if (p == NULL)
3143     precision = 0;
3144   else
3145     {
3146       if (p->expr_type != EXPR_CONSTANT
3147           || gfc_extract_int (p, &precision) != NULL)
3148         return NULL;
3149     }
3150
3151   if (q == NULL)
3152     range = 0;
3153   else
3154     {
3155       if (q->expr_type != EXPR_CONSTANT
3156           || gfc_extract_int (q, &range) != NULL)
3157         return NULL;
3158     }
3159
3160   kind = INT_MAX;
3161   found_precision = 0;
3162   found_range = 0;
3163
3164   for (i = 0; gfc_real_kinds[i].kind != 0; i++)
3165     {
3166       if (gfc_real_kinds[i].precision >= precision)
3167         found_precision = 1;
3168
3169       if (gfc_real_kinds[i].range >= range)
3170         found_range = 1;
3171
3172       if (gfc_real_kinds[i].precision >= precision
3173           && gfc_real_kinds[i].range >= range && gfc_real_kinds[i].kind < kind)
3174         kind = gfc_real_kinds[i].kind;
3175     }
3176
3177   if (kind == INT_MAX)
3178     {
3179       kind = 0;
3180
3181       if (!found_precision)
3182         kind = -1;
3183       if (!found_range)
3184         kind -= 2;
3185     }
3186
3187   result = gfc_int_expr (kind);
3188   result->where = (p != NULL) ? p->where : q->where;
3189
3190   return result;
3191 }
3192
3193
3194 gfc_expr *
3195 gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
3196 {
3197   gfc_expr *result;
3198   mpf_t i2, ln2, absv, lnx, pow2, frac;
3199   unsigned long exp2;
3200
3201   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
3202     return NULL;
3203
3204   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3205
3206   if (mpf_cmp (x->value.real, mpf_zero) == 0)
3207     {
3208       mpf_set (result->value.real, mpf_zero);
3209       return result;
3210     }
3211
3212   mpf_init_set_ui (i2, 2);
3213   mpf_init (ln2);
3214   mpf_init (absv);
3215   mpf_init (lnx);
3216   mpf_init (pow2);
3217   mpf_init (frac);
3218
3219   natural_logarithm (&i2, &ln2);
3220
3221   mpf_abs (absv, x->value.real);
3222   natural_logarithm (&absv, &lnx);
3223
3224   mpf_div (lnx, lnx, ln2);
3225   mpf_trunc (lnx, lnx);
3226   mpf_add_ui (lnx, lnx, 1);
3227
3228   /* Old exponent value, and fraction.  */
3229   exp2 = (unsigned long) mpf_get_d (lnx);
3230   mpf_pow_ui (pow2, i2, exp2);
3231
3232   mpf_div (frac, absv, pow2);
3233
3234   /* New exponent.  */
3235   exp2 = (unsigned long) mpz_get_d (i->value.integer);
3236   mpf_mul_2exp (result->value.real, frac, exp2);
3237
3238   mpf_clear (i2);
3239   mpf_clear (ln2);
3240   mpf_clear (absv);
3241   mpf_clear (lnx);
3242   mpf_clear (pow2);
3243   mpf_clear (frac);
3244
3245   return range_check (result, "SET_EXPONENT");
3246 }
3247
3248
3249 gfc_expr *
3250 gfc_simplify_shape (gfc_expr * source)
3251 {
3252   mpz_t shape[GFC_MAX_DIMENSIONS];
3253   gfc_expr *result, *e, *f;
3254   gfc_array_ref *ar;
3255   int n;
3256   try t;
3257
3258   result = gfc_start_constructor (BT_INTEGER, gfc_default_integer_kind (),
3259                                   &source->where);
3260
3261   if (source->rank == 0 || source->expr_type != EXPR_VARIABLE)
3262     return result;
3263
3264   ar = gfc_find_array_ref (source);
3265
3266   t = gfc_array_ref_shape (ar, shape);
3267
3268   for (n = 0; n < source->rank; n++)
3269     {
3270       e = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
3271                                &source->where);
3272
3273       if (t == SUCCESS)
3274         {
3275           mpz_set (e->value.integer, shape[n]);
3276           mpz_clear (shape[n]);
3277         }
3278       else
3279         {
3280           mpz_set_ui (e->value.integer, n + 1);
3281
3282           f = gfc_simplify_size (source, e);
3283           gfc_free_expr (e);
3284           if (f == NULL)
3285             {
3286               gfc_free_expr (result);
3287               return NULL;
3288             }
3289           else
3290             {
3291               e = f;
3292             }
3293         }
3294
3295       gfc_append_constructor (result, e);
3296     }
3297
3298   return result;
3299 }
3300
3301
3302 gfc_expr *
3303 gfc_simplify_size (gfc_expr * array, gfc_expr * dim)
3304 {
3305   mpz_t size;
3306   gfc_expr *result;
3307   int d;
3308
3309   if (dim == NULL)
3310     {
3311       if (gfc_array_size (array, &size) == FAILURE)
3312         return NULL;
3313     }
3314   else
3315     {
3316       if (dim->expr_type != EXPR_CONSTANT)
3317         return NULL;
3318
3319       d = mpz_get_ui (dim->value.integer) - 1;
3320       if (gfc_array_dimen_size (array, d, &size) == FAILURE)
3321         return NULL;
3322     }
3323
3324   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
3325                                 &array->where);
3326
3327   mpz_set (result->value.integer, size);
3328
3329   return result;
3330 }
3331
3332
3333 gfc_expr *
3334 gfc_simplify_sign (gfc_expr * x, gfc_expr * y)
3335 {
3336   gfc_expr *result;
3337
3338   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
3339     return NULL;
3340
3341   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3342
3343   switch (x->ts.type)
3344     {
3345     case BT_INTEGER:
3346       mpz_abs (result->value.integer, x->value.integer);
3347       if (mpz_sgn (y->value.integer) < 0)
3348         mpz_neg (result->value.integer, result->value.integer);
3349
3350       break;
3351
3352     case BT_REAL:
3353       /* TODO: Handle -0.0 and +0.0 correctly on machines that support
3354          it.  */
3355       mpf_abs (result->value.real, x->value.real);
3356       if (mpf_sgn (y->value.integer) < 0)
3357         mpf_neg (result->value.real, result->value.real);
3358
3359       break;
3360
3361     default:
3362       gfc_internal_error ("Bad type in gfc_simplify_sign");
3363     }
3364
3365   return result;
3366 }
3367
3368
3369 gfc_expr *
3370 gfc_simplify_sin (gfc_expr * x)
3371 {
3372   gfc_expr *result;
3373   mpf_t xp, xq;
3374
3375   if (x->expr_type != EXPR_CONSTANT)
3376     return NULL;
3377
3378   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3379
3380   switch (x->ts.type)
3381     {
3382     case BT_REAL:
3383       sine (&x->value.real, &result->value.real);
3384       break;
3385
3386     case BT_COMPLEX:
3387       mpf_init (xp);
3388       mpf_init (xq);
3389
3390       sine (&x->value.complex.r, &xp);
3391       hypercos (&x->value.complex.i, &xq);
3392       mpf_mul (result->value.complex.r, xp, xq);
3393
3394       cosine (&x->value.complex.r, &xp);
3395       hypersine (&x->value.complex.i, &xq);
3396       mpf_mul (result->value.complex.i, xp, xq);
3397
3398       mpf_clear (xp);
3399       mpf_clear (xq);
3400       break;
3401
3402     default:
3403       gfc_internal_error ("in gfc_simplify_sin(): Bad type");
3404     }
3405
3406   return range_check (result, "SIN");
3407 }
3408
3409
3410 gfc_expr *
3411 gfc_simplify_sinh (gfc_expr * x)
3412 {
3413   gfc_expr *result;
3414
3415   if (x->expr_type != EXPR_CONSTANT)
3416     return NULL;
3417
3418   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3419
3420   hypersine (&x->value.real, &result->value.real);
3421
3422   return range_check (result, "SINH");
3423 }
3424
3425
3426 /* The argument is always a double precision real that is converted to
3427    single precision.  TODO: Rounding!  */
3428
3429 gfc_expr *
3430 gfc_simplify_sngl (gfc_expr * a)
3431 {
3432   gfc_expr *result;
3433
3434   if (a->expr_type != EXPR_CONSTANT)
3435     return NULL;
3436
3437   result = gfc_real2real (a, gfc_default_real_kind ());
3438   return range_check (result, "SNGL");
3439 }
3440
3441
3442 gfc_expr *
3443 gfc_simplify_spacing (gfc_expr * x)
3444 {
3445   gfc_expr *result;
3446   mpf_t i1, i2, ln2, absv, lnx;
3447   long diff;
3448   unsigned long exp2;
3449   int i, p;
3450
3451   if (x->expr_type != EXPR_CONSTANT)
3452     return NULL;
3453
3454   i = gfc_validate_kind (x->ts.type, x->ts.kind);
3455   if (i == -1)
3456     gfc_internal_error ("gfc_simplify_spacing(): Bad kind");
3457
3458   p = gfc_real_kinds[i].digits;
3459
3460   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3461
3462   if (mpf_cmp (x->value.real, mpf_zero) == 0)
3463     {
3464       mpf_set (result->value.real, gfc_real_kinds[i].tiny);
3465       return result;
3466     }
3467
3468   mpf_init_set_ui (i1, 1);
3469   mpf_init_set_ui (i2, 2);
3470   mpf_init (ln2);
3471   mpf_init (absv);
3472   mpf_init (lnx);
3473
3474   natural_logarithm (&i2, &ln2);
3475
3476   mpf_abs (absv, x->value.real);
3477   natural_logarithm (&absv, &lnx);
3478
3479   mpf_div (lnx, lnx, ln2);
3480   mpf_trunc (lnx, lnx);
3481   mpf_add_ui (lnx, lnx, 1);
3482
3483   diff = (long) mpf_get_d (lnx) - (long) p;
3484   if (diff >= 0)
3485     {
3486       exp2 = (unsigned) diff;
3487       mpf_mul_2exp (result->value.real, i1, exp2);
3488     }
3489   else
3490     {
3491       diff = -diff;
3492       exp2 = (unsigned) diff;
3493       mpf_div_2exp (result->value.real, i1, exp2);
3494     }
3495
3496   mpf_clear (i1);
3497   mpf_clear (i2);
3498   mpf_clear (ln2);
3499   mpf_clear (absv);
3500   mpf_clear (lnx);
3501
3502   if (mpf_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
3503     mpf_set (result->value.real, gfc_real_kinds[i].tiny);
3504
3505   return range_check (result, "SPACING");
3506 }
3507
3508
3509 gfc_expr *
3510 gfc_simplify_sqrt (gfc_expr * e)
3511 {
3512   gfc_expr *result;
3513   mpf_t ac, ad, s, t, w;
3514
3515   if (e->expr_type != EXPR_CONSTANT)
3516     return NULL;
3517
3518   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
3519
3520   switch (e->ts.type)
3521     {
3522     case BT_REAL:
3523       if (mpf_cmp_si (e->value.real, 0) < 0)
3524         goto negative_arg;
3525       mpf_sqrt (result->value.real, e->value.real);
3526
3527       break;
3528
3529     case BT_COMPLEX:
3530       /* Formula taken from Numerical Recipes to avoid over- and
3531          underflow.  */
3532
3533       mpf_init (ac);
3534       mpf_init (ad);
3535       mpf_init (s);
3536       mpf_init (t);
3537       mpf_init (w);
3538
3539       if (mpf_cmp_ui (e->value.complex.r, 0) == 0
3540           && mpf_cmp_ui (e->value.complex.i, 0) == 0)
3541         {
3542
3543           mpf_set_ui (result->value.complex.r, 0);
3544           mpf_set_ui (result->value.complex.i, 0);
3545           break;
3546         }
3547
3548       mpf_abs (ac, e->value.complex.r);
3549       mpf_abs (ad, e->value.complex.i);
3550
3551       if (mpf_cmp (ac, ad) >= 0)
3552         {
3553           mpf_div (t, e->value.complex.i, e->value.complex.r);
3554           mpf_mul (t, t, t);
3555           mpf_add_ui (t, t, 1);
3556           mpf_sqrt (t, t);
3557           mpf_add_ui (t, t, 1);
3558           mpf_div_ui (t, t, 2);
3559           mpf_sqrt (t, t);
3560           mpf_sqrt (s, ac);
3561           mpf_mul (w, s, t);
3562         }
3563       else
3564         {
3565           mpf_div (s, e->value.complex.r, e->value.complex.i);
3566           mpf_mul (t, s, s);
3567           mpf_add_ui (t, t, 1);
3568           mpf_sqrt (t, t);
3569           mpf_abs (s, s);
3570           mpf_add (t, t, s);
3571           mpf_div_ui (t, t, 2);
3572           mpf_sqrt (t, t);
3573           mpf_sqrt (s, ad);
3574           mpf_mul (w, s, t);
3575         }
3576
3577       if (mpf_cmp_ui (w, 0) != 0 && mpf_cmp_ui (e->value.complex.r, 0) >= 0)
3578         {
3579           mpf_mul_ui (t, w, 2);
3580           mpf_div (result->value.complex.i, e->value.complex.i, t);
3581           mpf_set (result->value.complex.r, w);
3582         }
3583       else if (mpf_cmp_ui (w, 0) != 0
3584                && mpf_cmp_ui (e->value.complex.r, 0) < 0
3585                && mpf_cmp_ui (e->value.complex.i, 0) >= 0)
3586         {
3587           mpf_mul_ui (t, w, 2);
3588           mpf_div (result->value.complex.r, e->value.complex.i, t);
3589           mpf_set (result->value.complex.i, w);
3590         }
3591       else if (mpf_cmp_ui (w, 0) != 0
3592                && mpf_cmp_ui (e->value.complex.r, 0) < 0
3593                && mpf_cmp_ui (e->value.complex.i, 0) < 0)
3594         {
3595           mpf_mul_ui (t, w, 2);
3596           mpf_div (result->value.complex.r, ad, t);
3597           mpf_neg (w, w);
3598           mpf_set (result->value.complex.i, w);
3599         }
3600       else
3601         gfc_internal_error ("invalid complex argument of SQRT at %L",
3602                             &e->where);
3603
3604       mpf_clear (s);
3605       mpf_clear (t);
3606       mpf_clear (ac);
3607       mpf_clear (ad);
3608       mpf_clear (w);
3609
3610       break;
3611
3612     default:
3613       gfc_internal_error ("invalid argument of SQRT at %L", &e->where);
3614     }
3615
3616   return range_check (result, "SQRT");
3617
3618 negative_arg:
3619   gfc_free_expr (result);
3620   gfc_error ("Argument of SQRT at %L has a negative value", &e->where);
3621   return &gfc_bad_expr;
3622 }
3623
3624
3625 gfc_expr *
3626 gfc_simplify_tan (gfc_expr * x)
3627 {
3628   gfc_expr *result;
3629   mpf_t mpf_sin, mpf_cos, mag_cos;
3630   int i;
3631
3632   if (x->expr_type != EXPR_CONSTANT)
3633     return NULL;
3634
3635   i = gfc_validate_kind (BT_REAL, x->ts.kind);
3636   if (i == -1)
3637     gfc_internal_error ("gfc_simplify_tan(): Bad kind");
3638
3639   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3640
3641   mpf_init (mpf_sin);
3642   mpf_init (mpf_cos);
3643   mpf_init (mag_cos);
3644   sine (&x->value.real, &mpf_sin);
3645   cosine (&x->value.real, &mpf_cos);
3646   mpf_abs (mag_cos, mpf_cos);
3647   if (mpf_cmp_ui (mag_cos, 0) == 0)
3648     {
3649       gfc_error ("Tangent undefined at %L", &x->where);
3650       mpf_clear (mpf_sin);
3651       mpf_clear (mpf_cos);
3652       mpf_clear (mag_cos);
3653       gfc_free_expr (result);
3654       return &gfc_bad_expr;
3655     }
3656   else if (mpf_cmp (mag_cos, gfc_real_kinds[i].tiny) < 0)
3657     {
3658       gfc_error ("Tangent cannot be accurately evaluated at %L", &x->where);
3659       mpf_clear (mpf_sin);
3660       mpf_clear (mpf_cos);
3661       mpf_clear (mag_cos);
3662       gfc_free_expr (result);
3663       return &gfc_bad_expr;
3664     }
3665   else
3666     {
3667       mpf_div (result->value.real, mpf_sin, mpf_cos);
3668       mpf_clear (mpf_sin);
3669       mpf_clear (mpf_cos);
3670       mpf_clear (mag_cos);
3671     }
3672
3673   return range_check (result, "TAN");
3674 }
3675
3676
3677 gfc_expr *
3678 gfc_simplify_tanh (gfc_expr * x)
3679 {
3680   gfc_expr *result;
3681   mpf_t xp, xq;
3682
3683   if (x->expr_type != EXPR_CONSTANT)
3684     return NULL;
3685
3686   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3687
3688   mpf_init (xp);
3689   mpf_init (xq);
3690
3691   hypersine (&x->value.real, &xq);
3692   hypercos (&x->value.real, &xp);
3693
3694   mpf_div (result->value.real, xq, xp);
3695
3696   mpf_clear (xp);
3697   mpf_clear (xq);
3698
3699   return range_check (result, "TANH");
3700
3701 }
3702
3703
3704 gfc_expr *
3705 gfc_simplify_tiny (gfc_expr * e)
3706 {
3707   gfc_expr *result;
3708   int i;
3709
3710   i = gfc_validate_kind (BT_REAL, e->ts.kind);
3711   if (i == -1)
3712     gfc_internal_error ("gfc_simplify_error(): Bad kind");
3713
3714   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
3715   mpf_set (result->value.real, gfc_real_kinds[i].tiny);
3716
3717   return result;
3718 }
3719
3720
3721 gfc_expr *
3722 gfc_simplify_trim (gfc_expr * e)
3723 {
3724   gfc_expr *result;
3725   int count, i, len, lentrim;
3726
3727   if (e->expr_type != EXPR_CONSTANT)
3728     return NULL;
3729
3730   len = e->value.character.length;
3731
3732   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
3733
3734   for (count = 0, i = 1; i <= len; ++i)
3735     {
3736       if (e->value.character.string[len - i] == ' ')
3737         count++;
3738       else
3739         break;
3740     }
3741
3742   lentrim = len - count;
3743
3744   result->value.character.length = lentrim;
3745   result->value.character.string = gfc_getmem (lentrim + 1);
3746
3747   for (i = 0; i < lentrim; i++)
3748     result->value.character.string[i] = e->value.character.string[i];
3749
3750   result->value.character.string[lentrim] = '\0';       /* For debugger */
3751
3752   return result;
3753 }
3754
3755
3756 gfc_expr *
3757 gfc_simplify_ubound (gfc_expr * array, gfc_expr * dim)
3758 {
3759   return gfc_simplify_bound (array, dim, 1);
3760 }
3761
3762
3763 gfc_expr *
3764 gfc_simplify_verify (gfc_expr * s, gfc_expr * set, gfc_expr * b)
3765 {
3766   gfc_expr *result;
3767   int back;
3768   size_t index, len, lenset;
3769   size_t i;
3770
3771   if (s->expr_type != EXPR_CONSTANT || set->expr_type != EXPR_CONSTANT)
3772     return NULL;
3773
3774   if (b != NULL && b->value.logical != 0)
3775     back = 1;
3776   else
3777     back = 0;
3778
3779   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
3780                                 &s->where);
3781
3782   len = s->value.character.length;
3783   lenset = set->value.character.length;
3784
3785   if (len == 0)
3786     {
3787       mpz_set_ui (result->value.integer, 0);
3788       return result;
3789     }
3790
3791   if (back == 0)
3792     {
3793       if (lenset == 0)
3794         {
3795           mpz_set_ui (result->value.integer, len);
3796           return result;
3797         }
3798
3799       index =
3800         strspn (s->value.character.string, set->value.character.string) + 1;
3801       if (index > len)
3802         index = 0;
3803
3804     }
3805   else
3806     {
3807       if (lenset == 0)
3808         {
3809           mpz_set_ui (result->value.integer, 1);
3810           return result;
3811         }
3812       for (index = len; index > 0; index --)
3813         {
3814           for (i = 0; i < lenset; i++)
3815             {
3816               if (s->value.character.string[index - 1]
3817                     == set->value.character.string[i])
3818                 break;
3819             }
3820           if (i == lenset)
3821             break;
3822         }
3823     }
3824
3825   mpz_set_ui (result->value.integer, index);
3826   return result;
3827 }
3828
3829 /****************** Constant simplification *****************/
3830
3831 /* Master function to convert one constant to another.  While this is
3832    used as a simplification function, it requires the destination type
3833    and kind information which is supplied by a special case in
3834    do_simplify().  */
3835
3836 gfc_expr *
3837 gfc_convert_constant (gfc_expr * e, bt type, int kind)
3838 {
3839   gfc_expr *g, *result, *(*f) (gfc_expr *, int);
3840   gfc_constructor *head, *c, *tail = NULL;
3841
3842   switch (e->ts.type)
3843     {
3844     case BT_INTEGER:
3845       switch (type)
3846         {
3847         case BT_INTEGER:
3848           f = gfc_int2int;
3849           break;
3850         case BT_REAL:
3851           f = gfc_int2real;
3852           break;
3853         case BT_COMPLEX:
3854           f = gfc_int2complex;
3855           break;
3856         default:
3857           goto oops;
3858         }
3859       break;
3860
3861     case BT_REAL:
3862       switch (type)
3863         {
3864         case BT_INTEGER:
3865           f = gfc_real2int;
3866           break;
3867         case BT_REAL:
3868           f = gfc_real2real;
3869           break;
3870         case BT_COMPLEX:
3871           f = gfc_real2complex;
3872           break;
3873         default:
3874           goto oops;
3875         }
3876       break;
3877
3878     case BT_COMPLEX:
3879       switch (type)
3880         {
3881         case BT_INTEGER:
3882           f = gfc_complex2int;
3883           break;
3884         case BT_REAL:
3885           f = gfc_complex2real;
3886           break;
3887         case BT_COMPLEX:
3888           f = gfc_complex2complex;
3889           break;
3890
3891         default:
3892           goto oops;
3893         }
3894       break;
3895
3896     case BT_LOGICAL:
3897       if (type != BT_LOGICAL)
3898         goto oops;
3899       f = gfc_log2log;
3900       break;
3901
3902     default:
3903     oops:
3904       gfc_internal_error ("gfc_convert_constant(): Unexpected type");
3905     }
3906
3907   result = NULL;
3908
3909   switch (e->expr_type)
3910     {
3911     case EXPR_CONSTANT:
3912       result = f (e, kind);
3913       if (result == NULL)
3914         return &gfc_bad_expr;
3915       break;
3916
3917     case EXPR_ARRAY:
3918       if (!gfc_is_constant_expr (e))
3919         break;
3920
3921       head = NULL;
3922
3923       for (c = e->value.constructor; c; c = c->next)
3924         {
3925           if (head == NULL)
3926             head = tail = gfc_get_constructor ();
3927           else
3928             {
3929               tail->next = gfc_get_constructor ();
3930               tail = tail->next;
3931             }
3932
3933           tail->where = c->where;
3934
3935           if (c->iterator == NULL)
3936             tail->expr = f (c->expr, kind);
3937           else
3938             {
3939               g = gfc_convert_constant (c->expr, type, kind);
3940               if (g == &gfc_bad_expr)
3941                 return g;
3942               tail->expr = g;
3943             }
3944
3945           if (tail->expr == NULL)
3946             {
3947               gfc_free_constructor (head);
3948               return NULL;
3949             }
3950         }
3951
3952       result = gfc_get_expr ();
3953       result->ts.type = type;
3954       result->ts.kind = kind;
3955       result->expr_type = EXPR_ARRAY;
3956       result->value.constructor = head;
3957       result->shape = gfc_copy_shape (e->shape, e->rank);
3958       result->where = e->where;
3959       result->rank = e->rank;
3960       break;
3961
3962     default:
3963       break;
3964     }
3965
3966   return result;
3967 }
3968
3969
3970 /****************** Helper functions ***********************/
3971
3972 /* Given a collating table, create the inverse table.  */
3973
3974 static void
3975 invert_table (const int *table, int *xtable)
3976 {
3977   int i;
3978
3979   for (i = 0; i < 256; i++)
3980     xtable[i] = 0;
3981
3982   for (i = 0; i < 256; i++)
3983     xtable[table[i]] = i;
3984 }
3985
3986
3987 void
3988 gfc_simplify_init_1 (void)
3989 {
3990
3991   mpf_init_set_str (mpf_zero, "0.0", 10);
3992   mpf_init_set_str (mpf_half, "0.5", 10);
3993   mpf_init_set_str (mpf_one, "1.0", 10);
3994   mpz_init_set_str (mpz_zero, "0", 10);
3995
3996   invert_table (ascii_table, xascii_table);
3997 }
3998
3999
4000 void
4001 gfc_simplify_done_1 (void)
4002 {
4003
4004   mpf_clear (mpf_zero);
4005   mpf_clear (mpf_half);
4006   mpf_clear (mpf_one);
4007   mpz_clear (mpz_zero);
4008 }