Re: [gfortran,patch] Fix constant exponentiation with integer exponent

["François-Xavier Coudert" <fxcoudert at gmail dot com>]

> The patch comes with somewhat long testcases. Regtested on i686-
> linux, OK for mainline? Do you think it's OK for 4.2? (I'm not sure,
> it's rather invasive... and the cases where the compiler actually got
> it wrong were... rather rarely used)

After some more thinking (and some testing of the trans-const.c patch
for very large integers with -fno-range-check), here's an updated
patch. The changes are small, they only make sure that we don't reject
i**j when all of the following conditions are fulfilled:
 (a) i and j are integers
 (b) i is larger than huge(0_4)
 (c) i**j overflows (that is, j>0 and i>1)
 (d) -fno-range-check is specified.
The integer_exponentiation_5.F90 testcase, which uses
-fno-range-check, was updated to test this code path.

Regtested once more on i686-linux (it needs the patch for PR31262,
which can be found here:
http://gcc.gnu.org/ml/gcc-patches/2007-03/msg01332.html ; Jerry is
currently reviewing it, so it should be commited soon). OK for
mainline? Still same question as quoted above for 4.2.

Thanks,
FX

Attachment: pr30834_4.ChangeLog
Description: Binary data

Index: gcc/fortran/arith.c
===================================================================
--- gcc/fortran/arith.c	(revision 123071)
+++ gcc/fortran/arith.c	(working copy)
@@ -872,42 +872,69 @@
 }
 
 
-/* Raise a complex number to positive power.  */
+/* Raise a complex number to positive power (power > 0).
+   This function will modify the content of power.
 
+   Use Binary Method, which is not an optimal but a simple and reasonable
+   arithmetic. See section 4.6.3, "Evaluation of Powers" of Donald E. Knuth,
+   "Seminumerical Algorithms", Vol. 2, "The Art of Computer Programming",
+   3rd Edition, 1998.  */
+
 static void
-complex_pow_ui (gfc_expr *base, int power, gfc_expr *result)
+complex_pow (gfc_expr *result, gfc_expr *base, mpz_t power)
 {
-  mpfr_t re, im, a;
+  mpfr_t x_r, x_i, tmp, re, im;
 
   gfc_set_model (base->value.complex.r);
+  mpfr_init (x_r);
+  mpfr_init (x_i);
+  mpfr_init (tmp);
   mpfr_init (re);
   mpfr_init (im);
-  mpfr_init (a);
 
+  /* res = 1 */
   mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 
-  for (; power > 0; power--)
+  /* x = base */
+  mpfr_set (x_r, base->value.complex.r, GFC_RND_MODE);
+  mpfr_set (x_i, base->value.complex.i, GFC_RND_MODE);
+
+/* Macro for complex multiplication. We have to take care that
+   res_r/res_i and a_r/a_i can (and will) be the same variable.  */
+#define CMULT(res_r,res_i,a_r,a_i,b_r,b_i) \
+    mpfr_mul (re, a_r, b_r, GFC_RND_MODE), \
+    mpfr_mul (tmp, a_i, b_i, GFC_RND_MODE), \
+    mpfr_sub (re, re, tmp, GFC_RND_MODE), \
+    \
+    mpfr_mul (im, a_r, b_i, GFC_RND_MODE), \
+    mpfr_mul (tmp, a_i, b_r, GFC_RND_MODE), \
+    mpfr_add (res_i, im, tmp, GFC_RND_MODE), \
+    mpfr_set (res_r, re, GFC_RND_MODE)
+  
+#define res_r result->value.complex.r
+#define res_i result->value.complex.i
+
+  /* for (; power > 0; x *= x) */
+  for (; mpz_cmp_si (power, 0) > 0; CMULT(x_r,x_i,x_r,x_i,x_r,x_i))
     {
-      mpfr_mul (re, base->value.complex.r, result->value.complex.r,
-		GFC_RND_MODE);
-      mpfr_mul (a, base->value.complex.i, result->value.complex.i,
-		GFC_RND_MODE);
-      mpfr_sub (re, re, a, GFC_RND_MODE);
+      /* if (power & 1) res = res * x; */
+      if (mpz_congruent_ui_p (power, 1, 2))
+	CMULT(res_r,res_i,res_r,res_i,x_r,x_i);
 
-      mpfr_mul (im, base->value.complex.r, result->value.complex.i,
-		GFC_RND_MODE);
-      mpfr_mul (a, base->value.complex.i, result->value.complex.r,
-		GFC_RND_MODE);
-      mpfr_add (im, im, a, GFC_RND_MODE);
-
-      mpfr_set (result->value.complex.r, re, GFC_RND_MODE);
-      mpfr_set (result->value.complex.i, im, GFC_RND_MODE);
+      /* power /= 2; */
+      mpz_fdiv_q_ui (power, power, 2);
     }
 
+#undef res_r
+#undef res_i
+#undef CMULT
+
+  mpfr_clear (x_r);
+  mpfr_clear (x_i);
+  mpfr_clear (tmp);
   mpfr_clear (re);
   mpfr_clear (im);
-  mpfr_clear (a);
 }
 
 
