diff options
| -rw-r--r-- | gcc/fortran/ChangeLog | 14 | ||||
| -rw-r--r-- | gcc/fortran/arith.c | 285 | ||||
| -rw-r--r-- | gcc/fortran/expr.c | 25 | ||||
| -rw-r--r-- | gcc/fortran/gfortran.h | 2 | ||||
| -rw-r--r-- | gcc/testsuite/ChangeLog | 5 | ||||
| -rw-r--r-- | gcc/testsuite/gfortran.dg/power1.f90 | 58 |
6 files changed, 275 insertions, 114 deletions
diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index 3156df6..8bdf010 100644 --- a/gcc/fortran/ChangeLog +++ b/gcc/fortran/ChangeLog @@ -1,3 +1,17 @@ +2009-03-29 Steven G. Kargl <kargl@gcc.gnu.org> + + PR fortran/38823 + * gfortran.h: Add ARITH_PROHIBIT to arith enum. + expr.c (gfc_match_init_expr): Add global variable init_flag to + flag matching an initialization expression. + (check_intrinsic_op): Move no longer reachable error message to ... + * arith.c (arith_power): ... here. Remove gfc_ prefix in + gfc_arith_power. Use init_flag. Allow constant folding of x**y + when y is REAL or COMPLEX. + (eval_intrinsic): Remove restriction that y in x**y must be INTEGER + for constant folding. + * gfc_power: Update gfc_arith_power to arith_power + 2009-03-29 Daniel Kraft <d@domob.eu> PR fortran/37423 diff --git a/gcc/fortran/arith.c b/gcc/fortran/arith.c index 7440be3..17f2221 100644 --- a/gcc/fortran/arith.c +++ b/gcc/fortran/arith.c @@ -932,131 +932,213 @@ complex_pow (gfc_expr *result, gfc_expr *base, mpz_t power) } -/* Raise a number to an integer power. */ +/* Raise a number to a power. */ static arith -gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp) +arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp) { int power_sign; gfc_expr *result; arith rc; - - gcc_assert (op2->expr_type == EXPR_CONSTANT && op2->ts.type == BT_INTEGER); + extern bool init_flag; rc = ARITH_OK; result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where); - power_sign = mpz_sgn (op2->value.integer); - if (power_sign == 0) + switch (op2->ts.type) { - /* Handle something to the zeroth power. Since we're dealing - with integral exponents, there is no ambiguity in the - limiting procedure used to determine the value of 0**0. */ - switch (op1->ts.type) + case BT_INTEGER: + power_sign = mpz_sgn (op2->value.integer); + + if (power_sign == 0) { - case BT_INTEGER: - mpz_set_ui (result->value.integer, 1); - break; + /* Handle something to the zeroth power. Since we're dealing + with integral exponents, there is no ambiguity in the + limiting procedure used to determine the value of 0**0. */ + switch (op1->ts.type) + { + case BT_INTEGER: + mpz_set_ui (result->value.integer, 1); + break; - case BT_REAL: - mpfr_set_ui (result->value.real, 1, GFC_RND_MODE); - break; + case BT_REAL: + mpfr_set_ui (result->value.real, 1, GFC_RND_MODE); + break; - case BT_COMPLEX: - mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE); - mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE); - break; + case BT_COMPLEX: + mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE); + mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE); + break; - default: - gfc_internal_error ("gfc_arith_power(): Bad base"); + default: + gfc_internal_error ("arith_power(): Bad base"); + } } - } - else - { - switch (op1->ts.type) + else { - case BT_INTEGER: - { - int power; - - /* First, we simplify the cases of op1 == 1, 0 or -1. */ - if (mpz_cmp_si (op1->value.integer, 1) == 0) - { - /* 1**op2 == 1 */ - mpz_set_si (result->value.integer, 1); - } - else if (mpz_cmp_si (op1->value.integer, 0) == 0) - { - /* 0**op2 == 0, if op2 > 0 - 0**op2 overflow, if op2 < 0 ; in that case, we - set the result to 0 and return ARITH_DIV0. */ - mpz_set_si (result->value.integer, 0); - if (mpz_cmp_si (op2->value.integer, 0) < 0) - rc = ARITH_DIV0; - } - else if (mpz_cmp_si (op1->value.integer, -1) == 0) + switch (op1->ts.type) + { + case BT_INTEGER: { - /* (-1)**op2 == (-1)**(mod(op2,2)) */ - unsigned int odd = mpz_fdiv_ui (op2->value.integer, 2); - if (odd) - mpz_set_si (result->value.integer, -1); + int power; + + /* First, we simplify the cases of op1 == 1, 0 or -1. */ + if (mpz_cmp_si (op1->value.integer, 1) == 0) + { + /* 1**op2 == 1 */ + mpz_set_si (result->value.integer, 1); + } + else if (mpz_cmp_si (op1->value.integer, 0) == 0) + { + /* 0**op2 == 0, if op2 > 0 + 0**op2 overflow, if op2 < 0 ; in that case, we + set the result to 0 and return ARITH_DIV0. */ + mpz_set_si (result->value.integer, 0); + if (mpz_cmp_si (op2->value.integer, 0) < 0) + rc = ARITH_DIV0; + } + else if (mpz_cmp_si (op1->value.integer, -1) == 0) + { + /* (-1)**op2 == (-1)**(mod(op2,2)) */ + unsigned int odd = mpz_fdiv_ui (op2->value.integer, 2); + if (odd) + mpz_set_si (result->value.integer, -1); + else + mpz_set_si (result->value.integer, 1); + } + /* Then, we take care of op2 < 0. */ + else if (mpz_cmp_si (op2->value.integer, 0) < 0) + { + /* if op2 < 0, op1**op2 == 0 because abs(op1) > 1. */ + mpz_set_si (result->value.integer, 0); + } + else if (gfc_extract_int (op2, &power) != NULL) + { + /* If op2 doesn't fit in an int, the exponentiation will + overflow, because op2 > 0 and abs(op1) > 1. */ + mpz_t max; + int i; + i = gfc_validate_kind (BT_INTEGER, result->ts.kind, false); + + if (gfc_option.flag_range_check) + rc = ARITH_OVERFLOW; + + /* Still, we want to give the same value as the + processor. */ + mpz_init (max); + mpz_add_ui (max, gfc_integer_kinds[i].huge, 1); + mpz_mul_ui (max, max, 2); + mpz_powm (result->value.integer, op1->value.integer, + op2->value.integer, max); + mpz_clear (max); + } else - mpz_set_si (result->value.integer, 1); - } - /* Then, we take care of op2 < 0. */ - else if (mpz_cmp_si (op2->value.integer, 0) < 0) - { - /* if op2 < 0, op1**op2 == 0 because abs(op1) > 1. */ - mpz_set_si (result->value.integer, 0); + mpz_pow_ui (result->value.integer, op1->value.integer, + power); } - else if (gfc_extract_int (op2, &power) != NULL) + break; + + case BT_REAL: + mpfr_pow_z (result->value.real, op1->value.real, + op2->value.integer, GFC_RND_MODE); + break; + + case BT_COMPLEX: { - /* If op2 doesn't fit in an int, the exponentiation will - overflow, because op2 > 0 and abs(op1) > 1. */ - mpz_t max; - int i = gfc_validate_kind (BT_INTEGER, result->ts.kind, false); - - if (gfc_option.flag_range_check) - rc = ARITH_OVERFLOW; - - /* Still, we want to give the same value as the processor. */ - mpz_init (max); - mpz_add_ui (max, gfc_integer_kinds[i].huge, 1); - mpz_mul_ui (max, max, 2); - mpz_powm (result->value.integer, op1->value.integer, - op2->value.integer, max); - mpz_clear (max); + mpz_t apower; + + /* Compute op1**abs(op2) */ + mpz_init (apower); + mpz_abs (apower, op2->value.integer); + complex_pow (result, op1, apower); + mpz_clear (apower); + + /* If (op2 < 0), compute the inverse. */ + if (power_sign < 0) + complex_reciprocal (result); } - else - mpz_pow_ui (result->value.integer, op1->value.integer, power); - } - break; + break; - case BT_REAL: - mpfr_pow_z (result->value.real, op1->value.real, op2->value.integer, - GFC_RND_MODE); - break; + default: + break; + } + } + break; + + case BT_REAL: + + if (init_flag) + { + if (gfc_notify_std (GFC_STD_F2003,"Fortran 2003: Noninteger " + "exponent in an initialization " + "expression at %L", &op2->where) == FAILURE) + return ARITH_PROHIBIT; + } - case BT_COMPLEX: + if (mpfr_cmp_si (op1->value.real, 0) < 0) + { + gfc_error ("Raising a negative REAL at %L to " + "a REAL power is prohibited", &op1->where); + gfc_free (result); + return ARITH_PROHIBIT; + } + + mpfr_pow (result->value.real, op1->value.real, op2->value.real, + GFC_RND_MODE); + break; + + case BT_COMPLEX: + { + mpfr_t x, y, r, t; + + if (init_flag) { - mpz_t apower; + if (gfc_notify_std (GFC_STD_F2003,"Fortran 2003: Noninteger " + "exponent in an initialization " + "expression at %L", &op2->where) == FAILURE) + return ARITH_PROHIBIT; + } - /* Compute op1**abs(op2) */ - mpz_init (apower); - mpz_abs (apower, op2->value.integer); - complex_pow (result, op1, apower); - mpz_clear (apower); + gfc_set_model (op1->value.complex.r); - /* If (op2 < 0), compute the inverse. */ - if (power_sign < 0) - complex_reciprocal (result); + mpfr_init (r); + mpfr_hypot (r, op1->value.complex.r, op1->value.complex.i, + GFC_RND_MODE); + if (mpfr_cmp_si (r, 0) == 0) + { + mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE); + mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE); + mpfr_clear (r); break; } - - default: - break; - } + mpfr_log (r, r, GFC_RND_MODE); + + mpfr_init (t); + + mpfr_atan2 (t, op1->value.complex.i, op1->value.complex.r, + GFC_RND_MODE); + + mpfr_init (x); + mpfr_init (y); + + mpfr_mul (x, op2->value.complex.r, r, GFC_RND_MODE); + mpfr_mul (y, op2->value.complex.i, t, GFC_RND_MODE); + mpfr_sub (x, x, y, GFC_RND_MODE); + mpfr_exp (x, x, GFC_RND_MODE); + + mpfr_mul (y, op2->value.complex.r, t, GFC_RND_MODE); + mpfr_mul (t, op2->value.complex.i, r, GFC_RND_MODE); + mpfr_add (y, y, t, GFC_RND_MODE); + mpfr_cos (t, y, GFC_RND_MODE); + mpfr_sin (y, y, GFC_RND_MODE); + mpfr_mul (result->value.complex.r, x, t, GFC_RND_MODE); + mpfr_mul (result->value.complex.i, x, y, GFC_RND_MODE); + mpfr_clears (r, t, x, y, NULL); + } + break; + default: + gfc_internal_error ("arith_power(): unknown type"); } if (rc == ARITH_OK) @@ -1695,10 +1777,6 @@ eval_intrinsic (gfc_intrinsic_op op, gfc_internal_error ("eval_intrinsic(): Bad operator"); } - /* Try to combine the operators. */ - if (op == INTRINSIC_POWER && op2->ts.type != BT_INTEGER) - goto runtime; - if (op1->expr_type != EXPR_CONSTANT && (op1->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (op1) || !gfc_expanded_ac (op1))) @@ -1715,8 +1793,13 @@ eval_intrinsic (gfc_intrinsic_op op, else rc = reduce_binary (eval.f3, op1, op2, &result); + + /* Something went wrong. */ + if (op == INTRINSIC_POWER && rc == ARITH_PROHIBIT) + return NULL; + if (rc != ARITH_OK) - { /* Something went wrong. */ + { gfc_error (gfc_arith_error (rc), &op1->where); return NULL; } @@ -1908,7 +1991,7 @@ gfc_divide (gfc_expr *op1, gfc_expr *op2) gfc_expr * gfc_power (gfc_expr *op1, gfc_expr *op2) { - return eval_intrinsic_f3 (INTRINSIC_POWER, gfc_arith_power, op1, op2); + return eval_intrinsic_f3 (INTRINSIC_POWER, arith_power, op1, op2); } diff --git a/gcc/fortran/expr.c b/gcc/fortran/expr.c index 50444e4..8dec53f 100644 --- a/gcc/fortran/expr.c +++ b/gcc/fortran/expr.c @@ -1938,16 +1938,6 @@ check_intrinsic_op (gfc_expr *e, gfc_try (*check_function) (gfc_expr *)) if (!numeric_type (et0 (op1)) || !numeric_type (et0 (op2))) goto not_numeric; - if (e->value.op.op == INTRINSIC_POWER - && check_function == check_init_expr && et0 (op2) != BT_INTEGER) - { - if (gfc_notify_std (GFC_STD_F2003,"Fortran 2003: Noninteger " - "exponent in an initialization " - "expression at %L", &op2->where) - == FAILURE) - return FAILURE; - } - break; case INTRINSIC_CONCAT: @@ -2424,7 +2414,11 @@ gfc_reduce_init_expr (gfc_expr *expr) /* Match an initialization expression. We work by first matching an - expression, then reducing it to a constant. */ + expression, then reducing it to a constant. The reducing it to + constant part requires a global variable to flag the prohibition + of a non-integer exponent in -std=f95 mode. */ + +bool init_flag = false; match gfc_match_init_expr (gfc_expr **result) @@ -2435,18 +2429,25 @@ gfc_match_init_expr (gfc_expr **result) expr = NULL; + init_flag = true; + m = gfc_match_expr (&expr); if (m != MATCH_YES) - return m; + { + init_flag = false; + return m; + } t = gfc_reduce_init_expr (expr); if (t != SUCCESS) { gfc_free_expr (expr); + init_flag = false; return MATCH_ERROR; } *result = expr; + init_flag = false; return MATCH_YES; } diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h index 0bc5596..3a7f98a 100644 --- a/gcc/fortran/gfortran.h +++ b/gcc/fortran/gfortran.h @@ -199,7 +199,7 @@ gfc_intrinsic_op; /* Arithmetic results. */ typedef enum { ARITH_OK = 1, ARITH_OVERFLOW, ARITH_UNDERFLOW, ARITH_NAN, - ARITH_DIV0, ARITH_INCOMMENSURATE, ARITH_ASYMMETRIC + ARITH_DIV0, ARITH_INCOMMENSURATE, ARITH_ASYMMETRIC, ARITH_PROHIBIT } arith; diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog index 4588cab..6061a23 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,3 +1,8 @@ +2009-03-29 Steven G. Kargl <kargl@gcc.gnu.org> + + PR fortran/38823 + * gfortran.dg/power1.f90: New test. + 2009-03-29 Joseph Myers <joseph@codesourcery.com> PR c/456 diff --git a/gcc/testsuite/gfortran.dg/power1.f90 b/gcc/testsuite/gfortran.dg/power1.f90 new file mode 100644 index 0000000..50dbac2 --- /dev/null +++ b/gcc/testsuite/gfortran.dg/power1.f90 @@ -0,0 +1,58 @@ +! { dg-do run } +! Test fix for PR fortran/38823. +program power + + implicit none + + integer, parameter :: & + & s = kind(1.e0), & + & d = kind(1.d0), & + & e = max(selected_real_kind(precision(1.d0)+1), d) + + real(s), parameter :: ris = 2.e0_s**2 + real(d), parameter :: rid = 2.e0_d**2 + real(e), parameter :: rie = 2.e0_e**2 + complex(s), parameter :: cis = (2.e0_s,1.e0_s)**2 + complex(d), parameter :: cid = (2.e0_d,1.e0_d)**2 + complex(e), parameter :: cie = (2.e0_e,1.e0_e)**2 + + real(s), parameter :: rrs = 2.e0_s**2.e0 + real(d), parameter :: rrd = 2.e0_d**2.e0 + real(e), parameter :: rre = 2.e0_e**2.e0 + complex(s), parameter :: crs = (2.e0_s,1.e0_s)**2.e0 + complex(d), parameter :: crd = (2.e0_d,1.e0_d)**2.e0 + complex(e), parameter :: cre = (2.e0_e,1.e0_e)**2.e0 + + real(s), parameter :: rds = 2.e0_s**2.e0_d + real(d), parameter :: rdd = 2.e0_d**2.e0_d + real(e), parameter :: rde = 2.e0_e**2.e0_d + complex(s), parameter :: cds = (2.e0_s,1.e0_s)**2.e0_d + complex(d), parameter :: cdd = (2.e0_d,1.e0_d)**2.e0_d + complex(e), parameter :: cde = (2.e0_e,1.e0_e)**2.e0_d + + real(s), parameter :: eps_s = 1.e-5_s + real(d), parameter :: eps_d = 1.e-10_d + real(e), parameter :: eps_e = 1.e-10_e + + if (abs(ris - 4) > eps_s) call abort + if (abs(rid - 4) > eps_d) call abort + if (abs(rie - 4) > eps_e) call abort + if (abs(real(cis, s) - 3) > eps_s .or. abs(aimag(cis) - 4) > eps_s) call abort + if (abs(real(cid, d) - 3) > eps_d .or. abs(aimag(cid) - 4) > eps_d) call abort + if (abs(real(cie, e) - 3) > eps_e .or. abs(aimag(cie) - 4) > eps_e) call abort + + if (abs(rrs - 4) > eps_s) call abort + if (abs(rrd - 4) > eps_d) call abort + if (abs(rre - 4) > eps_e) call abort + if (abs(real(crs, s) - 3) > eps_s .or. abs(aimag(crs) - 4) > eps_s) call abort + if (abs(real(crd, d) - 3) > eps_d .or. abs(aimag(crd) - 4) > eps_d) call abort + if (abs(real(cre, e) - 3) > eps_e .or. abs(aimag(cre) - 4) > eps_e) call abort + + if (abs(rds - 4) > eps_s) call abort + if (abs(rdd - 4) > eps_d) call abort + if (abs(rde - 4) > eps_e) call abort + if (abs(real(cds, s) - 3) > eps_s .or. abs(aimag(cds) - 4) > eps_s) call abort + if (abs(real(cdd, d) - 3) > eps_d .or. abs(aimag(cdd) - 4) > eps_d) call abort + if (abs(real(cde, e) - 3) > eps_e .or. abs(aimag(cde) - 4) > eps_e) call abort + +end program power |