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

gcc/fortran/ChangeLog

+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

gcc/fortran/arith.c

@@ -932,19 +932,22 @@ 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);
+
+  switch (op2->ts.type)
+    {
+    case BT_INTEGER:
       power_sign = mpz_sgn (op2->value.integer);
 
       if (power_sign == 0)
@@ -968,7 +971,7 @@ gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 	      break;
 
 	    default:
-	  gfc_internal_error ("gfc_arith_power(): Bad base");
+	      gfc_internal_error ("arith_power(): Bad base");
 	    }
 	}
       else
@@ -1014,12 +1017,14 @@ gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 		    /* 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);
+		    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.  */
+		    /* 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);
@@ -1028,13 +1033,14 @@ gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 		    mpz_clear (max);
 		  }
 		else
-	      mpz_pow_ui (result->value.integer, op1->value.integer, power);
+		  mpz_pow_ui (result->value.integer, op1->value.integer,
+			      power);
 	      }
 	      break;
 
 	    case BT_REAL:
-	  mpfr_pow_z (result->value.real, op1->value.real, op2->value.integer,
-		      GFC_RND_MODE);
+	      mpfr_pow_z (result->value.real, op1->value.real,
+			  op2->value.integer, GFC_RND_MODE);
 	      break;
 
 	    case BT_COMPLEX:
@@ -1050,13 +1056,89 @@ gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 		/* If (op2 < 0), compute the inverse.  */
 		if (power_sign < 0)
 		  complex_reciprocal (result);
+	      }
+	      break;
 
+	    default:
 	      break;
 	    }
+	}
+      break;
 
-	default:
+    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:
+      {
+	mpfr_t x, y, r, t;
+
+	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;
 	  }
+
+	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);
+
+	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);
 }
 
 

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)
+    {
+      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;
 }

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;
 

gcc/testsuite/ChangeLog

+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

gcc/testsuite/gfortran.dg/power1.f90

+! { 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