1 /* Simplify intrinsic functions at compile-time.
2 Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009,
3 2010, 2011 Free Software Foundation, Inc.
4 Contributed by Andy Vaught & Katherine Holcomb
6 This file is part of GCC.
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
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
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/>. */
27 #include "intrinsic.h"
28 #include "target-memory.h"
29 #include "constructor.h"
30 #include "version.h" /* For version_string. */
33 gfc_expr gfc_bad_expr;
36 /* Note that 'simplification' is not just transforming expressions.
37 For functions that are not simplified at compile time, range
38 checking is done if possible.
40 The return convention is that each simplification function returns:
42 A new expression node corresponding to the simplified arguments.
43 The original arguments are destroyed by the caller, and must not
44 be a part of the new expression.
46 NULL pointer indicating that no simplification was possible and
47 the original expression should remain intact.
49 An expression pointer to gfc_bad_expr (a static placeholder)
50 indicating that some error has prevented simplification. The
51 error is generated within the function and should be propagated
54 By the time a simplification function gets control, it has been
55 decided that the function call is really supposed to be the
56 intrinsic. No type checking is strictly necessary, since only
57 valid types will be passed on. On the other hand, a simplification
58 subroutine may have to look at the type of an argument as part of
61 Array arguments are only passed to these subroutines that implement
62 the simplification of transformational intrinsics.
64 The functions in this file don't have much comment with them, but
65 everything is reasonably straight-forward. The Standard, chapter 13
66 is the best comment you'll find for this file anyway. */
68 /* Range checks an expression node. If all goes well, returns the
69 node, otherwise returns &gfc_bad_expr and frees the node. */
72 range_check (gfc_expr *result, const char *name)
77 if (result->expr_type != EXPR_CONSTANT)
80 switch (gfc_range_check (result))
86 gfc_error ("Result of %s overflows its kind at %L", name,
91 gfc_error ("Result of %s underflows its kind at %L", name,
96 gfc_error ("Result of %s is NaN at %L", name, &result->where);
100 gfc_error ("Result of %s gives range error for its kind at %L", name,
105 gfc_free_expr (result);
106 return &gfc_bad_expr;
110 /* A helper function that gets an optional and possibly missing
111 kind parameter. Returns the kind, -1 if something went wrong. */
114 get_kind (bt type, gfc_expr *k, const char *name, int default_kind)
121 if (k->expr_type != EXPR_CONSTANT)
123 gfc_error ("KIND parameter of %s at %L must be an initialization "
124 "expression", name, &k->where);
128 if (gfc_extract_int (k, &kind) != NULL
129 || gfc_validate_kind (type, kind, true) < 0)
131 gfc_error ("Invalid KIND parameter of %s at %L", name, &k->where);
139 /* Converts an mpz_t signed variable into an unsigned one, assuming
140 two's complement representations and a binary width of bitsize.
141 The conversion is a no-op unless x is negative; otherwise, it can
142 be accomplished by masking out the high bits. */
145 convert_mpz_to_unsigned (mpz_t x, int bitsize)
151 /* Confirm that no bits above the signed range are unset. */
152 gcc_assert (mpz_scan0 (x, bitsize-1) == ULONG_MAX);
154 mpz_init_set_ui (mask, 1);
155 mpz_mul_2exp (mask, mask, bitsize);
156 mpz_sub_ui (mask, mask, 1);
158 mpz_and (x, x, mask);
164 /* Confirm that no bits above the signed range are set. */
165 gcc_assert (mpz_scan1 (x, bitsize-1) == ULONG_MAX);
170 /* Converts an mpz_t unsigned variable into a signed one, assuming
171 two's complement representations and a binary width of bitsize.
172 If the bitsize-1 bit is set, this is taken as a sign bit and
173 the number is converted to the corresponding negative number. */
176 convert_mpz_to_signed (mpz_t x, int bitsize)
180 /* Confirm that no bits above the unsigned range are set. */
181 gcc_assert (mpz_scan1 (x, bitsize) == ULONG_MAX);
183 if (mpz_tstbit (x, bitsize - 1) == 1)
185 mpz_init_set_ui (mask, 1);
186 mpz_mul_2exp (mask, mask, bitsize);
187 mpz_sub_ui (mask, mask, 1);
189 /* We negate the number by hand, zeroing the high bits, that is
190 make it the corresponding positive number, and then have it
191 negated by GMP, giving the correct representation of the
194 mpz_add_ui (x, x, 1);
195 mpz_and (x, x, mask);
204 /* In-place convert BOZ to REAL of the specified kind. */
207 convert_boz (gfc_expr *x, int kind)
209 if (x && x->ts.type == BT_INTEGER && x->is_boz)
216 if (!gfc_convert_boz (x, &ts))
217 return &gfc_bad_expr;
224 /* Test that the expression is an constant array. */
227 is_constant_array_expr (gfc_expr *e)
234 if (e->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (e))
237 for (c = gfc_constructor_first (e->value.constructor);
238 c; c = gfc_constructor_next (c))
239 if (c->expr->expr_type != EXPR_CONSTANT
240 && c->expr->expr_type != EXPR_STRUCTURE)
247 /* Initialize a transformational result expression with a given value. */
250 init_result_expr (gfc_expr *e, int init, gfc_expr *array)
252 if (e && e->expr_type == EXPR_ARRAY)
254 gfc_constructor *ctor = gfc_constructor_first (e->value.constructor);
257 init_result_expr (ctor->expr, init, array);
258 ctor = gfc_constructor_next (ctor);
261 else if (e && e->expr_type == EXPR_CONSTANT)
263 int i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
270 e->value.logical = (init ? 1 : 0);
275 mpz_set (e->value.integer, gfc_integer_kinds[i].min_int);
276 else if (init == INT_MAX)
277 mpz_set (e->value.integer, gfc_integer_kinds[i].huge);
279 mpz_set_si (e->value.integer, init);
285 mpfr_set (e->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
286 mpfr_neg (e->value.real, e->value.real, GFC_RND_MODE);
288 else if (init == INT_MAX)
289 mpfr_set (e->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
291 mpfr_set_si (e->value.real, init, GFC_RND_MODE);
295 mpc_set_si (e->value.complex, init, GFC_MPC_RND_MODE);
301 gfc_expr *len = gfc_simplify_len (array, NULL);
302 gfc_extract_int (len, &length);
303 string = gfc_get_wide_string (length + 1);
304 gfc_wide_memset (string, 0, length);
306 else if (init == INT_MAX)
308 gfc_expr *len = gfc_simplify_len (array, NULL);
309 gfc_extract_int (len, &length);
310 string = gfc_get_wide_string (length + 1);
311 gfc_wide_memset (string, 255, length);
316 string = gfc_get_wide_string (1);
319 string[length] = '\0';
320 e->value.character.length = length;
321 e->value.character.string = string;
333 /* Helper function for gfc_simplify_dot_product() and gfc_simplify_matmul. */
336 compute_dot_product (gfc_expr *matrix_a, int stride_a, int offset_a,
337 gfc_expr *matrix_b, int stride_b, int offset_b)
339 gfc_expr *result, *a, *b;
341 result = gfc_get_constant_expr (matrix_a->ts.type, matrix_a->ts.kind,
343 init_result_expr (result, 0, NULL);
345 a = gfc_constructor_lookup_expr (matrix_a->value.constructor, offset_a);
346 b = gfc_constructor_lookup_expr (matrix_b->value.constructor, offset_b);
349 /* Copying of expressions is required as operands are free'd
350 by the gfc_arith routines. */
351 switch (result->ts.type)
354 result = gfc_or (result,
355 gfc_and (gfc_copy_expr (a),
362 result = gfc_add (result,
363 gfc_multiply (gfc_copy_expr (a),
371 offset_a += stride_a;
372 a = gfc_constructor_lookup_expr (matrix_a->value.constructor, offset_a);
374 offset_b += stride_b;
375 b = gfc_constructor_lookup_expr (matrix_b->value.constructor, offset_b);
382 /* Build a result expression for transformational intrinsics,
386 transformational_result (gfc_expr *array, gfc_expr *dim, bt type,
387 int kind, locus* where)
392 if (!dim || array->rank == 1)
393 return gfc_get_constant_expr (type, kind, where);
395 result = gfc_get_array_expr (type, kind, where);
396 result->shape = gfc_copy_shape_excluding (array->shape, array->rank, dim);
397 result->rank = array->rank - 1;
399 /* gfc_array_size() would count the number of elements in the constructor,
400 we have not built those yet. */
402 for (i = 0; i < result->rank; ++i)
403 nelem *= mpz_get_ui (result->shape[i]);
405 for (i = 0; i < nelem; ++i)
407 gfc_constructor_append_expr (&result->value.constructor,
408 gfc_get_constant_expr (type, kind, where),
416 typedef gfc_expr* (*transformational_op)(gfc_expr*, gfc_expr*);
418 /* Wrapper function, implements 'op1 += 1'. Only called if MASK
419 of COUNT intrinsic is .TRUE..
421 Interface and implimentation mimics arith functions as
422 gfc_add, gfc_multiply, etc. */
424 static gfc_expr* gfc_count (gfc_expr *op1, gfc_expr *op2)
428 gcc_assert (op1->ts.type == BT_INTEGER);
429 gcc_assert (op2->ts.type == BT_LOGICAL);
430 gcc_assert (op2->value.logical);
432 result = gfc_copy_expr (op1);
433 mpz_add_ui (result->value.integer, result->value.integer, 1);
441 /* Transforms an ARRAY with operation OP, according to MASK, to a
442 scalar RESULT. E.g. called if
444 REAL, PARAMETER :: array(n, m) = ...
445 REAL, PARAMETER :: s = SUM(array)
447 where OP == gfc_add(). */
450 simplify_transformation_to_scalar (gfc_expr *result, gfc_expr *array, gfc_expr *mask,
451 transformational_op op)
454 gfc_constructor *array_ctor, *mask_ctor;
456 /* Shortcut for constant .FALSE. MASK. */
458 && mask->expr_type == EXPR_CONSTANT
459 && !mask->value.logical)
462 array_ctor = gfc_constructor_first (array->value.constructor);
464 if (mask && mask->expr_type == EXPR_ARRAY)
465 mask_ctor = gfc_constructor_first (mask->value.constructor);
469 a = array_ctor->expr;
470 array_ctor = gfc_constructor_next (array_ctor);
472 /* A constant MASK equals .TRUE. here and can be ignored. */
476 mask_ctor = gfc_constructor_next (mask_ctor);
477 if (!m->value.logical)
481 result = op (result, gfc_copy_expr (a));
487 /* Transforms an ARRAY with operation OP, according to MASK, to an
488 array RESULT. E.g. called if
490 REAL, PARAMETER :: array(n, m) = ...
491 REAL, PARAMETER :: s(n) = PROD(array, DIM=1)
493 where OP == gfc_multiply(). The result might be post processed using post_op. */
496 simplify_transformation_to_array (gfc_expr *result, gfc_expr *array, gfc_expr *dim,
497 gfc_expr *mask, transformational_op op,
498 transformational_op post_op)
501 int done, i, n, arraysize, resultsize, dim_index, dim_extent, dim_stride;
502 gfc_expr **arrayvec, **resultvec, **base, **src, **dest;
503 gfc_constructor *array_ctor, *mask_ctor, *result_ctor;
505 int count[GFC_MAX_DIMENSIONS], extent[GFC_MAX_DIMENSIONS],
506 sstride[GFC_MAX_DIMENSIONS], dstride[GFC_MAX_DIMENSIONS],
507 tmpstride[GFC_MAX_DIMENSIONS];
509 /* Shortcut for constant .FALSE. MASK. */
511 && mask->expr_type == EXPR_CONSTANT
512 && !mask->value.logical)
515 /* Build an indexed table for array element expressions to minimize
516 linked-list traversal. Masked elements are set to NULL. */
517 gfc_array_size (array, &size);
518 arraysize = mpz_get_ui (size);
521 arrayvec = XCNEWVEC (gfc_expr*, arraysize);
523 array_ctor = gfc_constructor_first (array->value.constructor);
525 if (mask && mask->expr_type == EXPR_ARRAY)
526 mask_ctor = gfc_constructor_first (mask->value.constructor);
528 for (i = 0; i < arraysize; ++i)
530 arrayvec[i] = array_ctor->expr;
531 array_ctor = gfc_constructor_next (array_ctor);
535 if (!mask_ctor->expr->value.logical)
538 mask_ctor = gfc_constructor_next (mask_ctor);
542 /* Same for the result expression. */
543 gfc_array_size (result, &size);
544 resultsize = mpz_get_ui (size);
547 resultvec = XCNEWVEC (gfc_expr*, resultsize);
548 result_ctor = gfc_constructor_first (result->value.constructor);
549 for (i = 0; i < resultsize; ++i)
551 resultvec[i] = result_ctor->expr;
552 result_ctor = gfc_constructor_next (result_ctor);
555 gfc_extract_int (dim, &dim_index);
556 dim_index -= 1; /* zero-base index */
560 for (i = 0, n = 0; i < array->rank; ++i)
563 tmpstride[i] = (i == 0) ? 1 : tmpstride[i-1] * mpz_get_si (array->shape[i-1]);
566 dim_extent = mpz_get_si (array->shape[i]);
567 dim_stride = tmpstride[i];
571 extent[n] = mpz_get_si (array->shape[i]);
572 sstride[n] = tmpstride[i];
573 dstride[n] = (n == 0) ? 1 : dstride[n-1] * extent[n-1];
582 for (src = base, n = 0; n < dim_extent; src += dim_stride, ++n)
584 *dest = op (*dest, gfc_copy_expr (*src));
591 while (!done && count[n] == extent[n])
594 base -= sstride[n] * extent[n];
595 dest -= dstride[n] * extent[n];
598 if (n < result->rank)
609 /* Place updated expression in result constructor. */
610 result_ctor = gfc_constructor_first (result->value.constructor);
611 for (i = 0; i < resultsize; ++i)
614 result_ctor->expr = post_op (result_ctor->expr, resultvec[i]);
616 result_ctor->expr = resultvec[i];
617 result_ctor = gfc_constructor_next (result_ctor);
627 simplify_transformation (gfc_expr *array, gfc_expr *dim, gfc_expr *mask,
628 int init_val, transformational_op op)
632 if (!is_constant_array_expr (array)
633 || !gfc_is_constant_expr (dim))
637 && !is_constant_array_expr (mask)
638 && mask->expr_type != EXPR_CONSTANT)
641 result = transformational_result (array, dim, array->ts.type,
642 array->ts.kind, &array->where);
643 init_result_expr (result, init_val, NULL);
645 return !dim || array->rank == 1 ?
646 simplify_transformation_to_scalar (result, array, mask, op) :
647 simplify_transformation_to_array (result, array, dim, mask, op, NULL);
651 /********************** Simplification functions *****************************/
654 gfc_simplify_abs (gfc_expr *e)
658 if (e->expr_type != EXPR_CONSTANT)
664 result = gfc_get_constant_expr (BT_INTEGER, e->ts.kind, &e->where);
665 mpz_abs (result->value.integer, e->value.integer);
666 return range_check (result, "IABS");
669 result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
670 mpfr_abs (result->value.real, e->value.real, GFC_RND_MODE);
671 return range_check (result, "ABS");
674 gfc_set_model_kind (e->ts.kind);
675 result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
676 mpc_abs (result->value.real, e->value.complex, GFC_RND_MODE);
677 return range_check (result, "CABS");
680 gfc_internal_error ("gfc_simplify_abs(): Bad type");
686 simplify_achar_char (gfc_expr *e, gfc_expr *k, const char *name, bool ascii)
690 bool too_large = false;
692 if (e->expr_type != EXPR_CONSTANT)
695 kind = get_kind (BT_CHARACTER, k, name, gfc_default_character_kind);
697 return &gfc_bad_expr;
699 if (mpz_cmp_si (e->value.integer, 0) < 0)
701 gfc_error ("Argument of %s function at %L is negative", name,
703 return &gfc_bad_expr;
706 if (ascii && gfc_option.warn_surprising
707 && mpz_cmp_si (e->value.integer, 127) > 0)
708 gfc_warning ("Argument of %s function at %L outside of range [0,127]",
711 if (kind == 1 && mpz_cmp_si (e->value.integer, 255) > 0)
716 mpz_init_set_ui (t, 2);
717 mpz_pow_ui (t, t, 32);
718 mpz_sub_ui (t, t, 1);
719 if (mpz_cmp (e->value.integer, t) > 0)
726 gfc_error ("Argument of %s function at %L is too large for the "
727 "collating sequence of kind %d", name, &e->where, kind);
728 return &gfc_bad_expr;
731 result = gfc_get_character_expr (kind, &e->where, NULL, 1);
732 result->value.character.string[0] = mpz_get_ui (e->value.integer);
739 /* We use the processor's collating sequence, because all
740 systems that gfortran currently works on are ASCII. */
743 gfc_simplify_achar (gfc_expr *e, gfc_expr *k)
745 return simplify_achar_char (e, k, "ACHAR", true);
750 gfc_simplify_acos (gfc_expr *x)
754 if (x->expr_type != EXPR_CONSTANT)
760 if (mpfr_cmp_si (x->value.real, 1) > 0
761 || mpfr_cmp_si (x->value.real, -1) < 0)
763 gfc_error ("Argument of ACOS at %L must be between -1 and 1",
765 return &gfc_bad_expr;
767 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
768 mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
772 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
773 mpc_acos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
777 gfc_internal_error ("in gfc_simplify_acos(): Bad type");
780 return range_check (result, "ACOS");
784 gfc_simplify_acosh (gfc_expr *x)
788 if (x->expr_type != EXPR_CONSTANT)
794 if (mpfr_cmp_si (x->value.real, 1) < 0)
796 gfc_error ("Argument of ACOSH at %L must not be less than 1",
798 return &gfc_bad_expr;
801 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
802 mpfr_acosh (result->value.real, x->value.real, GFC_RND_MODE);
806 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
807 mpc_acosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
811 gfc_internal_error ("in gfc_simplify_acosh(): Bad type");
814 return range_check (result, "ACOSH");
818 gfc_simplify_adjustl (gfc_expr *e)
824 if (e->expr_type != EXPR_CONSTANT)
827 len = e->value.character.length;
829 for (count = 0, i = 0; i < len; ++i)
831 ch = e->value.character.string[i];
837 result = gfc_get_character_expr (e->ts.kind, &e->where, NULL, len);
838 for (i = 0; i < len - count; ++i)
839 result->value.character.string[i] = e->value.character.string[count + i];
846 gfc_simplify_adjustr (gfc_expr *e)
852 if (e->expr_type != EXPR_CONSTANT)
855 len = e->value.character.length;
857 for (count = 0, i = len - 1; i >= 0; --i)
859 ch = e->value.character.string[i];
865 result = gfc_get_character_expr (e->ts.kind, &e->where, NULL, len);
866 for (i = 0; i < count; ++i)
867 result->value.character.string[i] = ' ';
869 for (i = count; i < len; ++i)
870 result->value.character.string[i] = e->value.character.string[i - count];
877 gfc_simplify_aimag (gfc_expr *e)
881 if (e->expr_type != EXPR_CONSTANT)
884 result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
885 mpfr_set (result->value.real, mpc_imagref (e->value.complex), GFC_RND_MODE);
887 return range_check (result, "AIMAG");
892 gfc_simplify_aint (gfc_expr *e, gfc_expr *k)
894 gfc_expr *rtrunc, *result;
897 kind = get_kind (BT_REAL, k, "AINT", e->ts.kind);
899 return &gfc_bad_expr;
901 if (e->expr_type != EXPR_CONSTANT)
904 rtrunc = gfc_copy_expr (e);
905 mpfr_trunc (rtrunc->value.real, e->value.real);
907 result = gfc_real2real (rtrunc, kind);
909 gfc_free_expr (rtrunc);
911 return range_check (result, "AINT");
916 gfc_simplify_all (gfc_expr *mask, gfc_expr *dim)
918 return simplify_transformation (mask, dim, NULL, true, gfc_and);
923 gfc_simplify_dint (gfc_expr *e)
925 gfc_expr *rtrunc, *result;
927 if (e->expr_type != EXPR_CONSTANT)
930 rtrunc = gfc_copy_expr (e);
931 mpfr_trunc (rtrunc->value.real, e->value.real);
933 result = gfc_real2real (rtrunc, gfc_default_double_kind);
935 gfc_free_expr (rtrunc);
937 return range_check (result, "DINT");
942 gfc_simplify_dreal (gfc_expr *e)
944 gfc_expr *result = NULL;
946 if (e->expr_type != EXPR_CONSTANT)
949 result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
950 mpc_real (result->value.real, e->value.complex, GFC_RND_MODE);
952 return range_check (result, "DREAL");
957 gfc_simplify_anint (gfc_expr *e, gfc_expr *k)
962 kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
964 return &gfc_bad_expr;
966 if (e->expr_type != EXPR_CONSTANT)
969 result = gfc_get_constant_expr (e->ts.type, kind, &e->where);
970 mpfr_round (result->value.real, e->value.real);
972 return range_check (result, "ANINT");
977 gfc_simplify_and (gfc_expr *x, gfc_expr *y)
982 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
985 kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
990 result = gfc_get_constant_expr (BT_INTEGER, kind, &x->where);
991 mpz_and (result->value.integer, x->value.integer, y->value.integer);
992 return range_check (result, "AND");
995 return gfc_get_logical_expr (kind, &x->where,
996 x->value.logical && y->value.logical);
1005 gfc_simplify_any (gfc_expr *mask, gfc_expr *dim)
1007 return simplify_transformation (mask, dim, NULL, false, gfc_or);
1012 gfc_simplify_dnint (gfc_expr *e)
1016 if (e->expr_type != EXPR_CONSTANT)
1019 result = gfc_get_constant_expr (BT_REAL, gfc_default_double_kind, &e->where);
1020 mpfr_round (result->value.real, e->value.real);
1022 return range_check (result, "DNINT");
1027 gfc_simplify_asin (gfc_expr *x)
1031 if (x->expr_type != EXPR_CONSTANT)
1037 if (mpfr_cmp_si (x->value.real, 1) > 0
1038 || mpfr_cmp_si (x->value.real, -1) < 0)
1040 gfc_error ("Argument of ASIN at %L must be between -1 and 1",
1042 return &gfc_bad_expr;
1044 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1045 mpfr_asin (result->value.real, x->value.real, GFC_RND_MODE);
1049 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1050 mpc_asin (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1054 gfc_internal_error ("in gfc_simplify_asin(): Bad type");
1057 return range_check (result, "ASIN");
1062 gfc_simplify_asinh (gfc_expr *x)
1066 if (x->expr_type != EXPR_CONSTANT)
1069 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1074 mpfr_asinh (result->value.real, x->value.real, GFC_RND_MODE);
1078 mpc_asinh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1082 gfc_internal_error ("in gfc_simplify_asinh(): Bad type");
1085 return range_check (result, "ASINH");
1090 gfc_simplify_atan (gfc_expr *x)
1094 if (x->expr_type != EXPR_CONSTANT)
1097 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1102 mpfr_atan (result->value.real, x->value.real, GFC_RND_MODE);
1106 mpc_atan (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1110 gfc_internal_error ("in gfc_simplify_atan(): Bad type");
1113 return range_check (result, "ATAN");
1118 gfc_simplify_atanh (gfc_expr *x)
1122 if (x->expr_type != EXPR_CONSTANT)
1128 if (mpfr_cmp_si (x->value.real, 1) >= 0
1129 || mpfr_cmp_si (x->value.real, -1) <= 0)
1131 gfc_error ("Argument of ATANH at %L must be inside the range -1 "
1133 return &gfc_bad_expr;
1135 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1136 mpfr_atanh (result->value.real, x->value.real, GFC_RND_MODE);
1140 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1141 mpc_atanh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1145 gfc_internal_error ("in gfc_simplify_atanh(): Bad type");
1148 return range_check (result, "ATANH");
1153 gfc_simplify_atan2 (gfc_expr *y, gfc_expr *x)
1157 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1160 if (mpfr_sgn (y->value.real) == 0 && mpfr_sgn (x->value.real) == 0)
1162 gfc_error ("If first argument of ATAN2 %L is zero, then the "
1163 "second argument must not be zero", &x->where);
1164 return &gfc_bad_expr;
1167 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1168 mpfr_atan2 (result->value.real, y->value.real, x->value.real, GFC_RND_MODE);
1170 return range_check (result, "ATAN2");
1175 gfc_simplify_bessel_j0 (gfc_expr *x)
1179 if (x->expr_type != EXPR_CONSTANT)
1182 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1183 mpfr_j0 (result->value.real, x->value.real, GFC_RND_MODE);
1185 return range_check (result, "BESSEL_J0");
1190 gfc_simplify_bessel_j1 (gfc_expr *x)
1194 if (x->expr_type != EXPR_CONSTANT)
1197 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1198 mpfr_j1 (result->value.real, x->value.real, GFC_RND_MODE);
1200 return range_check (result, "BESSEL_J1");
1205 gfc_simplify_bessel_jn (gfc_expr *order, gfc_expr *x)
1210 if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
1213 n = mpz_get_si (order->value.integer);
1214 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1215 mpfr_jn (result->value.real, n, x->value.real, GFC_RND_MODE);
1217 return range_check (result, "BESSEL_JN");
1221 /* Simplify transformational form of JN and YN. */
1224 gfc_simplify_bessel_n2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x,
1231 mpfr_t x2rev, last1, last2;
1233 if (x->expr_type != EXPR_CONSTANT || order1->expr_type != EXPR_CONSTANT
1234 || order2->expr_type != EXPR_CONSTANT)
1237 n1 = mpz_get_si (order1->value.integer);
1238 n2 = mpz_get_si (order2->value.integer);
1239 result = gfc_get_array_expr (x->ts.type, x->ts.kind, &x->where);
1241 result->shape = gfc_get_shape (1);
1242 mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));
1247 /* Special case: x == 0; it is J0(0.0) == 1, JN(N > 0, 0.0) == 0; and
1248 YN(N, 0.0) = -Inf. */
1250 if (mpfr_cmp_ui (x->value.real, 0.0) == 0)
1252 if (!jn && gfc_option.flag_range_check)
1254 gfc_error ("Result of BESSEL_YN is -INF at %L", &result->where);
1255 gfc_free_expr (result);
1256 return &gfc_bad_expr;
1261 e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1262 mpfr_set_ui (e->value.real, 1, GFC_RND_MODE);
1263 gfc_constructor_append_expr (&result->value.constructor, e,
1268 for (i = n1; i <= n2; i++)
1270 e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1272 mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
1274 mpfr_set_inf (e->value.real, -1);
1275 gfc_constructor_append_expr (&result->value.constructor, e,
1282 /* Use the faster but more verbose recurrence algorithm. Bessel functions
1283 are stable for downward recursion and Neumann functions are stable
1284 for upward recursion. It is
1286 J(N-1, x) = x2rev * N * J(N, x) - J(N+1, x),
1287 Y(N+1, x) = x2rev * N * Y(N, x) - Y(N-1, x).
1288 Cf. http://dlmf.nist.gov/10.74#iv and http://dlmf.nist.gov/10.6#E1 */
1290 gfc_set_model_kind (x->ts.kind);
1292 /* Get first recursion anchor. */
1296 mpfr_jn (last1, n2, x->value.real, GFC_RND_MODE);
1298 mpfr_yn (last1, n1, x->value.real, GFC_RND_MODE);
1300 e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1301 mpfr_set (e->value.real, last1, GFC_RND_MODE);
1302 if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
1306 gfc_free_expr (result);
1307 return &gfc_bad_expr;
1309 gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
1317 /* Get second recursion anchor. */
1321 mpfr_jn (last2, n2-1, x->value.real, GFC_RND_MODE);
1323 mpfr_yn (last2, n1+1, x->value.real, GFC_RND_MODE);
1325 e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1326 mpfr_set (e->value.real, last2, GFC_RND_MODE);
1327 if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
1332 gfc_free_expr (result);
1333 return &gfc_bad_expr;
1336 gfc_constructor_insert_expr (&result->value.constructor, e, &x->where, -2);
1338 gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
1347 /* Start actual recursion. */
1350 mpfr_ui_div (x2rev, 2, x->value.real, GFC_RND_MODE);
1352 for (i = 2; i <= n2-n1; i++)
1354 e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1356 /* Special case: For YN, if the previous N gave -INF, set
1357 also N+1 to -INF. */
1358 if (!jn && !gfc_option.flag_range_check && mpfr_inf_p (last2))
1360 mpfr_set_inf (e->value.real, -1);
1361 gfc_constructor_append_expr (&result->value.constructor, e,
1366 mpfr_mul_si (e->value.real, x2rev, jn ? (n2-i+1) : (n1+i-1),
1368 mpfr_mul (e->value.real, e->value.real, last2, GFC_RND_MODE);
1369 mpfr_sub (e->value.real, e->value.real, last1, GFC_RND_MODE);
1371 if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
1375 gfc_constructor_insert_expr (&result->value.constructor, e, &x->where,
1378 gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
1380 mpfr_set (last1, last2, GFC_RND_MODE);
1381 mpfr_set (last2, e->value.real, GFC_RND_MODE);
1394 gfc_free_expr (result);
1395 return &gfc_bad_expr;
1400 gfc_simplify_bessel_jn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
1402 return gfc_simplify_bessel_n2 (order1, order2, x, true);
1407 gfc_simplify_bessel_y0 (gfc_expr *x)
1411 if (x->expr_type != EXPR_CONSTANT)
1414 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1415 mpfr_y0 (result->value.real, x->value.real, GFC_RND_MODE);
1417 return range_check (result, "BESSEL_Y0");
1422 gfc_simplify_bessel_y1 (gfc_expr *x)
1426 if (x->expr_type != EXPR_CONSTANT)
1429 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1430 mpfr_y1 (result->value.real, x->value.real, GFC_RND_MODE);
1432 return range_check (result, "BESSEL_Y1");
1437 gfc_simplify_bessel_yn (gfc_expr *order, gfc_expr *x)
1442 if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
1445 n = mpz_get_si (order->value.integer);
1446 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1447 mpfr_yn (result->value.real, n, x->value.real, GFC_RND_MODE);
1449 return range_check (result, "BESSEL_YN");
1454 gfc_simplify_bessel_yn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
1456 return gfc_simplify_bessel_n2 (order1, order2, x, false);
1461 gfc_simplify_bit_size (gfc_expr *e)
1463 int i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
1464 return gfc_get_int_expr (e->ts.kind, &e->where,
1465 gfc_integer_kinds[i].bit_size);
1470 gfc_simplify_btest (gfc_expr *e, gfc_expr *bit)
1474 if (e->expr_type != EXPR_CONSTANT || bit->expr_type != EXPR_CONSTANT)
1477 if (gfc_extract_int (bit, &b) != NULL || b < 0)
1478 return gfc_get_logical_expr (gfc_default_logical_kind, &e->where, false);
1480 return gfc_get_logical_expr (gfc_default_logical_kind, &e->where,
1481 mpz_tstbit (e->value.integer, b));
1486 compare_bitwise (gfc_expr *i, gfc_expr *j)
1491 gcc_assert (i->ts.type == BT_INTEGER);
1492 gcc_assert (j->ts.type == BT_INTEGER);
1494 mpz_init_set (x, i->value.integer);
1495 k = gfc_validate_kind (i->ts.type, i->ts.kind, false);
1496 convert_mpz_to_unsigned (x, gfc_integer_kinds[k].bit_size);
1498 mpz_init_set (y, j->value.integer);
1499 k = gfc_validate_kind (j->ts.type, j->ts.kind, false);
1500 convert_mpz_to_unsigned (y, gfc_integer_kinds[k].bit_size);
1502 res = mpz_cmp (x, y);
1510 gfc_simplify_bge (gfc_expr *i, gfc_expr *j)
1512 if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1515 return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1516 compare_bitwise (i, j) >= 0);
1521 gfc_simplify_bgt (gfc_expr *i, gfc_expr *j)
1523 if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1526 return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1527 compare_bitwise (i, j) > 0);
1532 gfc_simplify_ble (gfc_expr *i, gfc_expr *j)
1534 if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1537 return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1538 compare_bitwise (i, j) <= 0);
1543 gfc_simplify_blt (gfc_expr *i, gfc_expr *j)
1545 if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
1548 return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
1549 compare_bitwise (i, j) < 0);
1554 gfc_simplify_ceiling (gfc_expr *e, gfc_expr *k)
1556 gfc_expr *ceil, *result;
1559 kind = get_kind (BT_INTEGER, k, "CEILING", gfc_default_integer_kind);
1561 return &gfc_bad_expr;
1563 if (e->expr_type != EXPR_CONSTANT)
1566 ceil = gfc_copy_expr (e);
1567 mpfr_ceil (ceil->value.real, e->value.real);
1569 result = gfc_get_constant_expr (BT_INTEGER, kind, &e->where);
1570 gfc_mpfr_to_mpz (result->value.integer, ceil->value.real, &e->where);
1572 gfc_free_expr (ceil);
1574 return range_check (result, "CEILING");
1579 gfc_simplify_char (gfc_expr *e, gfc_expr *k)
1581 return simplify_achar_char (e, k, "CHAR", false);
1585 /* Common subroutine for simplifying CMPLX, COMPLEX and DCMPLX. */
1588 simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
1592 if (convert_boz (x, kind) == &gfc_bad_expr)
1593 return &gfc_bad_expr;
1595 if (convert_boz (y, kind) == &gfc_bad_expr)
1596 return &gfc_bad_expr;
1598 if (x->expr_type != EXPR_CONSTANT
1599 || (y != NULL && y->expr_type != EXPR_CONSTANT))
1602 result = gfc_get_constant_expr (BT_COMPLEX, kind, &x->where);
1607 mpc_set_z (result->value.complex, x->value.integer, GFC_MPC_RND_MODE);
1611 mpc_set_fr (result->value.complex, x->value.real, GFC_RND_MODE);
1615 mpc_set (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1619 gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (x)");
1623 return range_check (result, name);
1628 mpfr_set_z (mpc_imagref (result->value.complex),
1629 y->value.integer, GFC_RND_MODE);
1633 mpfr_set (mpc_imagref (result->value.complex),
1634 y->value.real, GFC_RND_MODE);
1638 gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (y)");
1641 return range_check (result, name);
1646 gfc_simplify_cmplx (gfc_expr *x, gfc_expr *y, gfc_expr *k)
1650 kind = get_kind (BT_REAL, k, "CMPLX", gfc_default_complex_kind);
1652 return &gfc_bad_expr;
1654 return simplify_cmplx ("CMPLX", x, y, kind);
1659 gfc_simplify_complex (gfc_expr *x, gfc_expr *y)
1663 if (x->ts.type == BT_INTEGER && y->ts.type == BT_INTEGER)
1664 kind = gfc_default_complex_kind;
1665 else if (x->ts.type == BT_REAL || y->ts.type == BT_INTEGER)
1667 else if (x->ts.type == BT_INTEGER || y->ts.type == BT_REAL)
1669 else if (x->ts.type == BT_REAL && y->ts.type == BT_REAL)
1670 kind = (x->ts.kind > y->ts.kind) ? x->ts.kind : y->ts.kind;
1674 return simplify_cmplx ("COMPLEX", x, y, kind);
1679 gfc_simplify_conjg (gfc_expr *e)
1683 if (e->expr_type != EXPR_CONSTANT)
1686 result = gfc_copy_expr (e);
1687 mpc_conj (result->value.complex, result->value.complex, GFC_MPC_RND_MODE);
1689 return range_check (result, "CONJG");
1694 gfc_simplify_cos (gfc_expr *x)
1698 if (x->expr_type != EXPR_CONSTANT)
1701 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1706 mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
1710 gfc_set_model_kind (x->ts.kind);
1711 mpc_cos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1715 gfc_internal_error ("in gfc_simplify_cos(): Bad type");
1718 return range_check (result, "COS");
1723 gfc_simplify_cosh (gfc_expr *x)
1727 if (x->expr_type != EXPR_CONSTANT)
1730 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1735 mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
1739 mpc_cosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
1746 return range_check (result, "COSH");
1751 gfc_simplify_count (gfc_expr *mask, gfc_expr *dim, gfc_expr *kind)
1755 if (!is_constant_array_expr (mask)
1756 || !gfc_is_constant_expr (dim)
1757 || !gfc_is_constant_expr (kind))
1760 result = transformational_result (mask, dim,
1762 get_kind (BT_INTEGER, kind, "COUNT",
1763 gfc_default_integer_kind),
1766 init_result_expr (result, 0, NULL);
1768 /* Passing MASK twice, once as data array, once as mask.
1769 Whenever gfc_count is called, '1' is added to the result. */
1770 return !dim || mask->rank == 1 ?
1771 simplify_transformation_to_scalar (result, mask, mask, gfc_count) :
1772 simplify_transformation_to_array (result, mask, dim, mask, gfc_count, NULL);
1777 gfc_simplify_dcmplx (gfc_expr *x, gfc_expr *y)
1779 return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind);
1784 gfc_simplify_dble (gfc_expr *e)
1786 gfc_expr *result = NULL;
1788 if (e->expr_type != EXPR_CONSTANT)
1791 if (convert_boz (e, gfc_default_double_kind) == &gfc_bad_expr)
1792 return &gfc_bad_expr;
1794 result = gfc_convert_constant (e, BT_REAL, gfc_default_double_kind);
1795 if (result == &gfc_bad_expr)
1796 return &gfc_bad_expr;
1798 return range_check (result, "DBLE");
1803 gfc_simplify_digits (gfc_expr *x)
1807 i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
1812 digits = gfc_integer_kinds[i].digits;
1817 digits = gfc_real_kinds[i].digits;
1824 return gfc_get_int_expr (gfc_default_integer_kind, NULL, digits);
1829 gfc_simplify_dim (gfc_expr *x, gfc_expr *y)
1834 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1837 kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
1838 result = gfc_get_constant_expr (x->ts.type, kind, &x->where);
1843 if (mpz_cmp (x->value.integer, y->value.integer) > 0)
1844 mpz_sub (result->value.integer, x->value.integer, y->value.integer);
1846 mpz_set_ui (result->value.integer, 0);
1851 if (mpfr_cmp (x->value.real, y->value.real) > 0)
1852 mpfr_sub (result->value.real, x->value.real, y->value.real,
1855 mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
1860 gfc_internal_error ("gfc_simplify_dim(): Bad type");
1863 return range_check (result, "DIM");
1868 gfc_simplify_dot_product (gfc_expr *vector_a, gfc_expr *vector_b)
1870 if (!is_constant_array_expr (vector_a)
1871 || !is_constant_array_expr (vector_b))
1874 gcc_assert (vector_a->rank == 1);
1875 gcc_assert (vector_b->rank == 1);
1876 gcc_assert (gfc_compare_types (&vector_a->ts, &vector_b->ts));
1878 return compute_dot_product (vector_a, 1, 0, vector_b, 1, 0);
1883 gfc_simplify_dprod (gfc_expr *x, gfc_expr *y)
1885 gfc_expr *a1, *a2, *result;
1887 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1890 a1 = gfc_real2real (x, gfc_default_double_kind);
1891 a2 = gfc_real2real (y, gfc_default_double_kind);
1893 result = gfc_get_constant_expr (BT_REAL, gfc_default_double_kind, &x->where);
1894 mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
1899 return range_check (result, "DPROD");
1904 simplify_dshift (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg,
1908 int i, k, size, shift;
1910 if (arg1->expr_type != EXPR_CONSTANT || arg2->expr_type != EXPR_CONSTANT
1911 || shiftarg->expr_type != EXPR_CONSTANT)
1914 k = gfc_validate_kind (BT_INTEGER, arg1->ts.kind, false);
1915 size = gfc_integer_kinds[k].bit_size;
1917 gfc_extract_int (shiftarg, &shift);
1919 /* DSHIFTR(I,J,SHIFT) = DSHIFTL(I,J,SIZE-SHIFT). */
1921 shift = size - shift;
1923 result = gfc_get_constant_expr (BT_INTEGER, arg1->ts.kind, &arg1->where);
1924 mpz_set_ui (result->value.integer, 0);
1926 for (i = 0; i < shift; i++)
1927 if (mpz_tstbit (arg2->value.integer, size - shift + i))
1928 mpz_setbit (result->value.integer, i);
1930 for (i = 0; i < size - shift; i++)
1931 if (mpz_tstbit (arg1->value.integer, i))
1932 mpz_setbit (result->value.integer, shift + i);
1934 /* Convert to a signed value. */
1935 convert_mpz_to_signed (result->value.integer, size);
1942 gfc_simplify_dshiftr (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg)
1944 return simplify_dshift (arg1, arg2, shiftarg, true);
1949 gfc_simplify_dshiftl (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg)
1951 return simplify_dshift (arg1, arg2, shiftarg, false);
1956 gfc_simplify_erf (gfc_expr *x)
1960 if (x->expr_type != EXPR_CONSTANT)
1963 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1964 mpfr_erf (result->value.real, x->value.real, GFC_RND_MODE);
1966 return range_check (result, "ERF");
1971 gfc_simplify_erfc (gfc_expr *x)
1975 if (x->expr_type != EXPR_CONSTANT)
1978 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
1979 mpfr_erfc (result->value.real, x->value.real, GFC_RND_MODE);
1981 return range_check (result, "ERFC");
1985 /* Helper functions to simplify ERFC_SCALED(x) = ERFC(x) * EXP(X**2). */
1987 #define MAX_ITER 200
1988 #define ARG_LIMIT 12
1990 /* Calculate ERFC_SCALED directly by its definition:
1992 ERFC_SCALED(x) = ERFC(x) * EXP(X**2)
1994 using a large precision for intermediate results. This is used for all
1995 but large values of the argument. */
1997 fullprec_erfc_scaled (mpfr_t res, mpfr_t arg)
2002 prec = mpfr_get_default_prec ();
2003 mpfr_set_default_prec (10 * prec);
2008 mpfr_set (a, arg, GFC_RND_MODE);
2009 mpfr_sqr (b, a, GFC_RND_MODE);
2010 mpfr_exp (b, b, GFC_RND_MODE);
2011 mpfr_erfc (a, a, GFC_RND_MODE);
2012 mpfr_mul (a, a, b, GFC_RND_MODE);
2014 mpfr_set (res, a, GFC_RND_MODE);
2015 mpfr_set_default_prec (prec);
2021 /* Calculate ERFC_SCALED using a power series expansion in 1/arg:
2023 ERFC_SCALED(x) = 1 / (x * sqrt(pi))
2024 * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
2027 This is used for large values of the argument. Intermediate calculations
2028 are performed with twice the precision. We don't do a fixed number of
2029 iterations of the sum, but stop when it has converged to the required
2032 asympt_erfc_scaled (mpfr_t res, mpfr_t arg)
2034 mpfr_t sum, x, u, v, w, oldsum, sumtrunc;
2039 prec = mpfr_get_default_prec ();
2040 mpfr_set_default_prec (2 * prec);
2050 mpfr_init (sumtrunc);
2051 mpfr_set_prec (oldsum, prec);
2052 mpfr_set_prec (sumtrunc, prec);
2054 mpfr_set (x, arg, GFC_RND_MODE);
2055 mpfr_set_ui (sum, 1, GFC_RND_MODE);
2056 mpz_set_ui (num, 1);
2058 mpfr_set (u, x, GFC_RND_MODE);
2059 mpfr_sqr (u, u, GFC_RND_MODE);
2060 mpfr_mul_ui (u, u, 2, GFC_RND_MODE);
2061 mpfr_pow_si (u, u, -1, GFC_RND_MODE);
2063 for (i = 1; i < MAX_ITER; i++)
2065 mpfr_set (oldsum, sum, GFC_RND_MODE);
2067 mpz_mul_ui (num, num, 2 * i - 1);
2070 mpfr_set (w, u, GFC_RND_MODE);
2071 mpfr_pow_ui (w, w, i, GFC_RND_MODE);
2073 mpfr_set_z (v, num, GFC_RND_MODE);
2074 mpfr_mul (v, v, w, GFC_RND_MODE);
2076 mpfr_add (sum, sum, v, GFC_RND_MODE);
2078 mpfr_set (sumtrunc, sum, GFC_RND_MODE);
2079 if (mpfr_cmp (sumtrunc, oldsum) == 0)
2083 /* We should have converged by now; otherwise, ARG_LIMIT is probably
2085 gcc_assert (i < MAX_ITER);
2087 /* Divide by x * sqrt(Pi). */
2088 mpfr_const_pi (u, GFC_RND_MODE);
2089 mpfr_sqrt (u, u, GFC_RND_MODE);
2090 mpfr_mul (u, u, x, GFC_RND_MODE);
2091 mpfr_div (sum, sum, u, GFC_RND_MODE);
2093 mpfr_set (res, sum, GFC_RND_MODE);
2094 mpfr_set_default_prec (prec);
2096 mpfr_clears (sum, x, u, v, w, oldsum, sumtrunc, NULL);
2102 gfc_simplify_erfc_scaled (gfc_expr *x)
2106 if (x->expr_type != EXPR_CONSTANT)
2109 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2110 if (mpfr_cmp_d (x->value.real, ARG_LIMIT) >= 0)
2111 asympt_erfc_scaled (result->value.real, x->value.real);
2113 fullprec_erfc_scaled (result->value.real, x->value.real);
2115 return range_check (result, "ERFC_SCALED");
2123 gfc_simplify_epsilon (gfc_expr *e)
2128 i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2130 result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
2131 mpfr_set (result->value.real, gfc_real_kinds[i].epsilon, GFC_RND_MODE);
2133 return range_check (result, "EPSILON");
2138 gfc_simplify_exp (gfc_expr *x)
2142 if (x->expr_type != EXPR_CONSTANT)
2145 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2150 mpfr_exp (result->value.real, x->value.real, GFC_RND_MODE);
2154 gfc_set_model_kind (x->ts.kind);
2155 mpc_exp (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
2159 gfc_internal_error ("in gfc_simplify_exp(): Bad type");
2162 return range_check (result, "EXP");
2167 gfc_simplify_exponent (gfc_expr *x)
2172 if (x->expr_type != EXPR_CONSTANT)
2175 result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
2178 gfc_set_model (x->value.real);
2180 if (mpfr_sgn (x->value.real) == 0)
2182 mpz_set_ui (result->value.integer, 0);
2186 i = (int) mpfr_get_exp (x->value.real);
2187 mpz_set_si (result->value.integer, i);
2189 return range_check (result, "EXPONENT");
2194 gfc_simplify_float (gfc_expr *a)
2198 if (a->expr_type != EXPR_CONSTANT)
2203 if (convert_boz (a, gfc_default_real_kind) == &gfc_bad_expr)
2204 return &gfc_bad_expr;
2206 result = gfc_copy_expr (a);
2209 result = gfc_int2real (a, gfc_default_real_kind);
2211 return range_check (result, "FLOAT");
2216 is_last_ref_vtab (gfc_expr *e)
2219 gfc_component *comp = NULL;
2221 if (e->expr_type != EXPR_VARIABLE)
2224 for (ref = e->ref; ref; ref = ref->next)
2225 if (ref->type == REF_COMPONENT)
2226 comp = ref->u.c.component;
2228 if (!e->ref || !comp)
2229 return e->symtree->n.sym->attr.vtab;
2231 if (comp->name[0] == '_' && strcmp (comp->name, "_vptr") == 0)
2239 gfc_simplify_extends_type_of (gfc_expr *a, gfc_expr *mold)
2241 /* Avoid simplification of resolved symbols. */
2242 if (is_last_ref_vtab (a) || is_last_ref_vtab (mold))
2245 if (a->ts.type == BT_DERIVED && mold->ts.type == BT_DERIVED)
2246 return gfc_get_logical_expr (gfc_default_logical_kind, &a->where,
2247 gfc_type_is_extension_of (mold->ts.u.derived,
2249 /* Return .false. if the dynamic type can never be the same. */
2250 if ((a->ts.type == BT_CLASS && mold->ts.type == BT_CLASS
2251 && !gfc_type_is_extension_of
2252 (mold->ts.u.derived->components->ts.u.derived,
2253 a->ts.u.derived->components->ts.u.derived)
2254 && !gfc_type_is_extension_of
2255 (a->ts.u.derived->components->ts.u.derived,
2256 mold->ts.u.derived->components->ts.u.derived))
2257 || (a->ts.type == BT_DERIVED && mold->ts.type == BT_CLASS
2258 && !gfc_type_is_extension_of
2260 mold->ts.u.derived->components->ts.u.derived)
2261 && !gfc_type_is_extension_of
2262 (mold->ts.u.derived->components->ts.u.derived,
2264 || (a->ts.type == BT_CLASS && mold->ts.type == BT_DERIVED
2265 && !gfc_type_is_extension_of
2266 (mold->ts.u.derived,
2267 a->ts.u.derived->components->ts.u.derived)))
2268 return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, false);
2270 if (mold->ts.type == BT_DERIVED
2271 && gfc_type_is_extension_of (mold->ts.u.derived,
2272 a->ts.u.derived->components->ts.u.derived))
2273 return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, true);
2280 gfc_simplify_same_type_as (gfc_expr *a, gfc_expr *b)
2282 /* Avoid simplification of resolved symbols. */
2283 if (is_last_ref_vtab (a) || is_last_ref_vtab (b))
2286 /* Return .false. if the dynamic type can never be the
2288 if ((a->ts.type == BT_CLASS || b->ts.type == BT_CLASS)
2289 && !gfc_type_compatible (&a->ts, &b->ts)
2290 && !gfc_type_compatible (&b->ts, &a->ts))
2291 return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, false);
2293 if (a->ts.type != BT_DERIVED || b->ts.type != BT_DERIVED)
2296 return gfc_get_logical_expr (gfc_default_logical_kind, &a->where,
2297 gfc_compare_derived_types (a->ts.u.derived,
2303 gfc_simplify_floor (gfc_expr *e, gfc_expr *k)
2309 kind = get_kind (BT_INTEGER, k, "FLOOR", gfc_default_integer_kind);
2311 gfc_internal_error ("gfc_simplify_floor(): Bad kind");
2313 if (e->expr_type != EXPR_CONSTANT)
2316 gfc_set_model_kind (kind);
2319 mpfr_floor (floor, e->value.real);
2321 result = gfc_get_constant_expr (BT_INTEGER, kind, &e->where);
2322 gfc_mpfr_to_mpz (result->value.integer, floor, &e->where);
2326 return range_check (result, "FLOOR");
2331 gfc_simplify_fraction (gfc_expr *x)
2334 mpfr_t absv, exp, pow2;
2336 if (x->expr_type != EXPR_CONSTANT)
2339 result = gfc_get_constant_expr (BT_REAL, x->ts.kind, &x->where);
2341 if (mpfr_sgn (x->value.real) == 0)
2343 mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2347 gfc_set_model_kind (x->ts.kind);
2352 mpfr_abs (absv, x->value.real, GFC_RND_MODE);
2353 mpfr_log2 (exp, absv, GFC_RND_MODE);
2355 mpfr_trunc (exp, exp);
2356 mpfr_add_ui (exp, exp, 1, GFC_RND_MODE);
2358 mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
2360 mpfr_div (result->value.real, absv, pow2, GFC_RND_MODE);
2362 mpfr_clears (exp, absv, pow2, NULL);
2364 return range_check (result, "FRACTION");
2369 gfc_simplify_gamma (gfc_expr *x)
2373 if (x->expr_type != EXPR_CONSTANT)
2376 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2377 mpfr_gamma (result->value.real, x->value.real, GFC_RND_MODE);
2379 return range_check (result, "GAMMA");
2384 gfc_simplify_huge (gfc_expr *e)
2389 i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2390 result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);
2395 mpz_set (result->value.integer, gfc_integer_kinds[i].huge);
2399 mpfr_set (result->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
2411 gfc_simplify_hypot (gfc_expr *x, gfc_expr *y)
2415 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2418 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2419 mpfr_hypot (result->value.real, x->value.real, y->value.real, GFC_RND_MODE);
2420 return range_check (result, "HYPOT");
2424 /* We use the processor's collating sequence, because all
2425 systems that gfortran currently works on are ASCII. */
2428 gfc_simplify_iachar (gfc_expr *e, gfc_expr *kind)
2434 if (e->expr_type != EXPR_CONSTANT)
2437 if (e->value.character.length != 1)
2439 gfc_error ("Argument of IACHAR at %L must be of length one", &e->where);
2440 return &gfc_bad_expr;
2443 index = e->value.character.string[0];
2445 if (gfc_option.warn_surprising && index > 127)
2446 gfc_warning ("Argument of IACHAR function at %L outside of range 0..127",
2449 k = get_kind (BT_INTEGER, kind, "IACHAR", gfc_default_integer_kind);
2451 return &gfc_bad_expr;
2453 result = gfc_get_int_expr (k, &e->where, index);
2455 return range_check (result, "IACHAR");
2460 do_bit_and (gfc_expr *result, gfc_expr *e)
2462 gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
2463 gcc_assert (result->ts.type == BT_INTEGER
2464 && result->expr_type == EXPR_CONSTANT);
2466 mpz_and (result->value.integer, result->value.integer, e->value.integer);
2472 gfc_simplify_iall (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
2474 return simplify_transformation (array, dim, mask, -1, do_bit_and);
2479 do_bit_ior (gfc_expr *result, gfc_expr *e)
2481 gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
2482 gcc_assert (result->ts.type == BT_INTEGER
2483 && result->expr_type == EXPR_CONSTANT);
2485 mpz_ior (result->value.integer, result->value.integer, e->value.integer);
2491 gfc_simplify_iany (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
2493 return simplify_transformation (array, dim, mask, 0, do_bit_ior);
2498 gfc_simplify_iand (gfc_expr *x, gfc_expr *y)
2502 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2505 result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
2506 mpz_and (result->value.integer, x->value.integer, y->value.integer);
2508 return range_check (result, "IAND");
2513 gfc_simplify_ibclr (gfc_expr *x, gfc_expr *y)
2518 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2521 gfc_extract_int (y, &pos);
2523 k = gfc_validate_kind (x->ts.type, x->ts.kind, false);
2525 result = gfc_copy_expr (x);
2527 convert_mpz_to_unsigned (result->value.integer,
2528 gfc_integer_kinds[k].bit_size);
2530 mpz_clrbit (result->value.integer, pos);
2532 convert_mpz_to_signed (result->value.integer,
2533 gfc_integer_kinds[k].bit_size);
2540 gfc_simplify_ibits (gfc_expr *x, gfc_expr *y, gfc_expr *z)
2547 if (x->expr_type != EXPR_CONSTANT
2548 || y->expr_type != EXPR_CONSTANT
2549 || z->expr_type != EXPR_CONSTANT)
2552 gfc_extract_int (y, &pos);
2553 gfc_extract_int (z, &len);
2555 k = gfc_validate_kind (BT_INTEGER, x->ts.kind, false);
2557 bitsize = gfc_integer_kinds[k].bit_size;
2559 if (pos + len > bitsize)
2561 gfc_error ("Sum of second and third arguments of IBITS exceeds "
2562 "bit size at %L", &y->where);
2563 return &gfc_bad_expr;
2566 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
2567 convert_mpz_to_unsigned (result->value.integer,
2568 gfc_integer_kinds[k].bit_size);
2570 bits = XCNEWVEC (int, bitsize);
2572 for (i = 0; i < bitsize; i++)
2575 for (i = 0; i < len; i++)
2576 bits[i] = mpz_tstbit (x->value.integer, i + pos);
2578 for (i = 0; i < bitsize; i++)
2581 mpz_clrbit (result->value.integer, i);
2582 else if (bits[i] == 1)
2583 mpz_setbit (result->value.integer, i);
2585 gfc_internal_error ("IBITS: Bad bit");
2590 convert_mpz_to_signed (result->value.integer,
2591 gfc_integer_kinds[k].bit_size);
2598 gfc_simplify_ibset (gfc_expr *x, gfc_expr *y)
2603 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2606 gfc_extract_int (y, &pos);
2608 k = gfc_validate_kind (x->ts.type, x->ts.kind, false);
2610 result = gfc_copy_expr (x);
2612 convert_mpz_to_unsigned (result->value.integer,
2613 gfc_integer_kinds[k].bit_size);
2615 mpz_setbit (result->value.integer, pos);
2617 convert_mpz_to_signed (result->value.integer,
2618 gfc_integer_kinds[k].bit_size);
2625 gfc_simplify_ichar (gfc_expr *e, gfc_expr *kind)
2631 if (e->expr_type != EXPR_CONSTANT)
2634 if (e->value.character.length != 1)
2636 gfc_error ("Argument of ICHAR at %L must be of length one", &e->where);
2637 return &gfc_bad_expr;
2640 index = e->value.character.string[0];
2642 k = get_kind (BT_INTEGER, kind, "ICHAR", gfc_default_integer_kind);
2644 return &gfc_bad_expr;
2646 result = gfc_get_int_expr (k, &e->where, index);
2648 return range_check (result, "ICHAR");
2653 gfc_simplify_ieor (gfc_expr *x, gfc_expr *y)
2657 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2660 result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
2661 mpz_xor (result->value.integer, x->value.integer, y->value.integer);
2663 return range_check (result, "IEOR");
2668 gfc_simplify_index (gfc_expr *x, gfc_expr *y, gfc_expr *b, gfc_expr *kind)
2671 int back, len, lensub;
2672 int i, j, k, count, index = 0, start;
2674 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT
2675 || ( b != NULL && b->expr_type != EXPR_CONSTANT))
2678 if (b != NULL && b->value.logical != 0)
2683 k = get_kind (BT_INTEGER, kind, "INDEX", gfc_default_integer_kind);
2685 return &gfc_bad_expr;
2687 result = gfc_get_constant_expr (BT_INTEGER, k, &x->where);
2689 len = x->value.character.length;
2690 lensub = y->value.character.length;
2694 mpz_set_si (result->value.integer, 0);
2702 mpz_set_si (result->value.integer, 1);
2705 else if (lensub == 1)
2707 for (i = 0; i < len; i++)
2709 for (j = 0; j < lensub; j++)
2711 if (y->value.character.string[j]
2712 == x->value.character.string[i])
2722 for (i = 0; i < len; i++)
2724 for (j = 0; j < lensub; j++)
2726 if (y->value.character.string[j]
2727 == x->value.character.string[i])
2732 for (k = 0; k < lensub; k++)
2734 if (y->value.character.string[k]
2735 == x->value.character.string[k + start])
2739 if (count == lensub)
2754 mpz_set_si (result->value.integer, len + 1);
2757 else if (lensub == 1)
2759 for (i = 0; i < len; i++)
2761 for (j = 0; j < lensub; j++)
2763 if (y->value.character.string[j]
2764 == x->value.character.string[len - i])
2766 index = len - i + 1;
2774 for (i = 0; i < len; i++)
2776 for (j = 0; j < lensub; j++)
2778 if (y->value.character.string[j]
2779 == x->value.character.string[len - i])
2782 if (start <= len - lensub)
2785 for (k = 0; k < lensub; k++)
2786 if (y->value.character.string[k]
2787 == x->value.character.string[k + start])
2790 if (count == lensub)
2807 mpz_set_si (result->value.integer, index);
2808 return range_check (result, "INDEX");
2813 simplify_intconv (gfc_expr *e, int kind, const char *name)
2815 gfc_expr *result = NULL;
2817 if (e->expr_type != EXPR_CONSTANT)
2820 result = gfc_convert_constant (e, BT_INTEGER, kind);
2821 if (result == &gfc_bad_expr)
2822 return &gfc_bad_expr;
2824 return range_check (result, name);
2829 gfc_simplify_int (gfc_expr *e, gfc_expr *k)
2833 kind = get_kind (BT_INTEGER, k, "INT", gfc_default_integer_kind);
2835 return &gfc_bad_expr;
2837 return simplify_intconv (e, kind, "INT");
2841 gfc_simplify_int2 (gfc_expr *e)
2843 return simplify_intconv (e, 2, "INT2");
2848 gfc_simplify_int8 (gfc_expr *e)
2850 return simplify_intconv (e, 8, "INT8");
2855 gfc_simplify_long (gfc_expr *e)
2857 return simplify_intconv (e, 4, "LONG");
2862 gfc_simplify_ifix (gfc_expr *e)
2864 gfc_expr *rtrunc, *result;
2866 if (e->expr_type != EXPR_CONSTANT)
2869 rtrunc = gfc_copy_expr (e);
2870 mpfr_trunc (rtrunc->value.real, e->value.real);
2872 result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
2874 gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real, &e->where);
2876 gfc_free_expr (rtrunc);
2878 return range_check (result, "IFIX");
2883 gfc_simplify_idint (gfc_expr *e)
2885 gfc_expr *rtrunc, *result;
2887 if (e->expr_type != EXPR_CONSTANT)
2890 rtrunc = gfc_copy_expr (e);
2891 mpfr_trunc (rtrunc->value.real, e->value.real);
2893 result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
2895 gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real, &e->where);
2897 gfc_free_expr (rtrunc);
2899 return range_check (result, "IDINT");
2904 gfc_simplify_ior (gfc_expr *x, gfc_expr *y)
2908 if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2911 result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
2912 mpz_ior (result->value.integer, x->value.integer, y->value.integer);
2914 return range_check (result, "IOR");
2919 do_bit_xor (gfc_expr *result, gfc_expr *e)
2921 gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
2922 gcc_assert (result->ts.type == BT_INTEGER
2923 && result->expr_type == EXPR_CONSTANT);
2925 mpz_xor (result->value.integer, result->value.integer, e->value.integer);
2931 gfc_simplify_iparity (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
2933 return simplify_transformation (array, dim, mask, 0, do_bit_xor);
2939 gfc_simplify_is_iostat_end (gfc_expr *x)
2941 if (x->expr_type != EXPR_CONSTANT)
2944 return gfc_get_logical_expr (gfc_default_logical_kind, &x->where,
2945 mpz_cmp_si (x->value.integer,
2946 LIBERROR_END) == 0);
2951 gfc_simplify_is_iostat_eor (gfc_expr *x)
2953 if (x->expr_type != EXPR_CONSTANT)
2956 return gfc_get_logical_expr (gfc_default_logical_kind, &x->where,
2957 mpz_cmp_si (x->value.integer,
2958 LIBERROR_EOR) == 0);
2963 gfc_simplify_isnan (gfc_expr *x)
2965 if (x->expr_type != EXPR_CONSTANT)
2968 return gfc_get_logical_expr (gfc_default_logical_kind, &x->where,
2969 mpfr_nan_p (x->value.real));
2973 /* Performs a shift on its first argument. Depending on the last
2974 argument, the shift can be arithmetic, i.e. with filling from the
2975 left like in the SHIFTA intrinsic. */
2977 simplify_shift (gfc_expr *e, gfc_expr *s, const char *name,
2978 bool arithmetic, int direction)
2981 int ashift, *bits, i, k, bitsize, shift;
2983 if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
2986 gfc_extract_int (s, &shift);
2988 k = gfc_validate_kind (BT_INTEGER, e->ts.kind, false);
2989 bitsize = gfc_integer_kinds[k].bit_size;
2991 result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);
2995 mpz_set (result->value.integer, e->value.integer);
2999 if (direction > 0 && shift < 0)
3001 /* Left shift, as in SHIFTL. */
3002 gfc_error ("Second argument of %s is negative at %L", name, &e->where);
3003 return &gfc_bad_expr;
3005 else if (direction < 0)
3007 /* Right shift, as in SHIFTR or SHIFTA. */
3010 gfc_error ("Second argument of %s is negative at %L",
3012 return &gfc_bad_expr;
3018 ashift = (shift >= 0 ? shift : -shift);
3020 if (ashift > bitsize)
3022 gfc_error ("Magnitude of second argument of %s exceeds bit size "
3023 "at %L", name, &e->where);
3024 return &gfc_bad_expr;
3027 bits = XCNEWVEC (int, bitsize);
3029 for (i = 0; i < bitsize; i++)
3030 bits[i] = mpz_tstbit (e->value.integer, i);
3035 for (i = 0; i < shift; i++)
3036 mpz_clrbit (result->value.integer, i);
3038 for (i = 0; i < bitsize - shift; i++)
3041 mpz_clrbit (result->value.integer, i + shift);
3043 mpz_setbit (result->value.integer, i + shift);
3049 if (arithmetic && bits[bitsize - 1])
3050 for (i = bitsize - 1; i >= bitsize - ashift; i--)
3051 mpz_setbit (result->value.integer, i);
3053 for (i = bitsize - 1; i >= bitsize - ashift; i--)
3054 mpz_clrbit (result->value.integer, i);
3056 for (i = bitsize - 1; i >= ashift; i--)
3059 mpz_clrbit (result->value.integer, i - ashift);
3061 mpz_setbit (result->value.integer, i - ashift);
3065 convert_mpz_to_signed (result->value.integer, bitsize);
3073 gfc_simplify_ishft (gfc_expr *e, gfc_expr *s)
3075 return simplify_shift (e, s, "ISHFT", false, 0);
3080 gfc_simplify_lshift (gfc_expr *e, gfc_expr *s)
3082 return simplify_shift (e, s, "LSHIFT", false, 1);
3087 gfc_simplify_rshift (gfc_expr *e, gfc_expr *s)
3089 return simplify_shift (e, s, "RSHIFT", true, -1);
3094 gfc_simplify_shifta (gfc_expr *e, gfc_expr *s)
3096 return simplify_shift (e, s, "SHIFTA", true, -1);
3101 gfc_simplify_shiftl (gfc_expr *e, gfc_expr *s)
3103 return simplify_shift (e, s, "SHIFTL", false, 1);
3108 gfc_simplify_shiftr (gfc_expr *e, gfc_expr *s)
3110 return simplify_shift (e, s, "SHIFTR", false, -1);
3115 gfc_simplify_ishftc (gfc_expr *e, gfc_expr *s, gfc_expr *sz)
3118 int shift, ashift, isize, ssize, delta, k;
3121 if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
3124 gfc_extract_int (s, &shift);
3126 k = gfc_validate_kind (e->ts.type, e->ts.kind, false);
3127 isize = gfc_integer_kinds[k].bit_size;
3131 if (sz->expr_type != EXPR_CONSTANT)
3134 gfc_extract_int (sz, &ssize);
3148 gfc_error ("Magnitude of second argument of ISHFTC exceeds "
3149 "BIT_SIZE of first argument at %L", &s->where);
3150 return &gfc_bad_expr;
3153 result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);
3155 mpz_set (result->value.integer, e->value.integer);
3160 convert_mpz_to_unsigned (result->value.integer, isize);
3162 bits = XCNEWVEC (int, ssize);
3164 for (i = 0; i < ssize; i++)
3165 bits[i] = mpz_tstbit (e->value.integer, i);
3167 delta = ssize - ashift;
3171 for (i = 0; i < delta; i++)
3174 mpz_clrbit (result->value.integer, i + shift);
3176 mpz_setbit (result->value.integer, i + shift);
3179 for (i = delta; i < ssize; i++)
3182 mpz_clrbit (result->value.integer, i - delta);
3184 mpz_setbit (result->value.integer, i - delta);
3189 for (i = 0; i < ashift; i++)
3192 mpz_clrbit (result->value.integer, i + delta);
3194 mpz_setbit (result->value.integer, i + delta);
3197 for (i = ashift; i < ssize; i++)
3200 mpz_clrbit (result->value.integer, i + shift);
3202 mpz_setbit (result->value.integer, i + shift);
3206 convert_mpz_to_signed (result->value.integer, isize);
3214 gfc_simplify_kind (gfc_expr *e)
3216 return gfc_get_int_expr (gfc_default_integer_kind, NULL, e->ts.kind);
3221 simplify_bound_dim (gfc_expr *array, gfc_expr *kind, int d, int upper,
3222 gfc_array_spec *as, gfc_ref *ref, bool coarray)
3224 gfc_expr *l, *u, *result;
3227 k = get_kind (BT_INTEGER, kind, upper ? "UBOUND" : "LBOUND",
3228 gfc_default_integer_kind);
3230 return &gfc_bad_expr;
3232 result = gfc_get_constant_expr (BT_INTEGER, k, &array->where);
3234 /* For non-variables, LBOUND(expr, DIM=n) = 1 and
3235 UBOUND(expr, DIM=n) = SIZE(expr, DIM=n). */
3236 if (!coarray && array->expr_type != EXPR_VARIABLE)
3240 gfc_expr* dim = result;
3241 mpz_set_si (dim->value.integer, d);
3243 result = gfc_simplify_size (array, dim, kind);
3244 gfc_free_expr (dim);
3249 mpz_set_si (result->value.integer, 1);
3254 /* Otherwise, we have a variable expression. */
3255 gcc_assert (array->expr_type == EXPR_VARIABLE);
3258 if (gfc_resolve_array_spec (as, 0) == FAILURE)
3261 /* The last dimension of an assumed-size array is special. */
3262 if ((!coarray && d == as->rank && as->type == AS_ASSUMED_SIZE && !upper)
3263 || (coarray && d == as->rank + as->corank
3264 && (!upper || gfc_option.coarray == GFC_FCOARRAY_SINGLE)))
3266 if (as->lower[d-1]->expr_type == EXPR_CONSTANT)
3268 gfc_free_expr (result);
3269 return gfc_copy_expr (as->lower[d-1]);
3275 result = gfc_get_constant_expr (BT_INTEGER, k, &array->where);
3277 /* Then, we need to know the extent of the given dimension. */
3278 if (coarray || ref->u.ar.type == AR_FULL)
3283 if (l->expr_type != EXPR_CONSTANT || u == NULL
3284 || u->expr_type != EXPR_CONSTANT)
3287 if (mpz_cmp (l->value.integer, u->value.integer) > 0)
3291 mpz_set_si (result->value.integer, 0);
3293 mpz_set_si (result->value.integer, 1);
3297 /* Nonzero extent. */
3299 mpz_set (result->value.integer, u->value.integer);
3301 mpz_set (result->value.integer, l->value.integer);
3308 if (gfc_ref_dimen_size (&ref->u.ar, d-1, &result->value.integer, NULL)
3313 mpz_set_si (result->value.integer, (long int) 1);
3317 return range_check (result, upper ? "UBOUND" : "LBOUND");
3320 gfc_free_expr (result);
3326 simplify_bound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
3332 if (array->ts.type == BT_CLASS)
3335 if (array->expr_type != EXPR_VARIABLE)
3342 /* Follow any component references. */
3343 as = array->symtree->n.sym->as;
3344 for (ref = array->ref; ref; ref = ref->next)
3349 switch (ref->u.ar.type)
3356 /* We're done because 'as' has already been set in the
3357 previous iteration. */
3374 as = ref->u.c.component->as;
3386 if (as && (as->type == AS_DEFERRED || as->type == AS_ASSUMED_SHAPE))
3391 /* Multi-dimensional bounds. */
3392 gfc_expr *bounds[GFC_MAX_DIMENSIONS];
3396 /* UBOUND(ARRAY) is not valid for an assumed-size array. */
3397 if (upper && as && as->type == AS_ASSUMED_SIZE)
3399 /* An error message will be emitted in
3400 check_assumed_size_reference (resolve.c). */
3401 return &gfc_bad_expr;
3404 /* Simplify the bounds for each dimension. */
3405 for (d = 0; d < array->rank; d++)
3407 bounds[d] = simplify_bound_dim (array, kind, d + 1, upper, as, ref,
3409 if (bounds[d] == NULL || bounds[d] == &gfc_bad_expr)
3413 for (j = 0; j < d; j++)
3414 gfc_free_expr (bounds[j]);
3419 /* Allocate the result expression. */
3420 k = get_kind (BT_INTEGER, kind, upper ? "UBOUND" : "LBOUND",
3421 gfc_default_integer_kind);
3423 return &gfc_bad_expr;
3425 e = gfc_get_array_expr (BT_INTEGER, k, &array->where);
3427 /* The result is a rank 1 array; its size is the rank of the first
3428 argument to {L,U}BOUND. */
3430 e->shape = gfc_get_shape (1);
3431 mpz_init_set_ui (e->shape[0], array->rank);
3433 /* Create the constructor for this array. */
3434 for (d = 0; d < array->rank; d++)
3435 gfc_constructor_append_expr (&e->value.constructor,
3436 bounds[d], &e->where);
3442 /* A DIM argument is specified. */
3443 if (dim->expr_type != EXPR_CONSTANT)
3446 d = mpz_get_si (dim->value.integer);
3448 if (d < 1 || d > array->rank
3449 || (d == array->rank && as && as->type == AS_ASSUMED_SIZE && upper))
3451 gfc_error ("DIM argument at %L is out of bounds", &dim->where);
3452 return &gfc_bad_expr;
3455 return simplify_bound_dim (array, kind, d, upper, as, ref, false);
3461 simplify_cobound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
3467 if (array->expr_type != EXPR_VARIABLE)
3470 /* Follow any component references. */
3471 as = (array->ts.type == BT_CLASS && array->ts.u.derived->components)
3472 ? array->ts.u.derived->components->as
3473 : array->symtree->n.sym->as;
3474 for (ref = array->ref; ref; ref = ref->next)
3479 switch (ref->u.ar.type)
3482 if (ref->u.ar.as->corank > 0)
3484 gcc_assert (as == ref->u.ar.as);
3491 /* We're done because 'as' has already been set in the
3492 previous iteration. */
3509 as = ref->u.c.component->as;
3522 if (as->cotype == AS_DEFERRED || as->cotype == AS_ASSUMED_SHAPE)
3527 /* Multi-dimensional cobounds. */
3528 gfc_expr *bounds[GFC_MAX_DIMENSIONS];
3532 /* Simplify the cobounds for each dimension. */
3533 for (d = 0; d < as->corank; d++)
3535 bounds[d] = simplify_bound_dim (array, kind, d + 1 + as->rank,
3536 upper, as, ref, true);
3537 if (bounds[d] == NULL || bounds[d] == &gfc_bad_expr)
3541 for (j = 0; j < d; j++)
3542 gfc_free_expr (bounds[j]);
3547 /* Allocate the result expression. */
3548 e = gfc_get_expr ();
3549 e->where = array->where;
3550 e->expr_type = EXPR_ARRAY;
3551 e->ts.type = BT_INTEGER;
3552 k = get_kind (BT_INTEGER, kind, upper ? "UCOBOUND" : "LCOBOUND",
3553 gfc_default_integer_kind);
3557 return &gfc_bad_expr;
3561 /* The result is a rank 1 array; its size is the rank of the first
3562 argument to {L,U}COBOUND. */
3564 e->shape = gfc_get_shape (1);
3565 mpz_init_set_ui (e->shape[0], as->corank);
3567 /* Create the constructor for this array. */
3568 for (d = 0; d < as->corank; d++)
3569 gfc_constructor_append_expr (&e->value.constructor,
3570 bounds[d], &e->where);
3575 /* A DIM argument is specified. */
3576 if (dim->expr_type != EXPR_CONSTANT)
3579 d = mpz_get_si (dim->value.integer);
3581 if (d < 1 || d > as->corank)
3583 gfc_error ("DIM argument at %L is out of bounds", &dim->where);
3584 return &gfc_bad_expr;
3587 return simplify_bound_dim (array, kind, d+as->rank, upper, as, ref, true);
3593 gfc_simplify_lbound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind)
3595 return simplify_bound (array, dim, kind, 0);
3600 gfc_simplify_lcobound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind)
3602 return simplify_cobound (array, dim, kind, 0);
3606 gfc_simplify_leadz (gfc_expr *e)
3608 unsigned long lz, bs;
3611 if (e->expr_type != EXPR_CONSTANT)
3614 i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
3615 bs = gfc_integer_kinds[i].bit_size;
3616 if (mpz_cmp_si (e->value.integer, 0) == 0)
3618 else if (mpz_cmp_si (e->value.integer, 0) < 0)
3621 lz = bs - mpz_sizeinbase (e->value.integer, 2);
3623 return gfc_get_int_expr (gfc_default_integer_kind, &e->where, lz);
3628 gfc_simplify_len (gfc_expr *e, gfc_expr *kind)
3631 int k = get_kind (BT_INTEGER, kind, "LEN", gfc_default_integer_kind);
3634 return &gfc_bad_expr;
3636 if (e->expr_type == EXPR_CONSTANT)
3638 result = gfc_get_constant_expr (BT_INTEGER, k, &e->where);
3639 mpz_set_si (result->value.integer, e->value.character.length);
3640 return range_check (result, "LEN");
3642 else if (e->ts.u.cl != NULL && e->ts.u.cl->length != NULL
3643 && e->ts.u.cl->length->expr_type == EXPR_CONSTANT
3644 && e->ts.u.cl->length->ts.type == BT_INTEGER)
3646 result = gfc_get_constant_expr (BT_INTEGER, k, &e->where);
3647 mpz_set (result->value.integer, e->ts.u.cl->length->value.integer);
3648 return range_check (result, "LEN");
3656 gfc_simplify_len_trim (gfc_expr *e, gfc_expr *kind)
3660 int k = get_kind (BT_INTEGER, kind, "LEN_TRIM", gfc_default_integer_kind);
3663 return &gfc_bad_expr;
3665 if (e->expr_type != EXPR_CONSTANT)
3668 len = e->value.character.length;
3669 for (count = 0, i = 1; i <= len; i++)
3670 if (e->value.character.string[len - i] == ' ')
3675 result = gfc_get_int_expr (k, &e->where, len - count);
3676 return range_check (result, "LEN_TRIM");
3680 gfc_simplify_lgamma (gfc_expr *x)
3685 if (x->expr_type != EXPR_CONSTANT)
3688 result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
3689 mpfr_lgamma (result->value.real, &sg, x->value.real, GFC_RND_MODE);
3691 return range_check (result, "LGAMMA");
3696 gfc_simplify_lge (gfc_expr *a, gfc_expr *b)
3698 if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
3701 return gfc_get_logical_expr (gfc_default_logical_kind, &a-