OSDN Git Service

2014-03-09 Janus Weil <janus@gcc.gnu.org>
[pf3gnuchains/gcc-fork.git] / gcc / fortran / simplify.c
1 /* Simplify intrinsic functions at compile-time.
2    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009,
3    2010, 2011 Free Software Foundation, 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 3, 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 COPYING3.  If not see
20 <http://www.gnu.org/licenses/>.  */
21
22 #include "config.h"
23 #include "system.h"
24 #include "flags.h"
25 #include "gfortran.h"
26 #include "arith.h"
27 #include "intrinsic.h"
28 #include "target-memory.h"
29 #include "constructor.h"
30 #include "version.h"  /* For version_string.  */
31
32
33 gfc_expr gfc_bad_expr;
34
35 static gfc_expr *simplify_size (gfc_expr *, gfc_expr *, int);
36
37
38 /* Note that 'simplification' is not just transforming expressions.
39    For functions that are not simplified at compile time, range
40    checking is done if possible.
41
42    The return convention is that each simplification function returns:
43
44      A new expression node corresponding to the simplified arguments.
45      The original arguments are destroyed by the caller, and must not
46      be a part of the new expression.
47
48      NULL pointer indicating that no simplification was possible and
49      the original expression should remain intact.
50
51      An expression pointer to gfc_bad_expr (a static placeholder)
52      indicating that some error has prevented simplification.  The
53      error is generated within the function and should be propagated
54      upwards
55
56    By the time a simplification function gets control, it has been
57    decided that the function call is really supposed to be the
58    intrinsic.  No type checking is strictly necessary, since only
59    valid types will be passed on.  On the other hand, a simplification
60    subroutine may have to look at the type of an argument as part of
61    its processing.
62
63    Array arguments are only passed to these subroutines that implement
64    the simplification of transformational intrinsics.
65
66    The functions in this file don't have much comment with them, but
67    everything is reasonably straight-forward.  The Standard, chapter 13
68    is the best comment you'll find for this file anyway.  */
69
70 /* Range checks an expression node.  If all goes well, returns the
71    node, otherwise returns &gfc_bad_expr and frees the node.  */
72
73 static gfc_expr *
74 range_check (gfc_expr *result, const char *name)
75 {
76   if (result == NULL)
77     return &gfc_bad_expr;
78
79   if (result->expr_type != EXPR_CONSTANT)
80     return result;
81
82   switch (gfc_range_check (result))
83     {
84       case ARITH_OK:
85         return result;
86  
87       case ARITH_OVERFLOW:
88         gfc_error ("Result of %s overflows its kind at %L", name,
89                    &result->where);
90         break;
91
92       case ARITH_UNDERFLOW:
93         gfc_error ("Result of %s underflows its kind at %L", name,
94                    &result->where);
95         break;
96
97       case ARITH_NAN:
98         gfc_error ("Result of %s is NaN at %L", name, &result->where);
99         break;
100
101       default:
102         gfc_error ("Result of %s gives range error for its kind at %L", name,
103                    &result->where);
104         break;
105     }
106
107   gfc_free_expr (result);
108   return &gfc_bad_expr;
109 }
110
111
112 /* A helper function that gets an optional and possibly missing
113    kind parameter.  Returns the kind, -1 if something went wrong.  */
114
115 static int
116 get_kind (bt type, gfc_expr *k, const char *name, int default_kind)
117 {
118   int kind;
119
120   if (k == NULL)
121     return default_kind;
122
123   if (k->expr_type != EXPR_CONSTANT)
124     {
125       gfc_error ("KIND parameter of %s at %L must be an initialization "
126                  "expression", name, &k->where);
127       return -1;
128     }
129
130   if (gfc_extract_int (k, &kind) != NULL
131       || gfc_validate_kind (type, kind, true) < 0)
132     {
133       gfc_error ("Invalid KIND parameter of %s at %L", name, &k->where);
134       return -1;
135     }
136
137   return kind;
138 }
139
140
141 /* Converts an mpz_t signed variable into an unsigned one, assuming
142    two's complement representations and a binary width of bitsize.
143    The conversion is a no-op unless x is negative; otherwise, it can
144    be accomplished by masking out the high bits.  */
145
146 static void
147 convert_mpz_to_unsigned (mpz_t x, int bitsize)
148 {
149   mpz_t mask;
150
151   if (mpz_sgn (x) < 0)
152     {
153       /* Confirm that no bits above the signed range are unset.  */
154       gcc_assert (mpz_scan0 (x, bitsize-1) == ULONG_MAX);
155
156       mpz_init_set_ui (mask, 1);
157       mpz_mul_2exp (mask, mask, bitsize);
158       mpz_sub_ui (mask, mask, 1);
159
160       mpz_and (x, x, mask);
161
162       mpz_clear (mask);
163     }
164   else
165     {
166       /* Confirm that no bits above the signed range are set.  */
167       gcc_assert (mpz_scan1 (x, bitsize-1) == ULONG_MAX);
168     }
169 }
170
171
172 /* Converts an mpz_t unsigned variable into a signed one, assuming
173    two's complement representations and a binary width of bitsize.
174    If the bitsize-1 bit is set, this is taken as a sign bit and
175    the number is converted to the corresponding negative number.  */
176
177 static void
178 convert_mpz_to_signed (mpz_t x, int bitsize)
179 {
180   mpz_t mask;
181
182   /* Confirm that no bits above the unsigned range are set.  */
183   gcc_assert (mpz_scan1 (x, bitsize) == ULONG_MAX);
184
185   if (mpz_tstbit (x, bitsize - 1) == 1)
186     {
187       mpz_init_set_ui (mask, 1);
188       mpz_mul_2exp (mask, mask, bitsize);
189       mpz_sub_ui (mask, mask, 1);
190
191       /* We negate the number by hand, zeroing the high bits, that is
192          make it the corresponding positive number, and then have it
193          negated by GMP, giving the correct representation of the
194          negative number.  */
195       mpz_com (x, x);
196       mpz_add_ui (x, x, 1);
197       mpz_and (x, x, mask);
198
199       mpz_neg (x, x);
200
201       mpz_clear (mask);
202     }
203 }
204
205
206 /* In-place convert BOZ to REAL of the specified kind.  */
207
208 static gfc_expr *
209 convert_boz (gfc_expr *x, int kind)
210 {
211   if (x && x->ts.type == BT_INTEGER && x->is_boz)
212     {
213       gfc_typespec ts;
214       gfc_clear_ts (&ts);
215       ts.type = BT_REAL;
216       ts.kind = kind;
217
218       if (!gfc_convert_boz (x, &ts))
219         return &gfc_bad_expr;
220     }
221
222   return x;
223 }
224
225
226 /* Test that the expression is an constant array.  */
227
228 static bool
229 is_constant_array_expr (gfc_expr *e)
230 {
231   gfc_constructor *c;
232
233   if (e == NULL)
234     return true;
235
236   if (e->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (e))
237     return false;
238
239   for (c = gfc_constructor_first (e->value.constructor);
240        c; c = gfc_constructor_next (c))
241     if (c->expr->expr_type != EXPR_CONSTANT
242           && c->expr->expr_type != EXPR_STRUCTURE)
243       return false;
244
245   return true;
246 }
247
248
249 /* Initialize a transformational result expression with a given value.  */
250
251 static void
252 init_result_expr (gfc_expr *e, int init, gfc_expr *array)
253 {
254   if (e && e->expr_type == EXPR_ARRAY)
255     {
256       gfc_constructor *ctor = gfc_constructor_first (e->value.constructor);
257       while (ctor)
258         {
259           init_result_expr (ctor->expr, init, array);
260           ctor = gfc_constructor_next (ctor);
261         }
262     }
263   else if (e && e->expr_type == EXPR_CONSTANT)
264     {
265       int i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
266       int length;
267       gfc_char_t *string;
268
269       switch (e->ts.type)
270         {
271           case BT_LOGICAL:
272             e->value.logical = (init ? 1 : 0);
273             break;
274
275           case BT_INTEGER:
276             if (init == INT_MIN)
277               mpz_set (e->value.integer, gfc_integer_kinds[i].min_int);
278             else if (init == INT_MAX)
279               mpz_set (e->value.integer, gfc_integer_kinds[i].huge);
280             else
281               mpz_set_si (e->value.integer, init);
282             break;
283
284           case BT_REAL:
285             if (init == INT_MIN)
286               {
287                 mpfr_set (e->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
288                 mpfr_neg (e->value.real, e->value.real, GFC_RND_MODE);
289               }
290             else if (init == INT_MAX)
291               mpfr_set (e->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
292             else
293               mpfr_set_si (e->value.real, init, GFC_RND_MODE);
294             break;
295
296           case BT_COMPLEX:
297             mpc_set_si (e->value.complex, init, GFC_MPC_RND_MODE);
298             break;
299
300           case BT_CHARACTER:
301             if (init == INT_MIN)
302               {
303                 gfc_expr *len = gfc_simplify_len (array, NULL);
304                 gfc_extract_int (len, &length);
305                 string = gfc_get_wide_string (length + 1);
306                 gfc_wide_memset (string, 0, length);
307               }
308             else if (init == INT_MAX)
309               {
310                 gfc_expr *len = gfc_simplify_len (array, NULL);
311                 gfc_extract_int (len, &length);
312                 string = gfc_get_wide_string (length + 1);
313                 gfc_wide_memset (string, 255, length);
314               }
315             else
316               {
317                 length = 0;
318                 string = gfc_get_wide_string (1);
319               }
320
321             string[length] = '\0';
322             e->value.character.length = length;
323             e->value.character.string = string;
324             break;
325
326           default:
327             gcc_unreachable();
328         }
329     }
330   else
331     gcc_unreachable();
332 }
333
334
335 /* Helper function for gfc_simplify_dot_product() and gfc_simplify_matmul;
336    if conj_a is true, the matrix_a is complex conjugated.  */
337
338 static gfc_expr *
339 compute_dot_product (gfc_expr *matrix_a, int stride_a, int offset_a,
340                      gfc_expr *matrix_b, int stride_b, int offset_b,
341                      bool conj_a)
342 {
343   gfc_expr *result, *a, *b, *c;
344
345   result = gfc_get_constant_expr (matrix_a->ts.type, matrix_a->ts.kind,
346                                   &matrix_a->where);
347   init_result_expr (result, 0, NULL);
348
349   a = gfc_constructor_lookup_expr (matrix_a->value.constructor, offset_a);
350   b = gfc_constructor_lookup_expr (matrix_b->value.constructor, offset_b);
351   while (a && b)
352     {
353       /* Copying of expressions is required as operands are free'd
354          by the gfc_arith routines.  */
355       switch (result->ts.type)
356         {
357           case BT_LOGICAL:
358             result = gfc_or (result,
359                              gfc_and (gfc_copy_expr (a),
360                                       gfc_copy_expr (b)));
361             break;
362
363           case BT_INTEGER:
364           case BT_REAL:
365           case BT_COMPLEX:
366             if (conj_a && a->ts.type == BT_COMPLEX)
367               c = gfc_simplify_conjg (a);
368             else
369               c = gfc_copy_expr (a);
370             result = gfc_add (result, gfc_multiply (c, gfc_copy_expr (b)));
371             break;
372
373           default:
374             gcc_unreachable();
375         }
376
377       offset_a += stride_a;
378       a = gfc_constructor_lookup_expr (matrix_a->value.constructor, offset_a);
379
380       offset_b += stride_b;
381       b = gfc_constructor_lookup_expr (matrix_b->value.constructor, offset_b);
382     }
383
384   return result;
385 }
386
387
388 /* Build a result expression for transformational intrinsics, 
389    depending on DIM. */
390
391 static gfc_expr *
392 transformational_result (gfc_expr *array, gfc_expr *dim, bt type,
393                          int kind, locus* where)
394 {
395   gfc_expr *result;
396   int i, nelem;
397
398   if (!dim || array->rank == 1)
399     return gfc_get_constant_expr (type, kind, where);
400
401   result = gfc_get_array_expr (type, kind, where);
402   result->shape = gfc_copy_shape_excluding (array->shape, array->rank, dim);
403   result->rank = array->rank - 1;
404
405   /* gfc_array_size() would count the number of elements in the constructor,
406      we have not built those yet.  */
407   nelem = 1;
408   for  (i = 0; i < result->rank; ++i)
409     nelem *= mpz_get_ui (result->shape[i]);
410
411   for (i = 0; i < nelem; ++i)
412     {
413       gfc_constructor_append_expr (&result->value.constructor,
414                                    gfc_get_constant_expr (type, kind, where),
415                                    NULL);
416     }
417
418   return result;
419 }
420
421
422 typedef gfc_expr* (*transformational_op)(gfc_expr*, gfc_expr*);
423
424 /* Wrapper function, implements 'op1 += 1'. Only called if MASK
425    of COUNT intrinsic is .TRUE..
426
427    Interface and implimentation mimics arith functions as
428    gfc_add, gfc_multiply, etc.  */
429
430 static gfc_expr* gfc_count (gfc_expr *op1, gfc_expr *op2)
431 {
432   gfc_expr *result;
433
434   gcc_assert (op1->ts.type == BT_INTEGER);
435   gcc_assert (op2->ts.type == BT_LOGICAL);
436   gcc_assert (op2->value.logical);
437
438   result = gfc_copy_expr (op1);
439   mpz_add_ui (result->value.integer, result->value.integer, 1);
440
441   gfc_free_expr (op1);
442   gfc_free_expr (op2);
443   return result;
444 }
445
446
447 /* Transforms an ARRAY with operation OP, according to MASK, to a
448    scalar RESULT. E.g. called if
449
450      REAL, PARAMETER :: array(n, m) = ...
451      REAL, PARAMETER :: s = SUM(array)
452
453   where OP == gfc_add().  */
454
455 static gfc_expr *
456 simplify_transformation_to_scalar (gfc_expr *result, gfc_expr *array, gfc_expr *mask,
457                                    transformational_op op)
458 {
459   gfc_expr *a, *m;
460   gfc_constructor *array_ctor, *mask_ctor;
461
462   /* Shortcut for constant .FALSE. MASK.  */
463   if (mask
464       && mask->expr_type == EXPR_CONSTANT
465       && !mask->value.logical)
466     return result;
467
468   array_ctor = gfc_constructor_first (array->value.constructor);
469   mask_ctor = NULL;
470   if (mask && mask->expr_type == EXPR_ARRAY)
471     mask_ctor = gfc_constructor_first (mask->value.constructor);
472
473   while (array_ctor)
474     {
475       a = array_ctor->expr;
476       array_ctor = gfc_constructor_next (array_ctor);
477
478       /* A constant MASK equals .TRUE. here and can be ignored.  */
479       if (mask_ctor)
480         {
481           m = mask_ctor->expr;
482           mask_ctor = gfc_constructor_next (mask_ctor);
483           if (!m->value.logical)
484             continue;
485         }
486
487       result = op (result, gfc_copy_expr (a));
488     }
489
490   return result;
491 }
492
493 /* Transforms an ARRAY with operation OP, according to MASK, to an
494    array RESULT. E.g. called if
495
496      REAL, PARAMETER :: array(n, m) = ...
497      REAL, PARAMETER :: s(n) = PROD(array, DIM=1)
498
499   where OP == gfc_multiply(). The result might be post processed using post_op. */ 
500
501 static gfc_expr *
502 simplify_transformation_to_array (gfc_expr *result, gfc_expr *array, gfc_expr *dim,
503                                   gfc_expr *mask, transformational_op op,
504                                   transformational_op post_op)
505 {
506   mpz_t size;
507   int done, i, n, arraysize, resultsize, dim_index, dim_extent, dim_stride;
508   gfc_expr **arrayvec, **resultvec, **base, **src, **dest;
509   gfc_constructor *array_ctor, *mask_ctor, *result_ctor;
510
511   int count[GFC_MAX_DIMENSIONS], extent[GFC_MAX_DIMENSIONS],
512       sstride[GFC_MAX_DIMENSIONS], dstride[GFC_MAX_DIMENSIONS],
513       tmpstride[GFC_MAX_DIMENSIONS];
514
515   /* Shortcut for constant .FALSE. MASK.  */
516   if (mask
517       && mask->expr_type == EXPR_CONSTANT
518       && !mask->value.logical)
519     return result;
520
521   /* Build an indexed table for array element expressions to minimize
522      linked-list traversal. Masked elements are set to NULL.  */
523   gfc_array_size (array, &size);
524   arraysize = mpz_get_ui (size);
525   mpz_clear (size);
526
527   arrayvec = XCNEWVEC (gfc_expr*, arraysize);
528
529   array_ctor = gfc_constructor_first (array->value.constructor);
530   mask_ctor = NULL;
531   if (mask && mask->expr_type == EXPR_ARRAY)
532     mask_ctor = gfc_constructor_first (mask->value.constructor);
533
534   for (i = 0; i < arraysize; ++i)
535     {
536       arrayvec[i] = array_ctor->expr;
537       array_ctor = gfc_constructor_next (array_ctor);
538
539       if (mask_ctor)
540         {
541           if (!mask_ctor->expr->value.logical)
542             arrayvec[i] = NULL;
543
544           mask_ctor = gfc_constructor_next (mask_ctor);
545         }
546     }
547
548   /* Same for the result expression.  */
549   gfc_array_size (result, &size);
550   resultsize = mpz_get_ui (size);
551   mpz_clear (size);
552
553   resultvec = XCNEWVEC (gfc_expr*, resultsize);
554   result_ctor = gfc_constructor_first (result->value.constructor);
555   for (i = 0; i < resultsize; ++i)
556     {
557       resultvec[i] = result_ctor->expr;
558       result_ctor = gfc_constructor_next (result_ctor);
559     }
560
561   gfc_extract_int (dim, &dim_index);
562   dim_index -= 1;               /* zero-base index */
563   dim_extent = 0;
564   dim_stride = 0;
565
566   for (i = 0, n = 0; i < array->rank; ++i)
567     {
568       count[i] = 0;
569       tmpstride[i] = (i == 0) ? 1 : tmpstride[i-1] * mpz_get_si (array->shape[i-1]);
570       if (i == dim_index)
571         {
572           dim_extent = mpz_get_si (array->shape[i]);
573           dim_stride = tmpstride[i];
574           continue;
575         }
576
577       extent[n] = mpz_get_si (array->shape[i]);
578       sstride[n] = tmpstride[i];
579       dstride[n] = (n == 0) ? 1 : dstride[n-1] * extent[n-1];
580       n += 1;
581     }
582
583   done = false;
584   base = arrayvec;
585   dest = resultvec;
586   while (!done)
587     {
588       for (src = base, n = 0; n < dim_extent; src += dim_stride, ++n)
589         if (*src)
590           *dest = op (*dest, gfc_copy_expr (*src));
591
592       count[0]++;
593       base += sstride[0];
594       dest += dstride[0];
595
596       n = 0;
597       while (!done && count[n] == extent[n])
598         {
599           count[n] = 0;
600           base -= sstride[n] * extent[n];
601           dest -= dstride[n] * extent[n];
602
603           n++;
604           if (n < result->rank)
605             {
606               count [n]++;
607               base += sstride[n];
608               dest += dstride[n];
609             }
610           else
611             done = true;
612        }
613     }
614
615   /* Place updated expression in result constructor.  */
616   result_ctor = gfc_constructor_first (result->value.constructor);
617   for (i = 0; i < resultsize; ++i)
618     {
619       if (post_op)
620         result_ctor->expr = post_op (result_ctor->expr, resultvec[i]);
621       else
622         result_ctor->expr = resultvec[i];
623       result_ctor = gfc_constructor_next (result_ctor);
624     }
625
626   free (arrayvec);
627   free (resultvec);
628   return result;
629 }
630
631
632 static gfc_expr *
633 simplify_transformation (gfc_expr *array, gfc_expr *dim, gfc_expr *mask,
634                          int init_val, transformational_op op)
635 {
636   gfc_expr *result;
637
638   if (!is_constant_array_expr (array)
639       || !gfc_is_constant_expr (dim))
640     return NULL;
641
642   if (mask
643       && !is_constant_array_expr (mask)
644       && mask->expr_type != EXPR_CONSTANT)
645     return NULL;
646
647   result = transformational_result (array, dim, array->ts.type,
648                                     array->ts.kind, &array->where);
649   init_result_expr (result, init_val, NULL);
650
651   return !dim || array->rank == 1 ?
652     simplify_transformation_to_scalar (result, array, mask, op) :
653     simplify_transformation_to_array (result, array, dim, mask, op, NULL);
654 }
655
656
657 /********************** Simplification functions *****************************/
658
659 gfc_expr *
660 gfc_simplify_abs (gfc_expr *e)
661 {
662   gfc_expr *result;
663
664   if (e->expr_type != EXPR_CONSTANT)
665     return NULL;
666
667   switch (e->ts.type)
668     {
669       case BT_INTEGER:
670         result = gfc_get_constant_expr (BT_INTEGER, e->ts.kind, &e->where);
671         mpz_abs (result->value.integer, e->value.integer);
672         return range_check (result, "IABS");
673
674       case BT_REAL:
675         result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
676         mpfr_abs (result->value.real, e->value.real, GFC_RND_MODE);
677         return range_check (result, "ABS");
678
679       case BT_COMPLEX:
680         gfc_set_model_kind (e->ts.kind);
681         result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
682         mpc_abs (result->value.real, e->value.complex, GFC_RND_MODE);
683         return range_check (result, "CABS");
684
685       default:
686         gfc_internal_error ("gfc_simplify_abs(): Bad type");
687     }
688 }
689
690
691 static gfc_expr *
692 simplify_achar_char (gfc_expr *e, gfc_expr *k, const char *name, bool ascii)
693 {
694   gfc_expr *result;
695   int kind;
696   bool too_large = false;
697
698   if (e->expr_type != EXPR_CONSTANT)
699     return NULL;
700
701   kind = get_kind (BT_CHARACTER, k, name, gfc_default_character_kind);
702   if (kind == -1)
703     return &gfc_bad_expr;
704
705   if (mpz_cmp_si (e->value.integer, 0) < 0)
706     {
707       gfc_error ("Argument of %s function at %L is negative", name,
708                  &e->where);
709       return &gfc_bad_expr;
710     }
711
712   if (ascii && gfc_option.warn_surprising
713       && mpz_cmp_si (e->value.integer, 127) > 0)
714     gfc_warning ("Argument of %s function at %L outside of range [0,127]",
715                  name, &e->where);
716
717   if (kind == 1 && mpz_cmp_si (e->value.integer, 255) > 0)
718     too_large = true;
719   else if (kind == 4)
720     {
721       mpz_t t;
722       mpz_init_set_ui (t, 2);
723       mpz_pow_ui (t, t, 32);
724       mpz_sub_ui (t, t, 1);
725       if (mpz_cmp (e->value.integer, t) > 0)
726         too_large = true;
727       mpz_clear (t);
728     }
729
730   if (too_large)
731     {
732       gfc_error ("Argument of %s function at %L is too large for the "
733                  "collating sequence of kind %d", name, &e->where, kind);
734       return &gfc_bad_expr;
735     }
736
737   result = gfc_get_character_expr (kind, &e->where, NULL, 1);
738   result->value.character.string[0] = mpz_get_ui (e->value.integer);
739
740   return result;
741 }
742
743
744
745 /* We use the processor's collating sequence, because all
746    systems that gfortran currently works on are ASCII.  */
747
748 gfc_expr *
749 gfc_simplify_achar (gfc_expr *e, gfc_expr *k)
750 {
751   return simplify_achar_char (e, k, "ACHAR", true);
752 }
753
754
755 gfc_expr *
756 gfc_simplify_acos (gfc_expr *x)
757 {
758   gfc_expr *result;
759
760   if (x->expr_type != EXPR_CONSTANT)
761     return NULL;
762
763   switch (x->ts.type)
764     {
765       case BT_REAL:
766         if (mpfr_cmp_si (x->value.real, 1) > 0
767             || mpfr_cmp_si (x->value.real, -1) < 0)
768           {
769             gfc_error ("Argument of ACOS at %L must be between -1 and 1",
770                        &x->where);
771             return &gfc_bad_expr;
772           }
773         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
774         mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
775         break;
776
777       case BT_COMPLEX:
778         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
779         mpc_acos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
780         break;
781
782       default:
783         gfc_internal_error ("in gfc_simplify_acos(): Bad type");
784     }
785
786   return range_check (result, "ACOS");
787 }
788
789 gfc_expr *
790 gfc_simplify_acosh (gfc_expr *x)
791 {
792   gfc_expr *result;
793
794   if (x->expr_type != EXPR_CONSTANT)
795     return NULL;
796
797   switch (x->ts.type)
798     {
799       case BT_REAL:
800         if (mpfr_cmp_si (x->value.real, 1) < 0)
801           {
802             gfc_error ("Argument of ACOSH at %L must not be less than 1",
803                        &x->where);
804             return &gfc_bad_expr;
805           }
806
807         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
808         mpfr_acosh (result->value.real, x->value.real, GFC_RND_MODE);
809         break;
810
811       case BT_COMPLEX:
812         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
813         mpc_acosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
814         break;
815
816       default:
817         gfc_internal_error ("in gfc_simplify_acosh(): Bad type");
818     }
819
820   return range_check (result, "ACOSH");
821 }
822
823 gfc_expr *
824 gfc_simplify_adjustl (gfc_expr *e)
825 {
826   gfc_expr *result;
827   int count, i, len;
828   gfc_char_t ch;
829
830   if (e->expr_type != EXPR_CONSTANT)
831     return NULL;
832
833   len = e->value.character.length;
834
835   for (count = 0, i = 0; i < len; ++i)
836     {
837       ch = e->value.character.string[i];
838       if (ch != ' ')
839         break;
840       ++count;
841     }
842
843   result = gfc_get_character_expr (e->ts.kind, &e->where, NULL, len);
844   for (i = 0; i < len - count; ++i)
845     result->value.character.string[i] = e->value.character.string[count + i];
846
847   return result;
848 }
849
850
851 gfc_expr *
852 gfc_simplify_adjustr (gfc_expr *e)
853 {
854   gfc_expr *result;
855   int count, i, len;
856   gfc_char_t ch;
857
858   if (e->expr_type != EXPR_CONSTANT)
859     return NULL;
860
861   len = e->value.character.length;
862
863   for (count = 0, i = len - 1; i >= 0; --i)
864     {
865       ch = e->value.character.string[i];
866       if (ch != ' ')
867         break;
868       ++count;
869     }
870
871   result = gfc_get_character_expr (e->ts.kind, &e->where, NULL, len);
872   for (i = 0; i < count; ++i)
873     result->value.character.string[i] = ' ';
874
875   for (i = count; i < len; ++i)
876     result->value.character.string[i] = e->value.character.string[i - count];
877
878   return result;
879 }
880
881
882 gfc_expr *
883 gfc_simplify_aimag (gfc_expr *e)
884 {
885   gfc_expr *result;
886
887   if (e->expr_type != EXPR_CONSTANT)
888     return NULL;
889
890   result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
891   mpfr_set (result->value.real, mpc_imagref (e->value.complex), GFC_RND_MODE);
892
893   return range_check (result, "AIMAG");
894 }
895
896
897 gfc_expr *
898 gfc_simplify_aint (gfc_expr *e, gfc_expr *k)
899 {
900   gfc_expr *rtrunc, *result;
901   int kind;
902
903   kind = get_kind (BT_REAL, k, "AINT", e->ts.kind);
904   if (kind == -1)
905     return &gfc_bad_expr;
906
907   if (e->expr_type != EXPR_CONSTANT)
908     return NULL;
909
910   rtrunc = gfc_copy_expr (e);
911   mpfr_trunc (rtrunc->value.real, e->value.real);
912
913   result = gfc_real2real (rtrunc, kind);
914
915   gfc_free_expr (rtrunc);
916
917   return range_check (result, "AINT");
918 }
919
920
921 gfc_expr *
922 gfc_simplify_all (gfc_expr *mask, gfc_expr *dim)
923 {
924   return simplify_transformation (mask, dim, NULL, true, gfc_and);
925 }
926
927
928 gfc_expr *
929 gfc_simplify_dint (gfc_expr *e)
930 {
931   gfc_expr *rtrunc, *result;
932
933   if (e->expr_type != EXPR_CONSTANT)
934     return NULL;
935
936   rtrunc = gfc_copy_expr (e);
937   mpfr_trunc (rtrunc->value.real, e->value.real);
938
939   result = gfc_real2real (rtrunc, gfc_default_double_kind);
940
941   gfc_free_expr (rtrunc);
942
943   return range_check (result, "DINT");
944 }
945
946
947 gfc_expr *
948 gfc_simplify_dreal (gfc_expr *e)
949 {
950   gfc_expr *result = NULL;
951
952   if (e->expr_type != EXPR_CONSTANT)
953     return NULL;
954
955   result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
956   mpc_real (result->value.real, e->value.complex, GFC_RND_MODE);
957
958   return range_check (result, "DREAL");
959 }
960
961
962 gfc_expr *
963 gfc_simplify_anint (gfc_expr *e, gfc_expr *k)
964 {
965   gfc_expr *result;
966   int kind;
967
968   kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
969   if (kind == -1)
970     return &gfc_bad_expr;
971
972   if (e->expr_type != EXPR_CONSTANT)
973     return NULL;
974
975   result = gfc_get_constant_expr (e->ts.type, kind, &e->where);
976   mpfr_round (result->value.real, e->value.real);
977
978   return range_check (result, "ANINT");
979 }
980
981
982 gfc_expr *
983 gfc_simplify_and (gfc_expr *x, gfc_expr *y)
984 {
985   gfc_expr *result;
986   int kind;
987
988   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
989     return NULL;
990
991   kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
992
993   switch (x->ts.type)
994     {
995       case BT_INTEGER:
996         result = gfc_get_constant_expr (BT_INTEGER, kind, &x->where);
997         mpz_and (result->value.integer, x->value.integer, y->value.integer);
998         return range_check (result, "AND");
999
1000       case BT_LOGICAL:
1001         return gfc_get_logical_expr (kind, &x->where,
1002                                      x->value.logical && y->value.logical);
1003
1004       default:
1005         gcc_unreachable ();
1006     }
1007 }
1008
1009
1010 gfc_expr *
1011 gfc_simplify_any (gfc_expr *mask, gfc_expr *dim)
1012 {
1013   return simplify_transformation (mask, dim, NULL, false, gfc_or);
1014 }
1015
1016
1017 gfc_expr *
1018 gfc_simplify_dnint (gfc_expr *e)
1019 {
1020   gfc_expr *result;
1021
1022   if (e->expr_type != EXPR_CONSTANT)
1023     return NULL;
1024
1025   result = gfc_get_constant_expr (BT_REAL, gfc_default_double_kind, &e->where);
1026   mpfr_round (result->value.real, e->value.real);
1027
1028   return range_check (result, "DNINT");
1029 }
1030
1031
1032 gfc_expr *
1033 gfc_simplify_asin (gfc_expr *x)
1034 {
1035   gfc_expr *result;
1036
1037   if (x->expr_type != EXPR_CONSTANT)
1038     return NULL;
1039
1040   switch (x->ts.type)
1041     {
1042       case BT_REAL:
1043         if (mpfr_cmp_si (x->value.real, 1) > 0
1044             || mpfr_cmp_si (x->value.real, -1) < 0)
1045           {
1046             gfc_error ("Argument of ASIN at %L must be between -1 and 1",
1047                        &x->where);
1048             return &gfc_bad_expr;
1049           }
1050         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1051         mpfr_asin (result->value.real, x->value.real, GFC_RND_MODE);
1052         break;
1053
1054       case BT_COMPLEX:
1055         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1056         mpc_asin (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1057         break;
1058
1059       default:
1060         gfc_internal_error ("in gfc_simplify_asin(): Bad type");
1061     }
1062
1063   return range_check (result, "ASIN");
1064 }
1065
1066
1067 gfc_expr *
1068 gfc_simplify_asinh (gfc_expr *x)
1069 {
1070   gfc_expr *result;
1071
1072   if (x->expr_type != EXPR_CONSTANT)
1073     return NULL;
1074
1075   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1076
1077   switch (x->ts.type)
1078     {
1079       case BT_REAL:
1080         mpfr_asinh (result->value.real, x->value.real, GFC_RND_MODE);
1081         break;
1082
1083       case BT_COMPLEX:
1084         mpc_asinh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1085         break;
1086
1087       default:
1088         gfc_internal_error ("in gfc_simplify_asinh(): Bad type");
1089     }
1090
1091   return range_check (result, "ASINH");
1092 }
1093
1094
1095 gfc_expr *
1096 gfc_simplify_atan (gfc_expr *x)
1097 {
1098   gfc_expr *result;
1099
1100   if (x->expr_type != EXPR_CONSTANT)
1101     return NULL;
1102
1103   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1104
1105   switch (x->ts.type)
1106     {
1107       case BT_REAL:
1108         mpfr_atan (result->value.real, x->value.real, GFC_RND_MODE);
1109         break;
1110
1111       case BT_COMPLEX:
1112         mpc_atan (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1113         break;
1114
1115       default:
1116         gfc_internal_error ("in gfc_simplify_atan(): Bad type");
1117     }
1118
1119   return range_check (result, "ATAN");
1120 }
1121
1122
1123 gfc_expr *
1124 gfc_simplify_atanh (gfc_expr *x)
1125 {
1126   gfc_expr *result;
1127
1128   if (x->expr_type != EXPR_CONSTANT)
1129     return NULL;
1130
1131   switch (x->ts.type)
1132     {
1133       case BT_REAL:
1134         if (mpfr_cmp_si (x->value.real, 1) >= 0
1135             || mpfr_cmp_si (x->value.real, -1) <= 0)
1136           {
1137             gfc_error ("Argument of ATANH at %L must be inside the range -1 "
1138                        "to 1", &x->where);
1139             return &gfc_bad_expr;
1140           }
1141         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1142         mpfr_atanh (result->value.real, x->value.real, GFC_RND_MODE);
1143         break;
1144
1145       case BT_COMPLEX:
1146         result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1147         mpc_atanh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1148         break;
1149
1150       default:
1151         gfc_internal_error ("in gfc_simplify_atanh(): Bad type");
1152     }
1153
1154   return range_check (result, "ATANH");
1155 }
1156
1157
1158 gfc_expr *
1159 gfc_simplify_atan2 (gfc_expr *y, gfc_expr *x)
1160 {
1161   gfc_expr *result;
1162
1163   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1164     return NULL;
1165
1166   if (mpfr_sgn (y->value.real) == 0 && mpfr_sgn (x->value.real) == 0)
1167     {
1168       gfc_error ("If first argument of ATAN2 %L is zero, then the "
1169                  "second argument must not be zero", &x->where);
1170       return &gfc_bad_expr;
1171     }
1172
1173   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1174   mpfr_atan2 (result->value.real, y->value.real, x->value.real, GFC_RND_MODE);
1175
1176   return range_check (result, "ATAN2");
1177 }
1178
1179
1180 gfc_expr *
1181 gfc_simplify_bessel_j0 (gfc_expr *x)
1182 {
1183   gfc_expr *result;
1184
1185   if (x->expr_type != EXPR_CONSTANT)
1186     return NULL;
1187
1188   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1189   mpfr_j0 (result->value.real, x->value.real, GFC_RND_MODE);
1190
1191   return range_check (result, "BESSEL_J0");
1192 }
1193
1194
1195 gfc_expr *
1196 gfc_simplify_bessel_j1 (gfc_expr *x)
1197 {
1198   gfc_expr *result;
1199
1200   if (x->expr_type != EXPR_CONSTANT)
1201     return NULL;
1202
1203   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1204   mpfr_j1 (result->value.real, x->value.real, GFC_RND_MODE);
1205
1206   return range_check (result, "BESSEL_J1");
1207 }
1208
1209
1210 gfc_expr *
1211 gfc_simplify_bessel_jn (gfc_expr *order, gfc_expr *x)
1212 {
1213   gfc_expr *result;
1214   long n;
1215
1216   if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
1217     return NULL;
1218
1219   n = mpz_get_si (order->value.integer);
1220   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1221   mpfr_jn (result->value.real, n, x->value.real, GFC_RND_MODE);
1222
1223   return range_check (result, "BESSEL_JN");
1224 }
1225
1226
1227 /* Simplify transformational form of JN and YN.  */
1228
1229 static gfc_expr *
1230 gfc_simplify_bessel_n2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x,
1231                         bool jn)
1232 {
1233   gfc_expr *result;
1234   gfc_expr *e;
1235   long n1, n2;
1236   int i;
1237   mpfr_t x2rev, last1, last2;
1238
1239   if (x->expr_type != EXPR_CONSTANT || order1->expr_type != EXPR_CONSTANT
1240       || order2->expr_type != EXPR_CONSTANT)
1241     return NULL;
1242
1243   n1 = mpz_get_si (order1->value.integer);
1244   n2 = mpz_get_si (order2->value.integer);
1245   result = gfc_get_array_expr (x->ts.type, x->ts.kind, &x->where);
1246   result->rank = 1;
1247   result->shape = gfc_get_shape (1);
1248   mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));
1249
1250   if (n2 < n1)
1251     return result;
1252
1253   /* Special case: x == 0; it is J0(0.0) == 1, JN(N > 0, 0.0) == 0; and
1254      YN(N, 0.0) = -Inf.  */
1255
1256   if (mpfr_cmp_ui (x->value.real, 0.0) == 0)
1257     {
1258       if (!jn && gfc_option.flag_range_check)
1259         {
1260           gfc_error ("Result of BESSEL_YN is -INF at %L", &result->where);
1261           gfc_free_expr (result);
1262           return &gfc_bad_expr;
1263         }
1264
1265       if (jn && n1 == 0)
1266         {
1267           e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1268           mpfr_set_ui (e->value.real, 1, GFC_RND_MODE);
1269           gfc_constructor_append_expr (&result->value.constructor, e,
1270                                        &x->where);
1271           n1++;
1272         }
1273
1274       for (i = n1; i <= n2; i++)
1275         {
1276           e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1277           if (jn)
1278             mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
1279           else
1280             mpfr_set_inf (e->value.real, -1);
1281           gfc_constructor_append_expr (&result->value.constructor, e,
1282                                        &x->where);
1283         }
1284
1285       return result;
1286     }
1287
1288   /* Use the faster but more verbose recurrence algorithm. Bessel functions
1289      are stable for downward recursion and Neumann functions are stable
1290      for upward recursion. It is
1291        x2rev = 2.0/x,
1292        J(N-1, x) = x2rev * N * J(N, x) - J(N+1, x),
1293        Y(N+1, x) = x2rev * N * Y(N, x) - Y(N-1, x).
1294      Cf. http://dlmf.nist.gov/10.74#iv and http://dlmf.nist.gov/10.6#E1  */
1295
1296   gfc_set_model_kind (x->ts.kind);
1297
1298   /* Get first recursion anchor.  */
1299
1300   mpfr_init (last1);
1301   if (jn)
1302     mpfr_jn (last1, n2, x->value.real, GFC_RND_MODE);
1303   else
1304     mpfr_yn (last1, n1, x->value.real, GFC_RND_MODE);
1305
1306   e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1307   mpfr_set (e->value.real, last1, GFC_RND_MODE);
1308   if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
1309     {
1310       mpfr_clear (last1);
1311       gfc_free_expr (e);
1312       gfc_free_expr (result);
1313       return &gfc_bad_expr;
1314     }
1315   gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
1316
1317   if (n1 == n2)
1318     {
1319       mpfr_clear (last1);
1320       return result;
1321     }
1322  
1323   /* Get second recursion anchor.  */
1324
1325   mpfr_init (last2);
1326   if (jn)
1327     mpfr_jn (last2, n2-1, x->value.real, GFC_RND_MODE);
1328   else
1329     mpfr_yn (last2, n1+1, x->value.real, GFC_RND_MODE);
1330
1331   e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1332   mpfr_set (e->value.real, last2, GFC_RND_MODE);
1333   if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
1334     {
1335       mpfr_clear (last1);
1336       mpfr_clear (last2);
1337       gfc_free_expr (e);
1338       gfc_free_expr (result);
1339       return &gfc_bad_expr;
1340     }
1341   if (jn)
1342     gfc_constructor_insert_expr (&result->value.constructor, e, &x->where, -2);
1343   else 
1344     gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
1345
1346   if (n1 + 1 == n2)
1347     {
1348       mpfr_clear (last1);
1349       mpfr_clear (last2);
1350       return result;
1351     }
1352
1353   /* Start actual recursion.  */
1354
1355   mpfr_init (x2rev);
1356   mpfr_ui_div (x2rev, 2, x->value.real, GFC_RND_MODE);
1357  
1358   for (i = 2; i <= n2-n1; i++)
1359     {
1360       e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1361
1362       /* Special case: For YN, if the previous N gave -INF, set
1363          also N+1 to -INF.  */
1364       if (!jn && !gfc_option.flag_range_check && mpfr_inf_p (last2))
1365         {
1366           mpfr_set_inf (e->value.real, -1);
1367           gfc_constructor_append_expr (&result->value.constructor, e,
1368                                        &x->where);
1369           continue;
1370         }
1371
1372       mpfr_mul_si (e->value.real, x2rev, jn ? (n2-i+1) : (n1+i-1),
1373                    GFC_RND_MODE);
1374       mpfr_mul (e->value.real, e->value.real, last2, GFC_RND_MODE);
1375       mpfr_sub (e->value.real, e->value.real, last1, GFC_RND_MODE);
1376
1377       if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
1378         goto error;
1379
1380       if (jn)
1381         gfc_constructor_insert_expr (&result->value.constructor, e, &x->where,
1382                                      -i-1);
1383       else
1384         gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
1385
1386       mpfr_set (last1, last2, GFC_RND_MODE);
1387       mpfr_set (last2, e->value.real, GFC_RND_MODE);
1388     }
1389
1390   mpfr_clear (last1);
1391   mpfr_clear (last2);
1392   mpfr_clear (x2rev);
1393   return result;
1394
1395 error:
1396   mpfr_clear (last1);
1397   mpfr_clear (last2);
1398   mpfr_clear (x2rev);
1399   gfc_free_expr (e);
1400   gfc_free_expr (result);
1401   return &gfc_bad_expr;
1402 }
1403
1404
1405 gfc_expr *
1406 gfc_simplify_bessel_jn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
1407 {
1408   return gfc_simplify_bessel_n2 (order1, order2, x, true);
1409 }
1410
1411
1412 gfc_expr *
1413 gfc_simplify_bessel_y0 (gfc_expr *x)
1414 {
1415   gfc_expr *result;
1416
1417   if (x->expr_type != EXPR_CONSTANT)
1418     return NULL;
1419
1420   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1421   mpfr_y0 (result->value.real, x->value.real, GFC_RND_MODE);
1422
1423   return range_check (result, "BESSEL_Y0");
1424 }
1425
1426
1427 gfc_expr *
1428 gfc_simplify_bessel_y1 (gfc_expr *x)
1429 {
1430   gfc_expr *result;
1431
1432   if (x->expr_type != EXPR_CONSTANT)
1433     return NULL;
1434
1435   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1436   mpfr_y1 (result->value.real, x->value.real, GFC_RND_MODE);
1437
1438   return range_check (result, "BESSEL_Y1");
1439 }
1440
1441
1442 gfc_expr *
1443 gfc_simplify_bessel_yn (gfc_expr *order, gfc_expr *x)
1444 {
1445   gfc_expr *result;
1446   long n;
1447
1448   if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
1449     return NULL;
1450
1451   n = mpz_get_si (order->value.integer);
1452   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1453   mpfr_yn (result->value.real, n, x->value.real, GFC_RND_MODE);
1454
1455   return range_check (result, "BESSEL_YN");
1456 }
1457
1458
1459 gfc_expr *
1460 gfc_simplify_bessel_yn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
1461 {
1462   return gfc_simplify_bessel_n2 (order1, order2, x, false);
1463 }
1464
1465
1466 gfc_expr *
1467 gfc_simplify_bit_size (gfc_expr *e)
1468 {
1469   int i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
1470   return gfc_get_int_expr (e->ts.kind, &e->where,
1471                            gfc_integer_kinds[i].bit_size);
1472 }
1473
1474
1475 gfc_expr *
1476 gfc_simplify_btest (gfc_expr *e, gfc_expr *bit)
1477 {
1478   int b;
1479
1480   if (e->expr_type != EXPR_CONSTANT || bit->expr_type != EXPR_CONSTANT)
1481     return NULL;
1482
1483   if (gfc_extract_int (bit, &b) != NULL || b < 0)
1484     return gfc_get_logical_expr (gfc_default_logical_kind, &e->where, false);
1485
1486   return gfc_get_logical_expr (gfc_default_logical_kind, &e->where,
1487                                mpz_tstbit (e->value.integer, b));
1488 }
1489
1490
1491 static int
1492 compare_bitwise (gfc_expr *i, gfc_expr *j)
1493 {
1494   mpz_t x, y;
1495   int k, res;
1496
1497   gcc_assert (i->ts.type == BT_INTEGER);
1498   gcc_assert (j->ts.type == BT_INTEGER);
1499
1500   mpz_init_set (x, i->value.integer);
1501   k = gfc_validate_kind (i->ts.type, i->ts.kind, false);
1502   convert_mpz_to_unsigned (x, gfc_integer_kinds[k].bit_size);
1503
1504   mpz_init_set (y, j->value.integer);
1505   k = gfc_validate_kind (j->ts.type, j->ts.kind, false);
1506   convert_mpz_to_unsigned (y, gfc_integer_kinds[k].bit_size);
1507
1508   res = mpz_cmp (x, y);
1509   mpz_clear (x);
1510   mpz_clear (y);
1511   return res;
1512 }
1513
1514
1515 gfc_expr *
1516 gfc_simplify_bge (gfc_expr *i, gfc_expr *j)
1517 {
1518   if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1519     return NULL;
1520
1521   return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1522                                compare_bitwise (i, j) >= 0);
1523 }
1524
1525
1526 gfc_expr *
1527 gfc_simplify_bgt (gfc_expr *i, gfc_expr *j)
1528 {
1529   if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1530     return NULL;
1531
1532   return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1533                                compare_bitwise (i, j) > 0);
1534 }
1535
1536
1537 gfc_expr *
1538 gfc_simplify_ble (gfc_expr *i, gfc_expr *j)
1539 {
1540   if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1541     return NULL;
1542
1543   return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1544                                compare_bitwise (i, j) <= 0);
1545 }
1546
1547
1548 gfc_expr *
1549 gfc_simplify_blt (gfc_expr *i, gfc_expr *j)
1550 {
1551   if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1552     return NULL;
1553
1554   return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1555                                compare_bitwise (i, j) < 0);
1556 }
1557
1558
1559 gfc_expr *
1560 gfc_simplify_ceiling (gfc_expr *e, gfc_expr *k)
1561 {
1562   gfc_expr *ceil, *result;
1563   int kind;
1564
1565   kind = get_kind (BT_INTEGER, k, "CEILING", gfc_default_integer_kind);
1566   if (kind == -1)
1567     return &gfc_bad_expr;
1568
1569   if (e->expr_type != EXPR_CONSTANT)
1570     return NULL;
1571
1572   ceil = gfc_copy_expr (e);
1573   mpfr_ceil (ceil->value.real, e->value.real);
1574
1575   result = gfc_get_constant_expr (BT_INTEGER, kind, &e->where);
1576   gfc_mpfr_to_mpz (result->value.integer, ceil->value.real, &e->where);
1577
1578   gfc_free_expr (ceil);
1579
1580   return range_check (result, "CEILING");
1581 }
1582
1583
1584 gfc_expr *
1585 gfc_simplify_char (gfc_expr *e, gfc_expr *k)
1586 {
1587   return simplify_achar_char (e, k, "CHAR", false);
1588 }
1589
1590
1591 /* Common subroutine for simplifying CMPLX, COMPLEX and DCMPLX.  */
1592
1593 static gfc_expr *
1594 simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
1595 {
1596   gfc_expr *result;
1597
1598   if (convert_boz (x, kind) == &gfc_bad_expr)
1599     return &gfc_bad_expr;
1600
1601   if (convert_boz (y, kind) == &gfc_bad_expr)
1602     return &gfc_bad_expr;
1603
1604   if (x->expr_type != EXPR_CONSTANT
1605       || (y != NULL && y->expr_type != EXPR_CONSTANT))
1606     return NULL;
1607
1608   result = gfc_get_constant_expr (BT_COMPLEX, kind, &x->where);
1609
1610   switch (x->ts.type)
1611     {
1612       case BT_INTEGER:
1613         mpc_set_z (result->value.complex, x->value.integer, GFC_MPC_RND_MODE);
1614         break;
1615
1616       case BT_REAL:
1617         mpc_set_fr (result->value.complex, x->value.real, GFC_RND_MODE);
1618         break;
1619
1620       case BT_COMPLEX:
1621         mpc_set (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1622         break;
1623
1624       default:
1625         gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (x)");
1626     }
1627
1628   if (!y)
1629     return range_check (result, name);
1630
1631   switch (y->ts.type)
1632     {
1633       case BT_INTEGER:
1634         mpfr_set_z (mpc_imagref (result->value.complex),
1635                     y->value.integer, GFC_RND_MODE);
1636         break;
1637
1638       case BT_REAL:
1639         mpfr_set (mpc_imagref (result->value.complex),
1640                   y->value.real, GFC_RND_MODE);
1641         break;
1642
1643       default:
1644         gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (y)");
1645     }
1646
1647   return range_check (result, name);
1648 }
1649
1650
1651 gfc_expr *
1652 gfc_simplify_cmplx (gfc_expr *x, gfc_expr *y, gfc_expr *k)
1653 {
1654   int kind;
1655
1656   kind = get_kind (BT_REAL, k, "CMPLX", gfc_default_complex_kind);
1657   if (kind == -1)
1658     return &gfc_bad_expr;
1659
1660   return simplify_cmplx ("CMPLX", x, y, kind);
1661 }
1662
1663
1664 gfc_expr *
1665 gfc_simplify_complex (gfc_expr *x, gfc_expr *y)
1666 {
1667   int kind;
1668
1669   if (x->ts.type == BT_INTEGER && y->ts.type == BT_INTEGER)
1670     kind = gfc_default_complex_kind;
1671   else if (x->ts.type == BT_REAL || y->ts.type == BT_INTEGER)
1672     kind = x->ts.kind;
1673   else if (x->ts.type == BT_INTEGER || y->ts.type == BT_REAL)
1674     kind = y->ts.kind;
1675   else if (x->ts.type == BT_REAL && y->ts.type == BT_REAL)
1676     kind = (x->ts.kind > y->ts.kind) ? x->ts.kind : y->ts.kind;
1677   else
1678     gcc_unreachable ();
1679
1680   return simplify_cmplx ("COMPLEX", x, y, kind);
1681 }
1682
1683
1684 gfc_expr *
1685 gfc_simplify_conjg (gfc_expr *e)
1686 {
1687   gfc_expr *result;
1688
1689   if (e->expr_type != EXPR_CONSTANT)
1690     return NULL;
1691
1692   result = gfc_copy_expr (e);
1693   mpc_conj (result->value.complex, result->value.complex, GFC_MPC_RND_MODE);
1694
1695   return range_check (result, "CONJG");
1696 }
1697
1698
1699 gfc_expr *
1700 gfc_simplify_cos (gfc_expr *x)
1701 {
1702   gfc_expr *result;
1703
1704   if (x->expr_type != EXPR_CONSTANT)
1705     return NULL;
1706
1707   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1708
1709   switch (x->ts.type)
1710     {
1711       case BT_REAL:
1712         mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
1713         break;
1714
1715       case BT_COMPLEX:
1716         gfc_set_model_kind (x->ts.kind);
1717         mpc_cos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1718         break;
1719
1720       default:
1721         gfc_internal_error ("in gfc_simplify_cos(): Bad type");
1722     }
1723
1724   return range_check (result, "COS");
1725 }
1726
1727
1728 gfc_expr *
1729 gfc_simplify_cosh (gfc_expr *x)
1730 {
1731   gfc_expr *result;
1732
1733   if (x->expr_type != EXPR_CONSTANT)
1734     return NULL;
1735
1736   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1737
1738   switch (x->ts.type)
1739     {
1740       case BT_REAL:
1741         mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
1742         break;
1743
1744       case BT_COMPLEX:
1745         mpc_cosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1746         break;
1747         
1748       default:
1749         gcc_unreachable ();
1750     }
1751
1752   return range_check (result, "COSH");
1753 }
1754
1755
1756 gfc_expr *
1757 gfc_simplify_count (gfc_expr *mask, gfc_expr *dim, gfc_expr *kind)
1758 {
1759   gfc_expr *result;
1760
1761   if (!is_constant_array_expr (mask)
1762       || !gfc_is_constant_expr (dim)
1763       || !gfc_is_constant_expr (kind))
1764     return NULL;
1765
1766   result = transformational_result (mask, dim,
1767                                     BT_INTEGER,
1768                                     get_kind (BT_INTEGER, kind, "COUNT",
1769                                               gfc_default_integer_kind),
1770                                     &mask->where);
1771
1772   init_result_expr (result, 0, NULL);
1773
1774   /* Passing MASK twice, once as data array, once as mask.
1775      Whenever gfc_count is called, '1' is added to the result.  */
1776   return !dim || mask->rank == 1 ?
1777     simplify_transformation_to_scalar (result, mask, mask, gfc_count) :
1778     simplify_transformation_to_array (result, mask, dim, mask, gfc_count, NULL);
1779 }
1780
1781
1782 gfc_expr *
1783 gfc_simplify_dcmplx (gfc_expr *x, gfc_expr *y)
1784 {
1785   return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind);
1786 }
1787
1788
1789 gfc_expr *
1790 gfc_simplify_dble (gfc_expr *e)
1791 {
1792   gfc_expr *result = NULL;
1793
1794   if (e->expr_type != EXPR_CONSTANT)
1795     return NULL;
1796
1797   if (convert_boz (e, gfc_default_double_kind) == &gfc_bad_expr)
1798     return &gfc_bad_expr;
1799
1800   result = gfc_convert_constant (e, BT_REAL, gfc_default_double_kind);
1801   if (result == &gfc_bad_expr)
1802     return &gfc_bad_expr;
1803
1804   return range_check (result, "DBLE");
1805 }
1806
1807
1808 gfc_expr *
1809 gfc_simplify_digits (gfc_expr *x)
1810 {
1811   int i, digits;
1812
1813   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
1814
1815   switch (x->ts.type)
1816     {
1817       case BT_INTEGER:
1818         digits = gfc_integer_kinds[i].digits;
1819         break;
1820
1821       case BT_REAL:
1822       case BT_COMPLEX:
1823         digits = gfc_real_kinds[i].digits;
1824         break;
1825
1826       default:
1827         gcc_unreachable ();
1828     }
1829
1830   return gfc_get_int_expr (gfc_default_integer_kind, NULL, digits);
1831 }
1832
1833
1834 gfc_expr *
1835 gfc_simplify_dim (gfc_expr *x, gfc_expr *y)
1836 {
1837   gfc_expr *result;
1838   int kind;
1839
1840   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1841     return NULL;
1842
1843   kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
1844   result = gfc_get_constant_expr (x->ts.type, kind, &x->where);
1845
1846   switch (x->ts.type)
1847     {
1848       case BT_INTEGER:
1849         if (mpz_cmp (x->value.integer, y->value.integer) > 0)
1850           mpz_sub (result->value.integer, x->value.integer, y->value.integer);
1851         else
1852           mpz_set_ui (result->value.integer, 0);
1853
1854         break;
1855
1856       case BT_REAL:
1857         if (mpfr_cmp (x->value.real, y->value.real) > 0)
1858           mpfr_sub (result->value.real, x->value.real, y->value.real,
1859                     GFC_RND_MODE);
1860         else
1861           mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
1862
1863         break;
1864
1865       default:
1866         gfc_internal_error ("gfc_simplify_dim(): Bad type");
1867     }
1868
1869   return range_check (result, "DIM");
1870 }
1871
1872
1873 gfc_expr*
1874 gfc_simplify_dot_product (gfc_expr *vector_a, gfc_expr *vector_b)
1875 {
1876   if (!is_constant_array_expr (vector_a)
1877       || !is_constant_array_expr (vector_b))
1878     return NULL;
1879
1880   gcc_assert (vector_a->rank == 1);
1881   gcc_assert (vector_b->rank == 1);
1882   gcc_assert (gfc_compare_types (&vector_a->ts, &vector_b->ts));
1883
1884   return compute_dot_product (vector_a, 1, 0, vector_b, 1, 0, true);
1885 }
1886
1887
1888 gfc_expr *
1889 gfc_simplify_dprod (gfc_expr *x, gfc_expr *y)
1890 {
1891   gfc_expr *a1, *a2, *result;
1892
1893   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1894     return NULL;
1895
1896   a1 = gfc_real2real (x, gfc_default_double_kind);
1897   a2 = gfc_real2real (y, gfc_default_double_kind);
1898
1899   result = gfc_get_constant_expr (BT_REAL, gfc_default_double_kind, &x->where);
1900   mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
1901
1902   gfc_free_expr (a2);
1903   gfc_free_expr (a1);
1904
1905   return range_check (result, "DPROD");
1906 }
1907
1908
1909 static gfc_expr *
1910 simplify_dshift (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg,
1911                       bool right)
1912 {
1913   gfc_expr *result;
1914   int i, k, size, shift;
1915
1916   if (arg1->expr_type != EXPR_CONSTANT || arg2->expr_type != EXPR_CONSTANT
1917       || shiftarg->expr_type != EXPR_CONSTANT)
1918     return NULL;
1919
1920   k = gfc_validate_kind (BT_INTEGER, arg1->ts.kind, false);
1921   size = gfc_integer_kinds[k].bit_size;
1922
1923   gfc_extract_int (shiftarg, &shift);
1924
1925   /* DSHIFTR(I,J,SHIFT) = DSHIFTL(I,J,SIZE-SHIFT).  */
1926   if (right)
1927     shift = size - shift;
1928
1929   result = gfc_get_constant_expr (BT_INTEGER, arg1->ts.kind, &arg1->where);
1930   mpz_set_ui (result->value.integer, 0);
1931
1932   for (i = 0; i < shift; i++)
1933     if (mpz_tstbit (arg2->value.integer, size - shift + i))
1934       mpz_setbit (result->value.integer, i);
1935
1936   for (i = 0; i < size - shift; i++)
1937     if (mpz_tstbit (arg1->value.integer, i))
1938       mpz_setbit (result->value.integer, shift + i);
1939
1940   /* Convert to a signed value.  */
1941   convert_mpz_to_signed (result->value.integer, size);
1942
1943   return result;
1944 }
1945
1946
1947 gfc_expr *
1948 gfc_simplify_dshiftr (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg)
1949 {
1950   return simplify_dshift (arg1, arg2, shiftarg, true);
1951 }
1952
1953
1954 gfc_expr *
1955 gfc_simplify_dshiftl (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg)
1956 {
1957   return simplify_dshift (arg1, arg2, shiftarg, false);
1958 }
1959
1960
1961 gfc_expr *
1962 gfc_simplify_erf (gfc_expr *x)
1963 {
1964   gfc_expr *result;
1965
1966   if (x->expr_type != EXPR_CONSTANT)
1967     return NULL;
1968
1969   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1970   mpfr_erf (result->value.real, x->value.real, GFC_RND_MODE);
1971
1972   return range_check (result, "ERF");
1973 }
1974
1975
1976 gfc_expr *
1977 gfc_simplify_erfc (gfc_expr *x)
1978 {
1979   gfc_expr *result;
1980
1981   if (x->expr_type != EXPR_CONSTANT)
1982     return NULL;
1983
1984   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1985   mpfr_erfc (result->value.real, x->value.real, GFC_RND_MODE);
1986
1987   return range_check (result, "ERFC");
1988 }
1989
1990
1991 /* Helper functions to simplify ERFC_SCALED(x) = ERFC(x) * EXP(X**2).  */
1992
1993 #define MAX_ITER 200
1994 #define ARG_LIMIT 12
1995
1996 /* Calculate ERFC_SCALED directly by its definition:
1997
1998      ERFC_SCALED(x) = ERFC(x) * EXP(X**2)
1999
2000    using a large precision for intermediate results.  This is used for all
2001    but large values of the argument.  */
2002 static void
2003 fullprec_erfc_scaled (mpfr_t res, mpfr_t arg)
2004 {
2005   mp_prec_t prec;
2006   mpfr_t a, b;
2007
2008   prec = mpfr_get_default_prec ();
2009   mpfr_set_default_prec (10 * prec);
2010
2011   mpfr_init (a);
2012   mpfr_init (b);
2013
2014   mpfr_set (a, arg, GFC_RND_MODE);
2015   mpfr_sqr (b, a, GFC_RND_MODE);
2016   mpfr_exp (b, b, GFC_RND_MODE);
2017   mpfr_erfc (a, a, GFC_RND_MODE);
2018   mpfr_mul (a, a, b, GFC_RND_MODE);
2019
2020   mpfr_set (res, a, GFC_RND_MODE);
2021   mpfr_set_default_prec (prec);
2022
2023   mpfr_clear (a);
2024   mpfr_clear (b);
2025 }
2026
2027 /* Calculate ERFC_SCALED using a power series expansion in 1/arg:
2028
2029     ERFC_SCALED(x) = 1 / (x * sqrt(pi))
2030                      * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
2031                                           / (2 * x**2)**n)
2032
2033   This is used for large values of the argument.  Intermediate calculations
2034   are performed with twice the precision.  We don't do a fixed number of
2035   iterations of the sum, but stop when it has converged to the required
2036   precision.  */
2037 static void
2038 asympt_erfc_scaled (mpfr_t res, mpfr_t arg)
2039 {
2040   mpfr_t sum, x, u, v, w, oldsum, sumtrunc;
2041   mpz_t num;
2042   mp_prec_t prec;
2043   unsigned i;
2044
2045   prec = mpfr_get_default_prec ();
2046   mpfr_set_default_prec (2 * prec);
2047
2048   mpfr_init (sum);
2049   mpfr_init (x);
2050   mpfr_init (u);
2051   mpfr_init (v);
2052   mpfr_init (w);
2053   mpz_init (num);
2054
2055   mpfr_init (oldsum);
2056   mpfr_init (sumtrunc);
2057   mpfr_set_prec (oldsum, prec);
2058   mpfr_set_prec (sumtrunc, prec);
2059
2060   mpfr_set (x, arg, GFC_RND_MODE);
2061   mpfr_set_ui (sum, 1, GFC_RND_MODE);
2062   mpz_set_ui (num, 1);
2063
2064   mpfr_set (u, x, GFC_RND_MODE);
2065   mpfr_sqr (u, u, GFC_RND_MODE);
2066   mpfr_mul_ui (u, u, 2, GFC_RND_MODE);
2067   mpfr_pow_si (u, u, -1, GFC_RND_MODE);
2068
2069   for (i = 1; i < MAX_ITER; i++)
2070   {
2071     mpfr_set (oldsum, sum, GFC_RND_MODE);
2072
2073     mpz_mul_ui (num, num, 2 * i - 1);
2074     mpz_neg (num, num);
2075
2076     mpfr_set (w, u, GFC_RND_MODE);
2077     mpfr_pow_ui (w, w, i, GFC_RND_MODE);
2078
2079     mpfr_set_z (v, num, GFC_RND_MODE);
2080     mpfr_mul (v, v, w, GFC_RND_MODE);
2081
2082     mpfr_add (sum, sum, v, GFC_RND_MODE);
2083
2084     mpfr_set (sumtrunc, sum, GFC_RND_MODE);
2085     if (mpfr_cmp (sumtrunc, oldsum) == 0)
2086       break;
2087   }
2088
2089   /* We should have converged by now; otherwise, ARG_LIMIT is probably
2090      set too low.  */
2091   gcc_assert (i < MAX_ITER);
2092
2093   /* Divide by x * sqrt(Pi).  */
2094   mpfr_const_pi (u, GFC_RND_MODE);
2095   mpfr_sqrt (u, u, GFC_RND_MODE);
2096   mpfr_mul (u, u, x, GFC_RND_MODE);
2097   mpfr_div (sum, sum, u, GFC_RND_MODE);
2098
2099   mpfr_set (res, sum, GFC_RND_MODE);
2100   mpfr_set_default_prec (prec);
2101
2102   mpfr_clears (sum, x, u, v, w, oldsum, sumtrunc, NULL);
2103   mpz_clear (num);
2104 }
2105
2106
2107 gfc_expr *
2108 gfc_simplify_erfc_scaled (gfc_expr *x)
2109 {
2110   gfc_expr *result;
2111
2112   if (x->expr_type != EXPR_CONSTANT)
2113     return NULL;
2114
2115   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2116   if (mpfr_cmp_d (x->value.real, ARG_LIMIT) >= 0)
2117     asympt_erfc_scaled (result->value.real, x->value.real);
2118   else
2119     fullprec_erfc_scaled (result->value.real, x->value.real);
2120
2121   return range_check (result, "ERFC_SCALED");
2122 }
2123
2124 #undef MAX_ITER
2125 #undef ARG_LIMIT
2126
2127
2128 gfc_expr *
2129 gfc_simplify_epsilon (gfc_expr *e)
2130 {
2131   gfc_expr *result;
2132   int i;
2133
2134   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2135
2136   result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
2137   mpfr_set (result->value.real, gfc_real_kinds[i].epsilon, GFC_RND_MODE);
2138
2139   return range_check (result, "EPSILON");
2140 }
2141
2142
2143 gfc_expr *
2144 gfc_simplify_exp (gfc_expr *x)
2145 {
2146   gfc_expr *result;
2147
2148   if (x->expr_type != EXPR_CONSTANT)
2149     return NULL;
2150
2151   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2152
2153   switch (x->ts.type)
2154     {
2155       case BT_REAL:
2156         mpfr_exp (result->value.real, x->value.real, GFC_RND_MODE);
2157         break;
2158
2159       case BT_COMPLEX:
2160         gfc_set_model_kind (x->ts.kind);
2161         mpc_exp (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
2162         break;
2163
2164       default:
2165         gfc_internal_error ("in gfc_simplify_exp(): Bad type");
2166     }
2167
2168   return range_check (result, "EXP");
2169 }
2170
2171
2172 gfc_expr *
2173 gfc_simplify_exponent (gfc_expr *x)
2174 {
2175   int i;
2176   gfc_expr *result;
2177
2178   if (x->expr_type != EXPR_CONSTANT)
2179     return NULL;
2180
2181   result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
2182                                   &x->where);
2183
2184   gfc_set_model (x->value.real);
2185
2186   if (mpfr_sgn (x->value.real) == 0)
2187     {
2188       mpz_set_ui (result->value.integer, 0);
2189       return result;
2190     }
2191
2192   i = (int) mpfr_get_exp (x->value.real);
2193   mpz_set_si (result->value.integer, i);
2194
2195   return range_check (result, "EXPONENT");
2196 }
2197
2198
2199 gfc_expr *
2200 gfc_simplify_float (gfc_expr *a)
2201 {
2202   gfc_expr *result;
2203
2204   if (a->expr_type != EXPR_CONSTANT)
2205     return NULL;
2206
2207   if (a->is_boz)
2208     {
2209       if (convert_boz (a, gfc_default_real_kind) == &gfc_bad_expr)
2210         return &gfc_bad_expr;
2211
2212       result = gfc_copy_expr (a);
2213     }
2214   else
2215     result = gfc_int2real (a, gfc_default_real_kind);
2216
2217   return range_check (result, "FLOAT");
2218 }
2219
2220
2221 static bool
2222 is_last_ref_vtab (gfc_expr *e)
2223 {
2224   gfc_ref *ref;
2225   gfc_component *comp = NULL;
2226
2227   if (e->expr_type != EXPR_VARIABLE)
2228     return false;
2229
2230   for (ref = e->ref; ref; ref = ref->next)
2231     if (ref->type == REF_COMPONENT)
2232       comp = ref->u.c.component;
2233
2234   if (!e->ref || !comp)
2235     return e->symtree->n.sym->attr.vtab;
2236
2237   if (comp->name[0] == '_' && strcmp (comp->name, "_vptr") == 0)
2238     return true;
2239
2240   return false;
2241 }
2242
2243
2244 gfc_expr *
2245 gfc_simplify_extends_type_of (gfc_expr *a, gfc_expr *mold)
2246 {
2247   /* Avoid simplification of resolved symbols.  */
2248   if (is_last_ref_vtab (a) || is_last_ref_vtab (mold))
2249     return NULL;
2250
2251   if (a->ts.type == BT_DERIVED && mold->ts.type == BT_DERIVED)
2252     return gfc_get_logical_expr (gfc_default_logical_kind, &a->where,
2253                                  gfc_type_is_extension_of (mold->ts.u.derived,
2254                                                            a->ts.u.derived));
2255   /* Return .false. if the dynamic type can never be the same.  */
2256   if ((a->ts.type == BT_CLASS && mold->ts.type == BT_CLASS
2257        && !gfc_type_is_extension_of
2258                         (mold->ts.u.derived->components->ts.u.derived,
2259                          a->ts.u.derived->components->ts.u.derived)
2260        && !gfc_type_is_extension_of
2261                         (a->ts.u.derived->components->ts.u.derived,
2262                          mold->ts.u.derived->components->ts.u.derived))
2263       || (a->ts.type == BT_DERIVED && mold->ts.type == BT_CLASS
2264           && !gfc_type_is_extension_of
2265                         (a->ts.u.derived,
2266                          mold->ts.u.derived->components->ts.u.derived)
2267           && !gfc_type_is_extension_of
2268                         (mold->ts.u.derived->components->ts.u.derived,
2269                          a->ts.u.derived))
2270       || (a->ts.type == BT_CLASS && mold->ts.type == BT_DERIVED
2271           && !gfc_type_is_extension_of
2272                         (mold->ts.u.derived,
2273                          a->ts.u.derived->components->ts.u.derived)))
2274     return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, false);
2275
2276   if (mold->ts.type == BT_DERIVED
2277       && gfc_type_is_extension_of (mold->ts.u.derived,
2278                                    a->ts.u.derived->components->ts.u.derived))
2279     return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, true);
2280
2281   return NULL;
2282 }
2283
2284
2285 gfc_expr *
2286 gfc_simplify_same_type_as (gfc_expr *a, gfc_expr *b)
2287 {
2288   /* Avoid simplification of resolved symbols.  */
2289   if (is_last_ref_vtab (a) || is_last_ref_vtab (b))
2290     return NULL;
2291
2292   /* Return .false. if the dynamic type can never be the
2293      same.  */
2294   if ((a->ts.type == BT_CLASS || b->ts.type == BT_CLASS)
2295       && !gfc_type_compatible (&a->ts, &b->ts)
2296       && !gfc_type_compatible (&b->ts, &a->ts))
2297     return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, false);
2298
2299   if (a->ts.type != BT_DERIVED || b->ts.type != BT_DERIVED)
2300      return NULL;
2301
2302   return gfc_get_logical_expr (gfc_default_logical_kind, &a->where,
2303                                gfc_compare_derived_types (a->ts.u.derived,
2304                                                           b->ts.u.derived));
2305 }
2306
2307
2308 gfc_expr *
2309 gfc_simplify_floor (gfc_expr *e, gfc_expr *k)
2310 {
2311   gfc_expr *result;
2312   mpfr_t floor;
2313   int kind;
2314
2315   kind = get_kind (BT_INTEGER, k, "FLOOR", gfc_default_integer_kind);
2316   if (kind == -1)
2317     gfc_internal_error ("gfc_simplify_floor(): Bad kind");
2318
2319   if (e->expr_type != EXPR_CONSTANT)
2320     return NULL;
2321
2322   gfc_set_model_kind (kind);
2323
2324   mpfr_init (floor);
2325   mpfr_floor (floor, e->value.real);
2326
2327   result = gfc_get_constant_expr (BT_INTEGER, kind, &e->where);
2328   gfc_mpfr_to_mpz (result->value.integer, floor, &e->where);
2329
2330   mpfr_clear (floor);
2331
2332   return range_check (result, "FLOOR");
2333 }
2334
2335
2336 gfc_expr *
2337 gfc_simplify_fraction (gfc_expr *x)
2338 {
2339   gfc_expr *result;
2340   mpfr_t absv, exp, pow2;
2341
2342   if (x->expr_type != EXPR_CONSTANT)
2343     return NULL;
2344
2345   result = gfc_get_constant_expr (BT_REAL, x->ts.kind, &x->where);
2346
2347   if (mpfr_sgn (x->value.real) == 0)
2348     {
2349       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2350       return result;
2351     }
2352
2353   gfc_set_model_kind (x->ts.kind);
2354   mpfr_init (exp);
2355   mpfr_init (absv);
2356   mpfr_init (pow2);
2357
2358   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
2359   mpfr_log2 (exp, absv, GFC_RND_MODE);
2360
2361   mpfr_trunc (exp, exp);
2362   mpfr_add_ui (exp, exp, 1, GFC_RND_MODE);
2363
2364   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
2365
2366   mpfr_div (result->value.real, absv, pow2, GFC_RND_MODE);
2367
2368   mpfr_clears (exp, absv, pow2, NULL);
2369
2370   return range_check (result, "FRACTION");
2371 }
2372
2373
2374 gfc_expr *
2375 gfc_simplify_gamma (gfc_expr *x)
2376 {
2377   gfc_expr *result;
2378
2379   if (x->expr_type != EXPR_CONSTANT)
2380     return NULL;
2381
2382   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2383   mpfr_gamma (result->value.real, x->value.real, GFC_RND_MODE);
2384
2385   return range_check (result, "GAMMA");
2386 }
2387
2388
2389 gfc_expr *
2390 gfc_simplify_huge (gfc_expr *e)
2391 {
2392   gfc_expr *result;
2393   int i;
2394
2395   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2396   result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);
2397
2398   switch (e->ts.type)
2399     {
2400       case BT_INTEGER:
2401         mpz_set (result->value.integer, gfc_integer_kinds[i].huge);
2402         break;
2403
2404       case BT_REAL:
2405         mpfr_set (result->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
2406         break;
2407
2408       default:
2409         gcc_unreachable ();
2410     }
2411
2412   return result;
2413 }
2414
2415
2416 gfc_expr *
2417 gfc_simplify_hypot (gfc_expr *x, gfc_expr *y)
2418 {
2419   gfc_expr *result;
2420
2421   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2422     return NULL;
2423
2424   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2425   mpfr_hypot (result->value.real, x->value.real, y->value.real, GFC_RND_MODE);
2426   return range_check (result, "HYPOT");
2427 }
2428
2429
2430 /* We use the processor's collating sequence, because all
2431    systems that gfortran currently works on are ASCII.  */
2432
2433 gfc_expr *
2434 gfc_simplify_iachar (gfc_expr *e, gfc_expr *kind)
2435 {
2436   gfc_expr *result;
2437   gfc_char_t index;
2438   int k;
2439
2440   if (e->expr_type != EXPR_CONSTANT)
2441     return NULL;
2442
2443   if (e->value.character.length != 1)
2444     {
2445       gfc_error ("Argument of IACHAR at %L must be of length one", &e->where);
2446       return &gfc_bad_expr;
2447     }
2448
2449   index = e->value.character.string[0];
2450
2451   if (gfc_option.warn_surprising && index > 127)
2452     gfc_warning ("Argument of IACHAR function at %L outside of range 0..127",
2453                  &e->where);
2454
2455   k = get_kind (BT_INTEGER, kind, "IACHAR", gfc_default_integer_kind);
2456   if (k == -1)
2457     return &gfc_bad_expr;
2458
2459   result = gfc_get_int_expr (k, &e->where, index);
2460
2461   return range_check (result, "IACHAR");
2462 }
2463
2464
2465 static gfc_expr *
2466 do_bit_and (gfc_expr *result, gfc_expr *e)
2467 {
2468   gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
2469   gcc_assert (result->ts.type == BT_INTEGER
2470               && result->expr_type == EXPR_CONSTANT);
2471
2472   mpz_and (result->value.integer, result->value.integer, e->value.integer);
2473   return result;
2474 }
2475
2476
2477 gfc_expr *
2478 gfc_simplify_iall (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
2479 {
2480   return simplify_transformation (array, dim, mask, -1, do_bit_and);
2481 }
2482
2483
2484 static gfc_expr *
2485 do_bit_ior (gfc_expr *result, gfc_expr *e)
2486 {
2487   gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
2488   gcc_assert (result->ts.type == BT_INTEGER
2489               && result->expr_type == EXPR_CONSTANT);
2490
2491   mpz_ior (result->value.integer, result->value.integer, e->value.integer);
2492   return result;
2493 }
2494
2495
2496 gfc_expr *
2497 gfc_simplify_iany (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
2498 {
2499   return simplify_transformation (array, dim, mask, 0, do_bit_ior);
2500 }
2501
2502
2503 gfc_expr *
2504 gfc_simplify_iand (gfc_expr *x, gfc_expr *y)
2505 {
2506   gfc_expr *result;
2507
2508   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2509     return NULL;
2510
2511   result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
2512   mpz_and (result->value.integer, x->value.integer, y->value.integer);
2513
2514   return range_check (result, "IAND");
2515 }
2516
2517
2518 gfc_expr *
2519 gfc_simplify_ibclr (gfc_expr *x, gfc_expr *y)
2520 {
2521   gfc_expr *result;
2522   int k, pos;
2523
2524   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2525     return NULL;
2526
2527   gfc_extract_int (y, &pos);
2528
2529   k = gfc_validate_kind (x->ts.type, x->ts.kind, false);
2530
2531   result = gfc_copy_expr (x);
2532
2533   convert_mpz_to_unsigned (result->value.integer,
2534                            gfc_integer_kinds[k].bit_size);
2535
2536   mpz_clrbit (result->value.integer, pos);
2537
2538   convert_mpz_to_signed (result->value.integer,
2539                          gfc_integer_kinds[k].bit_size);
2540
2541   return result;
2542 }
2543
2544
2545 gfc_expr *
2546 gfc_simplify_ibits (gfc_expr *x, gfc_expr *y, gfc_expr *z)
2547 {
2548   gfc_expr *result;
2549   int pos, len;
2550   int i, k, bitsize;
2551   int *bits;
2552
2553   if (x->expr_type != EXPR_CONSTANT
2554       || y->expr_type != EXPR_CONSTANT
2555       || z->expr_type != EXPR_CONSTANT)
2556     return NULL;
2557
2558   gfc_extract_int (y, &pos);
2559   gfc_extract_int (z, &len);
2560
2561   k = gfc_validate_kind (BT_INTEGER, x->ts.kind, false);
2562
2563   bitsize = gfc_integer_kinds[k].bit_size;
2564
2565   if (pos + len > bitsize)
2566     {
2567       gfc_error ("Sum of second and third arguments of IBITS exceeds "
2568                  "bit size at %L", &y->where);
2569       return &gfc_bad_expr;
2570     }
2571
2572   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2573   convert_mpz_to_unsigned (result->value.integer,
2574                            gfc_integer_kinds[k].bit_size);
2575
2576   bits = XCNEWVEC (int, bitsize);
2577
2578   for (i = 0; i < bitsize; i++)
2579     bits[i] = 0;
2580
2581   for (i = 0; i < len; i++)
2582     bits[i] = mpz_tstbit (x->value.integer, i + pos);
2583
2584   for (i = 0; i < bitsize; i++)
2585     {
2586       if (bits[i] == 0)
2587         mpz_clrbit (result->value.integer, i);
2588       else if (bits[i] == 1)
2589         mpz_setbit (result->value.integer, i);
2590       else
2591         gfc_internal_error ("IBITS: Bad bit");
2592     }
2593
2594   free (bits);
2595
2596   convert_mpz_to_signed (result->value.integer,
2597                          gfc_integer_kinds[k].bit_size);
2598
2599   return result;
2600 }
2601
2602
2603 gfc_expr *
2604 gfc_simplify_ibset (gfc_expr *x, gfc_expr *y)
2605 {
2606   gfc_expr *result;
2607   int k, pos;
2608
2609   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2610     return NULL;
2611
2612   gfc_extract_int (y, &pos);
2613
2614   k = gfc_validate_kind (x->ts.type, x->ts.kind, false);
2615
2616   result = gfc_copy_expr (x);
2617
2618   convert_mpz_to_unsigned (result->value.integer,
2619                            gfc_integer_kinds[k].bit_size);
2620
2621   mpz_setbit (result->value.integer, pos);
2622
2623   convert_mpz_to_signed (result->value.integer,
2624                          gfc_integer_kinds[k].bit_size);
2625
2626   return result;
2627 }
2628
2629
2630 gfc_expr *
2631 gfc_simplify_ichar (gfc_expr *e, gfc_expr *kind)
2632 {
2633   gfc_expr *result;
2634   gfc_char_t index;
2635   int k;
2636
2637   if (e->expr_type != EXPR_CONSTANT)
2638     return NULL;
2639
2640   if (e->value.character.length != 1)
2641     {
2642       gfc_error ("Argument of ICHAR at %L must be of length one", &e->where);
2643       return &gfc_bad_expr;
2644     }
2645
2646   index = e->value.character.string[0];
2647
2648   k = get_kind (BT_INTEGER, kind, "ICHAR", gfc_default_integer_kind);
2649   if (k == -1)
2650     return &gfc_bad_expr;
2651
2652   result = gfc_get_int_expr (k, &e->where, index);
2653
2654   return range_check (result, "ICHAR");
2655 }
2656
2657
2658 gfc_expr *
2659 gfc_simplify_ieor (gfc_expr *x, gfc_expr *y)
2660 {
2661   gfc_expr *result;
2662
2663   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2664     return NULL;
2665
2666   result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
2667   mpz_xor (result->value.integer, x->value.integer, y->value.integer);
2668
2669   return range_check (result, "IEOR");
2670 }
2671
2672
2673 gfc_expr *
2674 gfc_simplify_index (gfc_expr *x, gfc_expr *y, gfc_expr *b, gfc_expr *kind)
2675 {
2676   gfc_expr *result;
2677   int back, len, lensub;
2678   int i, j, k, count, index = 0, start;
2679
2680   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT 
2681       || ( b != NULL && b->expr_type !=  EXPR_CONSTANT))
2682     return NULL;
2683
2684   if (b != NULL && b->value.logical != 0)
2685     back = 1;
2686   else
2687     back = 0;
2688
2689   k = get_kind (BT_INTEGER, kind, "INDEX", gfc_default_integer_kind); 
2690   if (k == -1)
2691     return &gfc_bad_expr;
2692
2693   result = gfc_get_constant_expr (BT_INTEGER, k, &x->where);
2694
2695   len = x->value.character.length;
2696   lensub = y->value.character.length;
2697
2698   if (len < lensub)
2699     {
2700       mpz_set_si (result->value.integer, 0);
2701       return result;
2702     }
2703
2704   if (back == 0)
2705     {
2706       if (lensub == 0)
2707         {
2708           mpz_set_si (result->value.integer, 1);
2709           return result;
2710         }
2711       else if (lensub == 1)
2712         {
2713           for (i = 0; i < len; i++)
2714             {
2715               for (j = 0; j < lensub; j++)
2716                 {
2717                   if (y->value.character.string[j]
2718                       == x->value.character.string[i])
2719                     {
2720                       index = i + 1;
2721                       goto done;
2722                     }
2723                 }
2724             }
2725         }
2726       else
2727         {
2728           for (i = 0; i < len; i++)
2729             {
2730               for (j = 0; j < lensub; j++)
2731                 {
2732                   if (y->value.character.string[j]
2733                       == x->value.character.string[i])
2734                     {
2735                       start = i;
2736                       count = 0;
2737
2738                       for (k = 0; k < lensub; k++)
2739                         {
2740                           if (y->value.character.string[k]
2741                               == x->value.character.string[k + start])
2742                             count++;
2743                         }
2744
2745                       if (count == lensub)
2746                         {
2747                           index = start + 1;
2748                           goto done;
2749                         }
2750                     }
2751                 }
2752             }
2753         }
2754
2755     }
2756   else
2757     {
2758       if (lensub == 0)
2759         {
2760           mpz_set_si (result->value.integer, len + 1);
2761           return result;
2762         }
2763       else if (lensub == 1)
2764         {
2765           for (i = 0; i < len; i++)
2766             {
2767               for (j = 0; j < lensub; j++)
2768                 {
2769                   if (y->value.character.string[j]
2770                       == x->value.character.string[len - i])
2771                     {
2772                       index = len - i + 1;
2773                       goto done;
2774                     }
2775                 }
2776             }
2777         }
2778       else
2779         {
2780           for (i = 0; i < len; i++)
2781             {
2782               for (j = 0; j < lensub; j++)
2783                 {
2784                   if (y->value.character.string[j]
2785                       == x->value.character.string[len - i])
2786                     {
2787                       start = len - i;
2788                       if (start <= len - lensub)
2789                         {
2790                           count = 0;
2791                           for (k = 0; k < lensub; k++)
2792                             if (y->value.character.string[k]
2793                                 == x->value.character.string[k + start])
2794                               count++;
2795
2796                           if (count == lensub)
2797                             {
2798                               index = start + 1;
2799                               goto done;
2800                             }
2801                         }
2802                       else
2803                         {
2804                           continue;
2805                         }
2806                     }
2807                 }
2808             }
2809         }
2810     }
2811
2812 done:
2813   mpz_set_si (result->value.integer, index);
2814   return range_check (result, "INDEX");
2815 }
2816
2817
2818 static gfc_expr *
2819 simplify_intconv (gfc_expr *e, int kind, const char *name)
2820 {
2821   gfc_expr *result = NULL;
2822
2823   if (e->expr_type != EXPR_CONSTANT)
2824     return NULL;
2825
2826   result = gfc_convert_constant (e, BT_INTEGER, kind);
2827   if (result == &gfc_bad_expr)
2828     return &gfc_bad_expr;
2829
2830   return range_check (result, name);
2831 }
2832
2833
2834 gfc_expr *
2835 gfc_simplify_int (gfc_expr *e, gfc_expr *k)
2836 {
2837   int kind;
2838
2839   kind = get_kind (BT_INTEGER, k, "INT", gfc_default_integer_kind);
2840   if (kind == -1)
2841     return &gfc_bad_expr;
2842
2843   return simplify_intconv (e, kind, "INT");
2844 }
2845
2846 gfc_expr *
2847 gfc_simplify_int2 (gfc_expr *e)
2848 {
2849   return simplify_intconv (e, 2, "INT2");
2850 }
2851
2852
2853 gfc_expr *
2854 gfc_simplify_int8 (gfc_expr *e)
2855 {
2856   return simplify_intconv (e, 8, "INT8");
2857 }
2858
2859
2860 gfc_expr *
2861 gfc_simplify_long (gfc_expr *e)
2862 {
2863   return simplify_intconv (e, 4, "LONG");
2864 }
2865
2866
2867 gfc_expr *
2868 gfc_simplify_ifix (gfc_expr *e)
2869 {
2870   gfc_expr *rtrunc, *result;
2871
2872   if (e->expr_type != EXPR_CONSTANT)
2873     return NULL;
2874
2875   rtrunc = gfc_copy_expr (e);
2876   mpfr_trunc (rtrunc->value.real, e->value.real);
2877
2878   result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
2879                                   &e->where);
2880   gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real, &e->where);
2881
2882   gfc_free_expr (rtrunc);
2883
2884   return range_check (result, "IFIX");
2885 }
2886
2887
2888 gfc_expr *
2889 gfc_simplify_idint (gfc_expr *e)
2890 {
2891   gfc_expr *rtrunc, *result;
2892
2893   if (e->expr_type != EXPR_CONSTANT)
2894     return NULL;
2895
2896   rtrunc = gfc_copy_expr (e);
2897   mpfr_trunc (rtrunc->value.real, e->value.real);
2898
2899   result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
2900                                   &e->where);
2901   gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real, &e->where);
2902
2903   gfc_free_expr (rtrunc);
2904
2905   return range_check (result, "IDINT");
2906 }
2907
2908
2909 gfc_expr *
2910 gfc_simplify_ior (gfc_expr *x, gfc_expr *y)
2911 {
2912   gfc_expr *result;
2913
2914   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2915     return NULL;
2916
2917   result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
2918   mpz_ior (result->value.integer, x->value.integer, y->value.integer);
2919
2920   return range_check (result, "IOR");
2921 }
2922
2923
2924 static gfc_expr *
2925 do_bit_xor (gfc_expr *result, gfc_expr *e)
2926 {
2927   gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
2928   gcc_assert (result->ts.type == BT_INTEGER
2929               && result->expr_type == EXPR_CONSTANT);
2930
2931   mpz_xor (result->value.integer, result->value.integer, e->value.integer);
2932   return result;
2933 }
2934
2935
2936 gfc_expr *
2937 gfc_simplify_iparity (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
2938 {
2939   return simplify_transformation (array, dim, mask, 0, do_bit_xor);
2940 }
2941
2942
2943
2944 gfc_expr *
2945 gfc_simplify_is_iostat_end (gfc_expr *x)
2946 {
2947   if (x->expr_type != EXPR_CONSTANT)
2948     return NULL;
2949
2950   return gfc_get_logical_expr (gfc_default_logical_kind, &x->where,
2951                                mpz_cmp_si (x->value.integer,
2952                                            LIBERROR_END) == 0);
2953 }
2954
2955
2956 gfc_expr *
2957 gfc_simplify_is_iostat_eor (gfc_expr *x)
2958 {
2959   if (x->expr_type != EXPR_CONSTANT)
2960     return NULL;
2961
2962   return gfc_get_logical_expr (gfc_default_logical_kind, &x->where,
2963                                mpz_cmp_si (x->value.integer,
2964                                            LIBERROR_EOR) == 0);
2965 }
2966
2967
2968 gfc_expr *
2969 gfc_simplify_isnan (gfc_expr *x)
2970 {
2971   if (x->expr_type != EXPR_CONSTANT)
2972     return NULL;
2973
2974   return gfc_get_logical_expr (gfc_default_logical_kind, &x->where,
2975                                mpfr_nan_p (x->value.real));
2976 }
2977
2978
2979 /* Performs a shift on its first argument.  Depending on the last
2980    argument, the shift can be arithmetic, i.e. with filling from the
2981    left like in the SHIFTA intrinsic.  */
2982 static gfc_expr *
2983 simplify_shift (gfc_expr *e, gfc_expr *s, const char *name,
2984                 bool arithmetic, int direction)
2985 {
2986   gfc_expr *result;
2987   int ashift, *bits, i, k, bitsize, shift;
2988
2989   if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
2990     return NULL;
2991
2992   gfc_extract_int (s, &shift);
2993
2994   k = gfc_validate_kind (BT_INTEGER, e->ts.kind, false);
2995   bitsize = gfc_integer_kinds[k].bit_size;
2996
2997   result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);
2998
2999   if (shift == 0)
3000     {
3001       mpz_set (result->value.integer, e->value.integer);
3002       return result;
3003     }
3004
3005   if (direction > 0 && shift < 0)
3006     {
3007       /* Left shift, as in SHIFTL.  */
3008       gfc_error ("Second argument of %s is negative at %L", name, &e->where);
3009       return &gfc_bad_expr;
3010     }
3011   else if (direction < 0)
3012     {
3013       /* Right shift, as in SHIFTR or SHIFTA.  */
3014       if (shift < 0)
3015         {
3016           gfc_error ("Second argument of %s is negative at %L",
3017                      name, &e->where);
3018           return &gfc_bad_expr;
3019         }
3020
3021       shift = -shift;
3022     }
3023
3024   ashift = (shift >= 0 ? shift : -shift);
3025
3026   if (ashift > bitsize)
3027     {
3028       gfc_error ("Magnitude of second argument of %s exceeds bit size "
3029                  "at %L", name, &e->where);
3030       return &gfc_bad_expr;
3031     }
3032
3033   bits = XCNEWVEC (int, bitsize);
3034
3035   for (i = 0; i < bitsize; i++)
3036     bits[i] = mpz_tstbit (e->value.integer, i);
3037
3038   if (shift > 0)
3039     {
3040       /* Left shift.  */
3041       for (i = 0; i < shift; i++)
3042         mpz_clrbit (result->value.integer, i);
3043
3044       for (i = 0; i < bitsize - shift; i++)
3045         {
3046           if (bits[i] == 0)
3047             mpz_clrbit (result->value.integer, i + shift);
3048           else
3049             mpz_setbit (result->value.integer, i + shift);
3050         }
3051     }
3052   else
3053     {
3054       /* Right shift.  */
3055       if (arithmetic && bits[bitsize - 1])
3056         for (i = bitsize - 1; i >= bitsize - ashift; i--)
3057           mpz_setbit (result->value.integer, i);
3058       else
3059         for (i = bitsize - 1; i >= bitsize - ashift; i--)
3060           mpz_clrbit (result->value.integer, i);
3061
3062       for (i = bitsize - 1; i >= ashift; i--)
3063         {
3064           if (bits[i] == 0)
3065             mpz_clrbit (result->value.integer, i - ashift);
3066           else
3067             mpz_setbit (result->value.integer, i - ashift);
3068         }
3069     }
3070
3071   convert_mpz_to_signed (result->value.integer, bitsize);
3072   free (bits);
3073
3074   return result;
3075 }
3076
3077
3078 gfc_expr *
3079 gfc_simplify_ishft (gfc_expr *e, gfc_expr *s)
3080 {
3081   return simplify_shift (e, s, "ISHFT", false, 0);
3082 }
3083
3084
3085 gfc_expr *
3086 gfc_simplify_lshift (gfc_expr *e, gfc_expr *s)
3087 {
3088   return simplify_shift (e, s, "LSHIFT", false, 1);
3089 }
3090
3091
3092 gfc_expr *
3093 gfc_simplify_rshift (gfc_expr *e, gfc_expr *s)
3094 {
3095   return simplify_shift (e, s, "RSHIFT", true, -1);
3096 }
3097
3098
3099 gfc_expr *
3100 gfc_simplify_shifta (gfc_expr *e, gfc_expr *s)
3101 {
3102   return simplify_shift (e, s, "SHIFTA", true, -1);
3103 }
3104
3105
3106 gfc_expr *
3107 gfc_simplify_shiftl (gfc_expr *e, gfc_expr *s)
3108 {
3109   return simplify_shift (e, s, "SHIFTL", false, 1);
3110 }
3111
3112
3113 gfc_expr *
3114 gfc_simplify_shiftr (gfc_expr *e, gfc_expr *s)
3115 {
3116   return simplify_shift (e, s, "SHIFTR", false, -1);
3117 }
3118
3119
3120 gfc_expr *
3121 gfc_simplify_ishftc (gfc_expr *e, gfc_expr *s, gfc_expr *sz)
3122 {
3123   gfc_expr *result;
3124   int shift, ashift, isize, ssize, delta, k;
3125   int i, *bits;
3126
3127   if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
3128     return NULL;
3129
3130   gfc_extract_int (s, &shift);
3131
3132   k = gfc_validate_kind (e->ts.type, e->ts.kind, false);
3133   isize = gfc_integer_kinds[k].bit_size;
3134
3135   if (sz != NULL)
3136     {
3137       if (sz->expr_type != EXPR_CONSTANT)
3138         return NULL;
3139
3140       gfc_extract_int (sz, &ssize);
3141
3142     }
3143   else
3144     ssize = isize;
3145
3146   if (shift >= 0)
3147     ashift = shift;
3148   else
3149     ashift = -shift;
3150
3151   if (ashift > ssize)
3152     {
3153       if (sz == NULL)
3154         gfc_error ("Magnitude of second argument of ISHFTC exceeds "
3155                    "BIT_SIZE of first argument at %L", &s->where);
3156       return &gfc_bad_expr;
3157     }
3158
3159   result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);
3160
3161   mpz_set (result->value.integer, e->value.integer);
3162
3163   if (shift == 0)
3164     return result;
3165
3166   convert_mpz_to_unsigned (result->value.integer, isize);
3167
3168   bits = XCNEWVEC (int, ssize);
3169
3170   for (i = 0; i < ssize; i++)
3171     bits[i] = mpz_tstbit (e->value.integer, i);
3172
3173   delta = ssize - ashift;
3174
3175   if (shift > 0)
3176     {
3177       for (i = 0; i < delta; i++)
3178         {
3179           if (bits[i] == 0)
3180             mpz_clrbit (result->value.integer, i + shift);
3181           else
3182             mpz_setbit (result->value.integer, i + shift);
3183         }
3184
3185       for (i = delta; i < ssize; i++)
3186         {
3187           if (bits[i] == 0)
3188             mpz_clrbit (result->value.integer, i - delta);
3189           else
3190             mpz_setbit (result->value.integer, i - delta);
3191         }
3192     }
3193   else
3194     {
3195       for (i = 0; i < ashift; i++)
3196         {
3197           if (bits[i] == 0)
3198             mpz_clrbit (result->value.integer, i + delta);
3199           else
3200             mpz_setbit (result->value.integer, i + delta);
3201         }
3202
3203       for (i = ashift; i < ssize; i++)
3204         {
3205           if (bits[i] == 0)
3206             mpz_clrbit (result->value.integer, i + shift);
3207           else
3208             mpz_setbit (result->value.integer, i + shift);
3209         }
3210     }
3211
3212   convert_mpz_to_signed (result->value.integer, isize);
3213
3214   free (bits);
3215   return result;
3216 }
3217
3218
3219 gfc_expr *
3220 gfc_simplify_kind (gfc_expr *e)
3221 {
3222   return gfc_get_int_expr (gfc_default_integer_kind, NULL, e->ts.kind);
3223 }
3224
3225
3226 static gfc_expr *
3227 simplify_bound_dim (gfc_expr *array, gfc_expr *kind, int d, int upper,
3228                     gfc_array_spec *as, gfc_ref *ref, bool coarray)
3229 {
3230   gfc_expr *l, *u, *result;
3231   int k;
3232
3233   k = get_kind (BT_INTEGER, kind, upper ? "UBOUND" : "LBOUND",
3234                 gfc_default_integer_kind); 
3235   if (k == -1)
3236     return &gfc_bad_expr;
3237
3238   result = gfc_get_constant_expr (BT_INTEGER, k, &array->where);
3239
3240   /* For non-variables, LBOUND(expr, DIM=n) = 1 and
3241      UBOUND(expr, DIM=n) = SIZE(expr, DIM=n).  */
3242   if (!coarray && array->expr_type != EXPR_VARIABLE)
3243     {
3244       if (upper)
3245         {
3246           gfc_expr* dim = result;
3247           mpz_set_si (dim->value.integer, d);
3248
3249           result = simplify_size (array, dim, k);
3250           gfc_free_expr (dim);
3251           if (!result)
3252             goto returnNull;
3253         }
3254       else
3255         mpz_set_si (result->value.integer, 1);
3256
3257       goto done;
3258     }
3259
3260   /* Otherwise, we have a variable expression.  */
3261   gcc_assert (array->expr_type == EXPR_VARIABLE);
3262   gcc_assert (as);
3263
3264   if (gfc_resolve_array_spec (as, 0) == FAILURE)
3265     return NULL;
3266
3267   /* The last dimension of an assumed-size array is special.  */
3268   if ((!coarray && d == as->rank && as->type == AS_ASSUMED_SIZE && !upper)
3269       || (coarray && d == as->rank + as->corank
3270           && (!upper || gfc_option.coarray == GFC_FCOARRAY_SINGLE)))
3271     {
3272       if (as->lower[d-1]->expr_type == EXPR_CONSTANT)
3273         {
3274           gfc_free_expr (result);
3275           return gfc_copy_expr (as->lower[d-1]);
3276         }
3277
3278       goto returnNull;
3279     }
3280
3281   result = gfc_get_constant_expr (BT_INTEGER, k, &array->where);
3282
3283   /* Then, we need to know the extent of the given dimension.  */
3284   if (coarray || ref->u.ar.type == AR_FULL)
3285     {
3286       l = as->lower[d-1];
3287       u = as->upper[d-1];
3288
3289       if (l->expr_type != EXPR_CONSTANT || u == NULL
3290           || u->expr_type != EXPR_CONSTANT)
3291         goto returnNull;
3292
3293       if (mpz_cmp (l->value.integer, u->value.integer) > 0)
3294         {
3295           /* Zero extent.  */
3296           if (upper)
3297             mpz_set_si (result->value.integer, 0);
3298           else
3299             mpz_set_si (result->value.integer, 1);
3300         }
3301       else
3302         {
3303           /* Nonzero extent.  */
3304           if (upper)
3305             mpz_set (result->value.integer, u->value.integer);
3306           else
3307             mpz_set (result->value.integer, l->value.integer);
3308         }
3309     }
3310   else
3311     {
3312       if (upper)
3313         {
3314           if (gfc_ref_dimen_size (&ref->u.ar, d-1, &result->value.integer, NULL)
3315               != SUCCESS)
3316             goto returnNull;
3317         }
3318       else
3319         mpz_set_si (result->value.integer, (long int) 1);
3320     }
3321
3322 done:
3323   return range_check (result, upper ? "UBOUND" : "LBOUND");
3324
3325 returnNull:
3326   gfc_free_expr (result);
3327   return NULL;
3328 }
3329
3330
3331 static gfc_expr *
3332 simplify_bound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
3333 {
3334   gfc_ref *ref;
3335   gfc_array_spec *as;
3336   int d;
3337
3338   if (array->ts.type == BT_CLASS)
3339     return NULL;
3340
3341   if (array->expr_type != EXPR_VARIABLE)
3342     {
3343       as = NULL;
3344       ref = NULL;
3345       goto done;
3346     }
3347
3348   /* Follow any component references.  */
3349   as = array->symtree->n.sym->as;
3350   for (ref = array->ref; ref; ref = ref->next)
3351     {
3352       switch (ref->type)
3353         {
3354         case REF_ARRAY:
3355           switch (ref->u.ar.type)
3356             {
3357             case AR_ELEMENT:
3358               as = NULL;
3359               continue;
3360
3361             case AR_FULL:
3362               /* We're done because 'as' has already been set in the
3363                  previous iteration.  */
3364               if (!ref->next)
3365                 goto done;
3366
3367             /* Fall through.  */
3368
3369             case AR_UNKNOWN:
3370               return NULL;
3371
3372             case AR_SECTION:
3373               as = ref->u.ar.as;
3374               goto done;
3375             }
3376
3377           gcc_unreachable ();
3378
3379         case REF_COMPONENT:
3380           as = ref->u.c.component->as;
3381           continue;
3382
3383         case REF_SUBSTRING:
3384           continue;
3385         }
3386     }
3387
3388   gcc_unreachable ();
3389
3390  done:
3391
3392   if (as && (as->type == AS_DEFERRED || as->type == AS_ASSUMED_SHAPE))
3393     return NULL;
3394
3395   if (dim == NULL)
3396     {
3397       /* Multi-dimensional bounds.  */
3398       gfc_expr *bounds[GFC_MAX_DIMENSIONS];
3399       gfc_expr *e;
3400       int k;
3401
3402       /* UBOUND(ARRAY) is not valid for an assumed-size array.  */
3403       if (upper && as && as->type == AS_ASSUMED_SIZE)
3404         {
3405           /* An error message will be emitted in
3406              check_assumed_size_reference (resolve.c).  */
3407           return &gfc_bad_expr;
3408         }
3409
3410       /* Simplify the bounds for each dimension.  */
3411       for (d = 0; d < array->rank; d++)
3412         {
3413           bounds[d] = simplify_bound_dim (array, kind, d + 1, upper, as, ref,
3414                                           false);
3415           if (bounds[d] == NULL || bounds[d] == &gfc_bad_expr)
3416             {
3417               int j;
3418
3419               for (j = 0; j < d; j++)
3420                 gfc_free_expr (bounds[j]);
3421               return bounds[d];
3422             }
3423         }
3424
3425       /* Allocate the result expression.  */
3426       k = get_kind (BT_INTEGER, kind, upper ? "UBOUND" : "LBOUND",
3427                     gfc_default_integer_kind);
3428       if (k == -1)
3429         return &gfc_bad_expr;
3430
3431       e = gfc_get_array_expr (BT_INTEGER, k, &array->where);
3432
3433       /* The result is a rank 1 array; its size is the rank of the first
3434          argument to {L,U}BOUND.  */
3435       e->rank = 1;
3436       e->shape = gfc_get_shape (1);
3437       mpz_init_set_ui (e->shape[0], array->rank);
3438
3439       /* Create the constructor for this array.  */
3440       for (d = 0; d < array->rank; d++)
3441         gfc_constructor_append_expr (&e->value.constructor,
3442                                      bounds[d], &e->where);
3443
3444       return e;
3445     }
3446   else
3447     {
3448       /* A DIM argument is specified.  */
3449       if (dim->expr_type != EXPR_CONSTANT)
3450         return NULL;
3451
3452       d = mpz_get_si (dim->value.integer);
3453
3454       if (d < 1 || d > array->rank
3455           || (d == array->rank && as && as->type == AS_ASSUMED_SIZE && upper))
3456         {
3457           gfc_error ("DIM argument at %L is out of bounds", &dim->where);
3458           return &gfc_bad_expr;
3459         }
3460
3461       return simplify_bound_dim (array, kind, d, upper, as, ref, false);
3462     }
3463 }
3464
3465
3466 static gfc_expr *
3467 simplify_cobound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
3468 {
3469   gfc_ref *ref;
3470   gfc_array_spec *as;
3471   int d;
3472
3473   if (array->expr_type != EXPR_VARIABLE)
3474     return NULL;
3475
3476   /* Follow any component references.  */
3477   as = (array->ts.type == BT_CLASS && array->ts.u.derived->components)
3478        ? array->ts.u.derived->components->as
3479        : array->symtree->n.sym->as;
3480   for (ref = array->ref; ref; ref = ref->next)
3481     {
3482       switch (ref->type)
3483         {
3484         case REF_ARRAY:
3485           switch (ref->u.ar.type)
3486             {
3487             case AR_ELEMENT:
3488               if (ref->u.ar.as->corank > 0)
3489                 {
3490                   gcc_assert (as == ref->u.ar.as);
3491                   goto done;
3492                 }
3493               as = NULL;
3494               continue;
3495
3496             case AR_FULL:
3497               /* We're done because 'as' has already been set in the
3498                  previous iteration.  */
3499               if (!ref->next)
3500                 goto done;
3501
3502             /* Fall through.  */
3503
3504             case AR_UNKNOWN:
3505               return NULL;
3506
3507             case AR_SECTION:
3508               as = ref->u.ar.as;
3509               goto done;
3510             }
3511
3512           gcc_unreachable ();
3513
3514         case REF_COMPONENT:
3515           as = ref->u.c.component->as;
3516           continue;
3517
3518         case REF_SUBSTRING:
3519           continue;
3520         }
3521     }
3522
3523   if (!as)
3524     gcc_unreachable ();
3525
3526  done:
3527
3528   if (as->cotype == AS_DEFERRED || as->cotype == AS_ASSUMED_SHAPE)
3529     return NULL;
3530
3531   if (dim == NULL)
3532     {
3533       /* Multi-dimensional cobounds.  */
3534       gfc_expr *bounds[GFC_MAX_DIMENSIONS];
3535       gfc_expr *e;
3536       int k;
3537
3538       /* Simplify the cobounds for each dimension.  */
3539       for (d = 0; d < as->corank; d++)
3540         {
3541           bounds[d] = simplify_bound_dim (array, kind, d + 1 + as->rank,
3542                                           upper, as, ref, true);
3543           if (bounds[d] == NULL || bounds[d] == &gfc_bad_expr)
3544             {
3545               int j;
3546
3547               for (j = 0; j < d; j++)
3548                 gfc_free_expr (bounds[j]);
3549               return bounds[d];
3550             }
3551         }
3552
3553       /* Allocate the result expression.  */
3554       e = gfc_get_expr ();
3555       e->where = array->where;
3556       e->expr_type = EXPR_ARRAY;
3557       e->ts.type = BT_INTEGER;
3558       k = get_kind (BT_INTEGER, kind, upper ? "UCOBOUND" : "LCOBOUND",
3559                     gfc_default_integer_kind); 
3560       if (k == -1)
3561         {
3562           gfc_free_expr (e);
3563           return &gfc_bad_expr;
3564         }
3565       e->ts.kind = k;
3566
3567       /* The result is a rank 1 array; its size is the rank of the first
3568          argument to {L,U}COBOUND.  */
3569       e->rank = 1;
3570       e->shape = gfc_get_shape (1);
3571       mpz_init_set_ui (e->shape[0], as->corank);
3572
3573       /* Create the constructor for this array.  */
3574       for (d = 0; d < as->corank; d++)
3575         gfc_constructor_append_expr (&e->value.constructor,
3576                                      bounds[d], &e->where);
3577       return e;
3578     }
3579   else
3580     {
3581       /* A DIM argument is specified.  */
3582       if (dim->expr_type != EXPR_CONSTANT)
3583         return NULL;
3584
3585       d = mpz_get_si (dim->value.integer);
3586
3587       if (d < 1 || d > as->corank)
3588         {
3589           gfc_error ("DIM argument at %L is out of bounds", &dim->where);
3590           return &gfc_bad_expr;
3591         }
3592
3593       return simplify_bound_dim (array, kind, d+as->rank, upper, as, ref, true);
3594     }
3595 }
3596
3597
3598 gfc_expr *
3599 gfc_simplify_lbound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind)
3600 {
3601   return simplify_bound (array, dim, kind, 0);
3602 }
3603
3604
3605 gfc_expr *
3606 gfc_simplify_lcobound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind)
3607 {
3608   return simplify_cobound (array, dim, kind, 0);
3609 }
3610
3611 gfc_expr *
3612 gfc_simplify_leadz (gfc_expr *e)
3613 {
3614   unsigned long lz, bs;
3615   int i;
3616
3617   if (e->expr_type != EXPR_CONSTANT)
3618     return NULL;
3619
3620   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
3621   bs = gfc_integer_kinds[i].bit_size;
3622   if (mpz_cmp_si (e->value.integer, 0) == 0)
3623     lz = bs;
3624   else if (mpz_cmp_si (e->value.integer, 0) < 0)
3625     lz = 0;
3626   else
3627     lz = bs - mpz_sizeinbase (e->value.integer, 2);
3628
3629   return gfc_get_int_expr (gfc_default_integer_kind, &e->where, lz);
3630 }
3631
3632
3633 gfc_expr *
3634 gfc_simplify_len (gfc_expr *e, gfc_expr *kind)
3635 {
3636   gfc_expr *result;
3637   int k = get_kind (BT_INTEGER, kind, "LEN", gfc_default_integer_kind);
3638
3639   if (k == -1)
3640     return &gfc_bad_expr;
3641
3642   if (e->expr_type == EXPR_CONSTANT)
3643     {
3644       result = gfc_get_constant_expr (BT_INTEGER, k, &e->where);
3645       mpz_set_si (result->value.integer, e->value.character.length);
3646       return range_check (result, "LEN");
3647     }
3648   else if (e->ts.u.cl != NULL && e->ts.u.cl->length != NULL
3649            && e->ts.u.cl->length->expr_type == EXPR_CONSTANT
3650            && e->ts.u.cl->length->ts.type == BT_INTEGER)
3651     {
3652       result = gfc_get_constant_expr (BT_INTEGER, k, &e->where);
3653       mpz_set (result->value.integer, e->ts.u.cl->length->value.integer);
3654       return range_check (result, "LEN");
3655     }
3656   else
3657     return NULL;
3658 }
3659
3660
3661 gfc_expr *
3662 gfc_simplify_len_trim (gfc_expr *e, gfc_expr *kind)
3663 {
3664   gfc_expr *result;
3665   int count, len, i;
3666   int k = get_kind (BT_INTEGER, kind, "LEN_TRIM", gfc_default_integer_kind);
3667
3668   if (k == -1)
3669     return &gfc_bad_expr;
3670
3671   if (e->expr_type != EXPR_CONSTANT)
3672     return NULL;
3673
3674   len = e->value.character.length;
3675   for (count = 0, i = 1; i <= len; i++)
3676     if (e->value.character.string[len - i] == ' ')
3677       count++;
3678     else
3679       break;
3680
3681   result = gfc_get_int_expr (k, &e->where, len - count);
3682   return range_check (result, "LEN_TRIM");
3683 }
3684
3685 gfc_expr *
3686 gfc_simplify_lgamma (gfc_expr *x)
3687 {
3688   gfc_expr *result;
3689   int sg;
3690
3691   if (x->expr_type != EXPR_CONSTANT)
3692     return NULL;
3693
3694   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
3695   mpfr_lgamma (result->value.real, &sg, x->value.real, GFC_RND_MODE);
3696
3697   return range_check (result, "LGAMMA");
3698 }
3699
3700