@@ -916,20 +943,17 @@
 static arith
 gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
-  int power, apower;
+  int power_sign;
   gfc_expr *result;
-  mpz_t unity_z;
-  mpfr_t unity_f;
   arith rc;
 
+  gcc_assert (op2->expr_type == EXPR_CONSTANT && op2->ts.type == BT_INTEGER);
+
   rc = ARITH_OK;
-
-  if (gfc_extract_int (op2, &power) != NULL)
-    gfc_internal_error ("gfc_arith_power(): Bad exponent");
-
   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+  power_sign = mpz_sgn (op2->value.integer);
 
-  if (power == 0)
+  if (power_sign == 0)
     {
       /* Handle something to the zeroth power.  Since we're dealing
 	 with integral exponents, there is no ambiguity in the
@@ -955,45 +979,87 @@
     }
   else
     {
-      apower = power;
-      if (power < 0)
-	apower = -power;
-
       switch (op1->ts.type)
 	{
 	case BT_INTEGER:
-	  mpz_pow_ui (result->value.integer, op1->value.integer, apower);
+	  {
+	    int power;
 
-	  if (power < 0)
-	    {
-	      mpz_init_set_ui (unity_z, 1);
-	      mpz_tdiv_q (result->value.integer, unity_z,
-			  result->value.integer);
-	      mpz_clear (unity_z);
-	    }
+	    /* 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 = gfc_validate_kind (BT_INTEGER, result->ts.kind, false);
+
+		if (gfc_option.flag_range_check)
+		  rc = ARITH_OVERFLOW;
+
+		/* Sill, 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_pow_ui (result->value.integer, op1->value.integer, power);
+	  }
 	  break;
 
 	case BT_REAL:
-	  mpfr_pow_ui (result->value.real, op1->value.real, apower,
-		       GFC_RND_MODE);
-
-	  if (power < 0)
-	    {
-	      gfc_set_model (op1->value.real);
-	      mpfr_init (unity_f);
-	      mpfr_set_ui (unity_f, 1, GFC_RND_MODE);
-	      mpfr_div (result->value.real, unity_f, result->value.real,
-			GFC_RND_MODE);
-	      mpfr_clear (unity_f);
-	    }
+	  mpfr_pow_z (result->value.real, op1->value.real, op2->value.integer,
+		      GFC_RND_MODE);
 	  break;
 
 	case BT_COMPLEX:
-	  complex_pow_ui (op1, apower, result);
-	  if (power < 0)
-	    complex_reciprocal (result);
-	  break;
+	  {
+	    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);
+
+	    break;
+	  }
+
 	default:
 	  break;
 	}
Index: gcc/testsuite/gfortran.dg/integer_exponentiation_3.F90
===================================================================
--- gcc/testsuite/gfortran.dg/integer_exponentiation_3.F90	(revision 0)
+++ gcc/testsuite/gfortran.dg/integer_exponentiation_3.F90	(revision 0)
@@ -0,0 +1,201 @@
+! { dg-do run }
+! { dg-options "" }
+module mod_check
+  implicit none
+
+  interface check
+    module procedure check_i8
+    module procedure check_i4
+    module procedure check_r8
+    module procedure check_r4
+    module procedure check_c8
+    module procedure check_c4
+  end interface check
+
+  interface acheck
+    module procedure acheck_c8
+    module procedure acheck_c4
+  end interface acheck
+
+contains
+
+  subroutine check_i8 (a, b)
+    integer(kind=8), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_i8
+
+  subroutine check_i4 (a, b)
+    integer(kind=4), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_i4
+
+  subroutine check_r8 (a, b)
+    real(kind=8), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_r8
+
+  subroutine check_r4 (a, b)
+    real(kind=4), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_r4
+
+  subroutine check_c8 (a, b)
+    complex(kind=8), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_c8
+
+  subroutine check_c4 (a, b)
+    complex(kind=4), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_c4
+
+  subroutine acheck_c8 (a, b)
+    complex(kind=8), intent(in) :: a, b
+    if (abs(a-b) > 1.d-9 * min(abs(a),abs(b))) call abort()
+  end subroutine acheck_c8
+
+  subroutine acheck_c4 (a, b)
+    complex(kind=4), intent(in) :: a, b
+    if (abs(a-b) > 1.e-5 * min(abs(a),abs(b))) call abort()
+  end subroutine acheck_c4
+
+end module mod_check
+
+program test
+  use mod_check
+  implicit none
+
+  integer(kind=4) :: i4
+  integer(kind=8) :: i8
+  real(kind=4) :: r4
+  real(kind=8) :: r8
+  complex(kind=4) :: c4
+  complex(kind=8) :: c8
+
+#define TEST(base,exp,var) var = base; call check((var)**(exp),(base)**(exp))
+#define ATEST(base,exp,var) var = base; call acheck((var)**(exp),(base)**(exp))
+
+!!!!! INTEGER BASE !!!!!
+  TEST(0,0,i4)
+  TEST(0_8,0_8,i8)
+  TEST(1,0,i4)
+  TEST(1_8,0_8,i8)
+  TEST(-1,0,i4)
+  TEST(-1_8,0_8,i8)
+  TEST(huge(0),0,i4)
+  TEST(huge(0_8),0_8,i8)
+  TEST(-huge(0)-1,0,i4)
+  TEST(-huge(0_8)-1_8,0_8,i8)
+
+  TEST(1,1,i4)
+  TEST(1_8,1_8,i8)
+  TEST(1,2,i4)
+  TEST(1_8,2_8,i8)
+  TEST(1,-1,i4)
+  TEST(1_8,-1_8,i8)
+  TEST(1,-2,i4)
+  TEST(1_8,-2_8,i8)
+  TEST(1,huge(0),i4)
+  TEST(1_8,huge(0_8),i8)
+  TEST(1,-huge(0)-1,i4)
+  TEST(1_8,-huge(0_8)-1_8,i8)
+
+  TEST(-1,1,i4)
+  TEST(-1_8,1_8,i8)
+  TEST(-1,2,i4)
+  TEST(-1_8,2_8,i8)
+  TEST(-1,-1,i4)
+  TEST(-1_8,-1_8,i8)
+  TEST(-1,-2,i4)
+  TEST(-1_8,-2_8,i8)
+  TEST(-1,huge(0),i4)
+  TEST(-1_8,huge(0_8),i8)
+  TEST(-1,-huge(0)-1,i4)
+  TEST(-1_8,-huge(0_8)-1_8,i8)
+
+  TEST(2,9,i4)
+  TEST(2_8,9_8,i8)
+  TEST(-2,9,i4)
+  TEST(-2_8,9_8,i8)
+  TEST(2,-9,i4)
+  TEST(2_8,-9_8,i8)
+  TEST(-2,-9,i4)
+  TEST(-2_8,-9_8,i8)
+
+!!!!! REAL BASE !!!!!
+  TEST(0.0,0,r4)
+  TEST(0.0,1,r4)
+  TEST(0.0,huge(0),r4)
+  TEST(0.0,0_8,r4)
+  TEST(0.0,1_8,r4)
+  TEST(0.0,huge(0_8),r4)
+
+  TEST(1.0,0,r4)
+  TEST(1.0,1,r4)
+  TEST(1.0,-1,r4)
+  TEST(1.0,huge(0),r4)
+  TEST(1.0,-huge(0)-1,r4)
+  TEST(1.0,0_8,r4)
+  TEST(1.0,1_8,r4)
+  TEST(1.0,-1_8,r4)
+  TEST(1.0,huge(0_8),r4)
+  TEST(1.0,-huge(0_8)-1_8,r4)
+
+  TEST(-1.0,0,r4)
+  TEST(-1.0,1,r4)
+  TEST(-1.0,-1,r4)
+  TEST(-1.0,huge(0),r4)
+  TEST(-1.0,-huge(0)-1,r4)
+  TEST(-1.0,0_8,r4)
+  TEST(-1.0,1_8,r4)
+  TEST(-1.0,-1_8,r4)
+  TEST(-1.0,huge(0_8),r4)
+  TEST(-1.0,-huge(0_8)-1_8,r4)
+
+  TEST(2.0,0,r4)
+  TEST(2.0,1,r4)
+  TEST(2.0,-1,r4)
+  TEST(2.0,3,r4)
+  TEST(2.0,-3,r4)
+  TEST(2.0,0_8,r4)
+  TEST(2.0,1_8,r4)
+  TEST(2.0,-1_8,r4)
+  TEST(2.0,3_8,r4)
+  TEST(2.0,-3_8,r4)
+
+  TEST(nearest(1.0,-1.0),0,r4)
+  TEST(nearest(1.0,-1.0),huge(0),r4) ! { dg-warning "Arithmetic underflow" }
+  TEST(nearest(1.0,-1.0),0_8,r4)
+  TEST(nearest(1.0_8,-1.0),huge(0_8),r8) ! { dg-warning "Arithmetic underflow" }
+
+  TEST(nearest(1.0,-1.0),107,r4)
+  TEST(nearest(1.0,1.0),107,r4)
+
+!!!!! COMPLEX BASE !!!!!
+  TEST((1.0,0.2),0,c4)
+  TEST((1.0,0.2),1,c4)
+  TEST((1.0,0.2),2,c4)
+  TEST((1.0,0.2),9,c4)
+  ATEST((1.0,0.2),-1,c4)
+  ATEST((1.0,0.2),-2,c4)
+  ATEST((1.0,0.2),-9,c4)
+
+  TEST((0.0,0.2),0,c4)
+  TEST((0.0,0.2),1,c4)
+  TEST((0.0,0.2),2,c4)
+  TEST((0.0,0.2),9,c4)
+  ATEST((0.0,0.2),-1,c4)
+  ATEST((0.0,0.2),-2,c4)
+  ATEST((0.0,0.2),-9,c4)
+
+  TEST((1.0,0.),0,c4)
+  TEST((1.0,0.),1,c4)
+  TEST((1.0,0.),2,c4)
+  TEST((1.0,0.),9,c4)
+  ATEST((1.0,0.),-1,c4)
+  ATEST((1.0,0.),-2,c4)
+  ATEST((1.0,0.),-9,c4)
+
+end program test
+
+! { dg-final { cleanup-modules "mod_check" } }
Index: gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90
===================================================================
--- gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90	(revision 0)
+++ gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90	(revision 0)
@@ -0,0 +1,44 @@
+! { dg-do compile }
+! { dg-options "" }
+program test
+  implicit none
+
+!!!!!! INTEGER BASE !!!!!!
+  print *, 0**0
+  print *, 0**1
+  print *, 0**(-1) ! { dg-error "Division by zero" }
+  print *, 0**(huge(0))
+  print *, 0**(-huge(0)-1) ! { dg-error "Division by zero" }
+  print *, 0**(2_8**32)
+  print *, 0**(-(2_8**32)) ! { dg-error "Division by zero" }
+
+  print *, 1**huge(0)
+  print *, 1**(-huge(0)-1)
+  print *, 1**huge(0_8)
+  print *, 1**(-huge(0_8)-1_8)
+  print *, (-1)**huge(0)
+  print *, (-1)**(-huge(0)-1)
+  print *, (-1)**huge(0_8)
+  print *, (-1)**(-huge(0_8)-1_8)
+
+  print *, 2**huge(0) ! { dg-error "Arithmetic overflow" }
+  print *, 2**huge(0_8) ! { dg-error "Arithmetic overflow" }
+  print *, (-2)**huge(0) ! { dg-error "Arithmetic overflow" }
+  print *, (-2)**huge(0_8) ! { dg-error "Arithmetic overflow" }
+
+  print *, 2**(-huge(0)-1)
+  print *, 2**(-huge(0_8)-1_8)
+  print *, (-2)**(-huge(0)-1)
+  print *, (-2)**(-huge(0_8)-1_8)
+
+!!!!!! REAL BASE !!!!!!
+  print *, 0.0**(-1) ! { dg-error "Arithmetic overflow" }
+  print *, 0.0**(-huge(0)-1) ! { dg-error "Arithmetic overflow" }
+  print *, 2.0**huge(0) ! { dg-error "Arithmetic overflow" }
+  print *, nearest(1.0,-1.0)**(-huge(0)) ! { dg-error "Arithmetic overflow" }
+
+!!!!!! COMPLEX BASE !!!!!!
+  print *, (2.0,-4.3)**huge(0) ! { dg-error "Arithmetic NaN" }
+  print *, (2.0,-4.3)**(-huge(0)) ! { dg-error "Arithmetic NaN" }
+
+end program test
Index: gcc/testsuite/gfortran.dg/integer_exponentiation_5.F90
===================================================================
--- gcc/testsuite/gfortran.dg/integer_exponentiation_5.F90	(revision 0)
+++ gcc/testsuite/gfortran.dg/integer_exponentiation_5.F90	(revision 0)
@@ -0,0 +1,78 @@
+! { dg-do run }
+! { dg-options "-fno-range-check" }
+module mod_check
+  implicit none
+
+  interface check
+    module procedure check_i8
+    module procedure check_i4
+    module procedure check_r8
+    module procedure check_r4
+    module procedure check_c8
+    module procedure check_c4
+  end interface check
+
+contains
+
+  subroutine check_i8 (a, b)
+    integer(kind=8), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_i8
+
+  subroutine check_i4 (a, b)
+    integer(kind=4), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_i4
+
+  subroutine check_r8 (a, b)
+    real(kind=8), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_r8
+
+  subroutine check_r4 (a, b)
+    real(kind=4), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_r4
+
+  subroutine check_c8 (a, b)
+    complex(kind=8), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_c8
+
+  subroutine check_c4 (a, b)
+    complex(kind=4), intent(in) :: a, b
+    if (a /= b) call abort()
+  end subroutine check_c4
+
+end module mod_check
+
+program test
+  use mod_check
+  implicit none
+
+  integer(kind=4) :: i4
+  integer(kind=8) :: i8
+  real(kind=4) :: r4
+  real(kind=8) :: r8
+  complex(kind=4) :: c4
+  complex(kind=8) :: c8
+
+#define TEST(base,exp,var) var = base; call check((var)**(exp),(base)**(exp))
+
+!!!!! INTEGER BASE !!!!!
+  TEST(3,23,i4)
+  TEST(-3,23,i4)
+  TEST(3_8,43_8,i8)
+  TEST(-3_8,43_8,i8)
+
+  TEST(17_8,int(huge(0),kind=8)+1,i8)
+
+!!!!! REAL BASE !!!!!
+  TEST(0.0,-1,r4)
+  TEST(0.0,-huge(0)-1,r4)
+  TEST(2.0,huge(0),r4)
+  TEST(nearest(1.0,-1.0),-huge(0),r4)
+
+end program test
+
+! { dg-final { cleanup-modules "mod_check" } }