[gfortran] Don't error out for integral zeroth power of zero wasRe: [OT] __builtin_cpow((0,0),(0,0))

[Tobias Schlüter <tobias dot schlueter at physik dot uni-muenchen dot de>]

Joost VandeVondele wrote:
> not quite that silly, as you can see running this program with and without 
> the line commented.
> 
> integer :: i=0
> real :: a=0.0
> !write(6,*) 0.0**0
> write(6,*) a**i
> END

I had missed the fact that your exponent is an integer, not a real.  In this
case the result is indeed well-defined, as there's only one limit involved, as
the exponentiation operation is then defined as
 x**n = lim_{y->x} exp (n * ln (y))
(with only the branch cut ambiguity remaining) This doesn't define the case of
an integer base though, as continuity doesn't make much sense for functions
defined over the set of integers.  Symmetry between positive and negative
exponents together with the standard's definition of integral exponentiation
supports 0**0 == 1, but it is undefined in the standard).

(OTOH with x and y both continuously variable, for x == y == 0, x**y is not
uniquely defined, as
 x**y := lim_{w->x} lim_{u->y} exp (u * ln (w))
and
 x**y := lim_{u->y} lim_{w->x} exp (u * ln (w))
will yield different results in the pathological case we're discussing.  But
this doesn't matter for the patch, as we only constant-fold integral exponents
during compilation, and leave this case to the library)

I propose to either change arith.c to not error out at all in this special
case, and make 0**0 == 1, or to only give the error we're currently giving
("Indeterminate form 0 ** 0 at") in the case of an integral base.  A patch for
the former is below.

Remarks?

- Tobi

2005-03-09  Tobias Schl"uter  <tobias.schlueter@physik.uni-muenchen.de>

	* arith.c (gfc_arith_error): Remove case ARITH_0TO0.
	(gfc_arith_power): Consistently simplify to one for integral
	zero exponent.
	* gfortran.h (arith): Remove ARITH_0TO0.

Index: arith.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/fortran/arith.c,v
retrieving revision 1.23
diff -u -p -r1.23 arith.c
--- arith.c     27 Feb 2005 17:32:26 -0000      1.23
+++ arith.c     9 Mar 2005 19:29:44 -0000
@@ -152,9 +152,6 @@ gfc_arith_error (arith code)
     case ARITH_DIV0:
       p = "Division by zero";
       break;
-    case ARITH_0TO0:
-      p = "Indeterminate form 0 ** 0";
-      break;
     case ARITH_INCOMMENSURATE:
       p = "Array operands are incommensurate";
       break;
@@ -918,33 +915,23 @@ gfc_arith_power (gfc_expr * op1, gfc_exp
   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

   if (power == 0)
-    {                          /* Handle something to the zeroth power */
+    {
+      /* 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:
-         if (mpz_sgn (op1->value.integer) == 0)
-           rc = ARITH_0TO0;
-         else
-           mpz_set_ui (result->value.integer, 1);
+         mpz_set_ui (result->value.integer, 1);
          break;

        case BT_REAL:
-         if (mpfr_sgn (op1->value.real) == 0)
-           rc = ARITH_0TO0;
-         else
-           mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
+         mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
          break;

        case BT_COMPLEX:
-         if (mpfr_sgn (op1->value.complex.r) == 0
-             && mpfr_sgn (op1->value.complex.i) == 0)
-           rc = ARITH_0TO0;
-         else
-           {
-             mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
-             mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
-           }
-
+         mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
+         mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
          break;

        default:
Index: gfortran.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/fortran/gfortran.h,v
retrieving revision 1.59
diff -u -p -r1.59 gfortran.h
--- gfortran.h  28 Feb 2005 00:38:12 -0000      1.59
+++ gfortran.h  9 Mar 2005 19:39:18 -0000
@@ -181,7 +181,7 @@ extern mstring intrinsic_operators[];
 /* Arithmetic results.  */
 typedef enum
 { ARITH_OK = 1, ARITH_OVERFLOW, ARITH_UNDERFLOW, ARITH_NAN,
-  ARITH_DIV0, ARITH_0TO0, ARITH_INCOMMENSURATE, ARITH_ASYMMETRIC
+  ARITH_DIV0, ARITH_INCOMMENSURATE, ARITH_ASYMMETRIC
 }
 arith;