[PATCH]: Use mpc_pow in the fortran frontend

["Kaveh R. GHAZI" <ghazi at caip dot rutgers dot edu>]

This patch makes use of mpc_pow in the fortran frontend.  Note mpc_pow
is only available in the svn repo for MPC.  I've guarded the uses of
mpc_pow with the appropriate version checks.

The call to mpc_pow replaces two separate areas, a general pow, and
the integer exponent specialization.  The existing fortran integer
exponent code had a problem immortalized in a testcase, and as
discussed here:
http://gcc.gnu.org/ml/fortran/2009-06/msg00267.html
I'm commenting out some lines in one testcase that were incorrect.

I'll reinsert the corrected lines once we've made MPC a
hard-requirement, and assuming this patch is accepted I'll put a note
in PR40302 so I don't forget.

Bootstrapped on x86_64-unknown-linux-gnu with MPC svn (which has
mpc_pow), with mpc-0.6 (no pow), and without MPC.  All three scenarios
had no regressions.

Okay for mainline?

		Thanks,
		--Kaveh


2009-06-26  Kaveh R. Ghazi  <ghazi@caip.rutgers.edu>

gcc/fortran:

	* gfortran.h: Define HAVE_mpc_pow.
	* arith.c (complex_reciprocal, complex_pow): If HAVE_mpc_pow,
	don't define these functions.
	(arith_power): If HAVE_mpc_pow, use mpc_pow.

gcc/testsuite:
	* gfortran.dg/integer_exponentiation_4.f90: Temporarily
	comment out some values and add some cases.

diff -rup orig/egcc-SVN20090625/gcc/fortran/gfortran.h egcc-SVN20090625/gcc/fortran/gfortran.h
--- orig/egcc-SVN20090625/gcc/fortran/gfortran.h	2009-06-25 02:01:27.000000000 +0200
+++ egcc-SVN20090625/gcc/fortran/gfortran.h	2009-06-25 22:54:40.000000000 +0200
@@ -1558,6 +1558,9 @@ gfc_intrinsic_sym;
 #include <mpfr.h>
 #ifdef HAVE_mpc
 #include <mpc.h>
+# if MPC_VERSION >= MPC_VERSION_NUM(0,6,1)
+#  define HAVE_mpc_pow
+# endif
 #else
 #define mpc_realref(X) ((X).r)
 #define mpc_imagref(X) ((X).i)
diff -rup orig/egcc-SVN20090625/gcc/fortran/arith.c egcc-SVN20090625/gcc/fortran/arith.c
--- orig/egcc-SVN20090625/gcc/fortran/arith.c	2009-06-19 17:35:32.000000000 +0200
+++ egcc-SVN20090625/gcc/fortran/arith.c	2009-06-25 22:54:40.000000000 +0200
@@ -896,6 +896,7 @@ gfc_arith_divide (gfc_expr *op1, gfc_exp

 /* Compute the reciprocal of a complex number (guaranteed nonzero).  */

+#if ! defined(HAVE_mpc_pow)
 static void
 complex_reciprocal (gfc_expr *op)
 {
@@ -922,6 +923,7 @@ complex_reciprocal (gfc_expr *op)
   }
 #endif
 }
+#endif /* ! HAVE_mpc_pow */


 /* Raise a complex number to positive power (power > 0).
@@ -932,6 +934,7 @@ complex_reciprocal (gfc_expr *op)
    "Seminumerical Algorithms", Vol. 2, "The Art of Computer Programming",
    3rd Edition, 1998.  */

+#if ! defined(HAVE_mpc_pow)
 static void
 complex_pow (gfc_expr *result, gfc_expr *base, mpz_t power)
 {
@@ -988,6 +991,7 @@ complex_pow (gfc_expr *result, gfc_expr

   mpfr_clears (x_r, x_i, tmp, re, im, NULL);
 }
+#endif /* ! HAVE_mpc_pow */


 /* Raise a number to a power.  */
@@ -1107,6 +1111,15 @@ arith_power (gfc_expr *op1, gfc_expr *op

 	    case BT_COMPLEX:
 	      {
+#ifdef HAVE_mpc_pow
+		mpc_t apower;
+		gfc_set_model (mpc_realref (op1->value.complex));
+		mpc_init2 (apower, mpfr_get_default_prec());
+		mpc_set_z (apower, op2->value.integer, GFC_MPC_RND_MODE);
+		mpc_pow(result->value.complex, op1->value.complex, apower,
+			GFC_MPC_RND_MODE);
+		mpc_clear (apower);
+#else
 		mpz_t apower;

 		/* Compute op1**abs(op2)  */
@@ -1118,6 +1131,7 @@ arith_power (gfc_expr *op1, gfc_expr *op
 		/* If (op2 < 0), compute the inverse.  */
 		if (power_sign < 0)
 		  complex_reciprocal (result);
+#endif /* HAVE_mpc_pow */
 	      }
 	      break;

@@ -1159,6 +1173,10 @@ arith_power (gfc_expr *op1, gfc_expr *op
 	      return ARITH_PROHIBIT;
 	  }

+#ifdef HAVE_mpc_pow
+	mpc_pow (result->value.complex, op1->value.complex,
+		 op2->value.complex, GFC_MPC_RND_MODE);
+#else
 	{
 	mpfr_t x, y, r, t;

@@ -1211,6 +1229,7 @@ arith_power (gfc_expr *op1, gfc_expr *op
 	mpfr_mul (mpc_imagref (result->value.complex), x, y, GFC_RND_MODE);
 	mpfr_clears (r, t, x, y, NULL);
 	}
+#endif /* HAVE_mpc_pow */
       }
       break;
     default:
diff -rup orig/egcc-SVN20090625/gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90 egcc-SVN20090625/gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90
--- orig/egcc-SVN20090625/gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90	2008-03-14 00:37:50.000000000 +0100
+++ egcc-SVN20090625/gcc/testsuite/gfortran.dg/integer_exponentiation_4.f90	2009-06-27 02:16:40.000000000 +0200
@@ -38,7 +38,10 @@ program test
   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" }
+! Put these lines back in (and "no-" -> "dg-") prior to gcc-4.5.
+!  print *, (2.0,-4.3)**huge(0) ! { no-error "Arithmetic overflow" }
+!  print *, (2.0,-4.3)**huge(0_8) ! { no-error "Arithmetic overflow" }
+!  print *, (2.0,-4.3)**(-huge(0))
+!  print *, (2.0,-4.3)**(-huge(0_8))

 end program test