This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[PATCH, fortran, 4.5] Fix constant folding of noninteger exponents.
- From: Steve Kargl <sgk at troutmask dot apl dot washington dot edu>
- To: fortran at gcc dot gnu dot org, gcc-patches at gcc dot gnu dot org
- Date: Sun, 18 Jan 2009 16:31:40 -0800
- Subject: [PATCH, fortran, 4.5] Fix constant folding of noninteger exponents.
The attached patch has been regression tested on i386-*-freebsd
and x86_64-*-freebsd. There were no regressions.
This patch is not obvious. :-O
The original PR fortran/388232 was for a program of the form
program z
x = (-2.e0)**3.e0
end program z
where the Fortran 95 (page 97) and 2003 standards have
"Raising a negative-valued primary of type real to a real power is
prohibited."
The patch morphed into something that (partially) also addressed
PR fortran/38822. I'll get back to that in a moment.
Constant folding occurs before the checking of an initialization
expression for conformance, and AFAICT there is no way to tell
if gfortran is currently parsing an initialization expression.
The crux of the problem comes from page 94 of the F95 standard:
"An initialization expression is a constant expression in which
the exponentiation operation is permitted only with an integer
power, and ..."
The "in which ..." clause has been removed in Fortran 2003. gfortran
currently only enforces this checking with -std=f95.
kargl[11] cat > a.f90
program z
real :: x = 2.e0**3.e0
end program z
kargl[12] gfc43 -o z a.f90
kargl[13] gfc43 -o z -std=f95 a.f90
a.f90:2.25:
real :: x = 2.e0**3.e0
1
Error: Fortran 2003: Noninteger exponent in an initialization expression at (1)
So, gfortran appears to be doing what one would want. Unfortnately, this
isn't the case. Consider
kargl[19] cat a.f90
program z
real :: x = transfer(2.e0**3.e0, 1)
end program z
kargl[20] gfc43 -o -std=f95 z a.f90
f951: internal compiler error: in gfc_target_encode_expr, at
fortran/target-memory.c:235
There are clearly 2 problems here. First, there is an ICE. Second,
we do not get a warning/error for the non-integer exponent. The cause
of the ICE is that the gfortran frontend does not constant fold the
above initialization express. It builds a runtime expression, which
the middle-end would eventually fold if the 2.e0**3.e0 wasn't the
argument to TRANSFER. The presence of the runtime expression is
the origins of the ICE. arith.c(eval_intrinsics) contains
/* Try to combine the operators. */
if (op == INTRINSIC_POWER && op2->ts.type != BT_INTEGER)
goto runtime;
If we remove this restriction to allow the FE to constant fold, we
open up a Pandora's box because the F95 restrictions can only be
detected at a lower level in arith.c(gfc_arith_power). Note, I've
renamed gfc_arith_power to arith_pow. This function normally returns
an error code that is used to report errors. None of the current
error codes adequately describe the error. So, I've introduced
ARITH_PROHIBIT. I've used ARITH_PROHIBIT to detect both restrictions.
kargl[25] gfc4x -o z -std=f95 a.f90
a.f90:2.34:
real :: x = transfer(2.e0**3.e0, 1)
1
Error: Fortran 2003: Noninteger exponent in an initialization expression at (1)
a.f90:3.8:
y = (-2.0)**3.e0
1
Error: Raising a negative REAL at (1) to a REAL power is prohibited
The careful reader will note that expression for y is not flagged
as having an non-integer exponent. To accomplish, I've added a global
flag init_flag that is set to true when gfc_match_init_expr is entered
and cleared before returning.
Now back to PR fortran/38822, a very bare bone code fragment is
kargl[36] gfc4x -o z -std=f95 a.f90
kargl[37] cat a.f90
program z
real x(0*transfer(2.e0**3.e0, 1) + 2)
end program z
kargl[38] gfc4x -o z -std=f95 a.f90
Notice that the noninteger exponent restriction also applied to
specification statements. I haven't found a patch to fix this
yet.
Enjoy.
2009-01-17 Steven G. Kargl <kargl@gcc.gnu.org>
PR fortran/38823
* gfortran.dg/power1.f90: New test.
2009-01-17 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
--
Steve
Index: gcc/testsuite/gfortran.dg/power1.f90
===================================================================
--- gcc/testsuite/gfortran.dg/power1.f90 (revision 0)
+++ gcc/testsuite/gfortran.dg/power1.f90 (revision 0)
@@ -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
Index: gcc/fortran/gfortran.h
===================================================================
--- gcc/fortran/gfortran.h (revision 143457)
+++ gcc/fortran/gfortran.h (working copy)
@@ -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;
Index: gcc/fortran/expr.c
===================================================================
--- gcc/fortran/expr.c (revision 143457)
+++ gcc/fortran/expr.c (working copy)
@@ -1938,16 +1938,6 @@ check_intrinsic_op (gfc_expr *e, gfc_try
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;
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;
}
Index: gcc/fortran/arith.c
===================================================================
--- gcc/fortran/arith.c (revision 143457)
+++ gcc/fortran/arith.c (working copy)
@@ -932,131 +932,213 @@ complex_pow (gfc_expr *result, gfc_expr
}
-/* 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;
+ }
+
+ 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:
+ 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);
-
- /* If (op2 < 0), compute the inverse. */
- if (power_sign < 0)
- complex_reciprocal (result);
+ gfc_set_model (op1->value.complex.r);
+ 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;
}
+ mpfr_log (r, r, GFC_RND_MODE);
- default:
- break;
- }
+ 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);
}