OSDN Git Service

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