re PR fortran/38823 (Diagnose and treat (-2.0)**2.0 properly)

2009-03-29  Steven G. Kargl  <kargl@gcc.gnu.org>

        PR fortran/38823
        * gfortran.dg/power1.f90: New test.

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

From-SVN: r145261

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 3156df6..8bdf010 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
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