re PR fortran/33197 (Fortran 2008: math functions)

	PR fortran/33197

	* intrinsic.c (add_functions): Use ERFC_SCALED simplification.
	* intrinsic.h (gfc_simplify_erfc_scaled): New prototype.
	* simplify.c (fullprec_erfc_scaled, asympt_erfc_scaled,
	gfc_simplify_erfc_scaled): New functions.

	* gfortran.dg/erf_2.F90: New test.
	* gfortran.dg/erfc_scaled_2.f90: New test.

From-SVN: r147621

gcc/fortran/ChangeLog

 2009-05-16  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
 
+	PR fortran/33197
+	* intrinsic.c (add_functions): Use ERFC_SCALED simplification.
+	* intrinsic.h (gfc_simplify_erfc_scaled): New prototype.
+	* simplify.c (fullprec_erfc_scaled, asympt_erfc_scaled,
+	gfc_simplify_erfc_scaled): New functions.
+
+2009-05-16  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
+
 	PR fortran/31243
 	* resolve.c (resolve_substring): Don't allow too large substring
 	indexes.

gcc/fortran/intrinsic.c

@@ -1431,8 +1431,9 @@ add_functions (void)
   make_generic ("erfc", GFC_ISYM_ERFC, GFC_STD_F2008);
 
   add_sym_1 ("erfc_scaled", GFC_ISYM_ERFC_SCALED, CLASS_ELEMENTAL, ACTUAL_NO,
-	     BT_REAL, dr, GFC_STD_F2008, gfc_check_fn_r, NULL,
-	     gfc_resolve_g77_math1, x, BT_REAL, dr, REQUIRED);
+	     BT_REAL, dr, GFC_STD_F2008, gfc_check_fn_r,
+	     gfc_simplify_erfc_scaled, gfc_resolve_g77_math1, x, BT_REAL,
+	     dr, REQUIRED);
 
   make_generic ("erfc_scaled", GFC_ISYM_ERFC_SCALED, GFC_STD_F2008);
 

gcc/fortran/intrinsic.h

@@ -232,6 +232,7 @@ gfc_expr *gfc_simplify_dprod (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_epsilon (gfc_expr *);
 gfc_expr *gfc_simplify_erf (gfc_expr *);
 gfc_expr *gfc_simplify_erfc (gfc_expr *);
+gfc_expr *gfc_simplify_erfc_scaled (gfc_expr *);
 gfc_expr *gfc_simplify_exp (gfc_expr *);
 gfc_expr *gfc_simplify_exponent (gfc_expr *);
 gfc_expr *gfc_simplify_float (gfc_expr *);

gcc/fortran/simplify.c

@@ -1213,6 +1213,143 @@ gfc_simplify_erfc (gfc_expr *x)
 }
 
 
+/* Helper functions to simplify ERFC_SCALED(x) = ERFC(x) * EXP(X**2).  */
+
+#define MAX_ITER 200
+#define ARG_LIMIT 12
+
+/* Calculate ERFC_SCALED directly by its definition:
+
+     ERFC_SCALED(x) = ERFC(x) * EXP(X**2)
+
+   using a large precision for intermediate results.  This is used for all
+   but large values of the argument.  */
+static void
+fullprec_erfc_scaled (mpfr_t res, mpfr_t arg)
+{
+  mp_prec_t prec;
+  mpfr_t a, b;
+
+  prec = mpfr_get_default_prec ();
+  mpfr_set_default_prec (10 * prec);
+
+  mpfr_init (a);
+  mpfr_init (b);
+
+  mpfr_set (a, arg, GFC_RND_MODE);
+  mpfr_sqr (b, a, GFC_RND_MODE);
+  mpfr_exp (b, b, GFC_RND_MODE);
+  mpfr_erfc (a, a, GFC_RND_MODE);
+  mpfr_mul (a, a, b, GFC_RND_MODE);
+
+  mpfr_set (res, a, GFC_RND_MODE);
+  mpfr_set_default_prec (prec);
+
+  mpfr_clear (a);
+  mpfr_clear (b);
+}
+
+/* Calculate ERFC_SCALED using a power series expansion in 1/arg:
+
+    ERFC_SCALED(x) = 1 / (x * sqrt(pi))
+                     * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
+                                          / (2 * x**2)**n)
+
+  This is used for large values of the argument.  Intermediate calculations
+  are performed with twice the precision.  We don't do a fixed number of
+  iterations of the sum, but stop when it has converged to the required
+  precision.  */
+static void
+asympt_erfc_scaled (mpfr_t res, mpfr_t arg)
+{
+  mpfr_t sum, x, u, v, w, oldsum, sumtrunc;
+  mpz_t num;
+  mp_prec_t prec;
+  unsigned i;
+
+  prec = mpfr_get_default_prec ();
+  mpfr_set_default_prec (2 * prec);
+
+  mpfr_init (sum);
+  mpfr_init (x);
+  mpfr_init (u);
+  mpfr_init (v);
+  mpfr_init (w);
+  mpz_init (num);
+
+  mpfr_init (oldsum);
+  mpfr_init (sumtrunc);
+  mpfr_set_prec (oldsum, prec);
+  mpfr_set_prec (sumtrunc, prec);
+
+  mpfr_set (x, arg, GFC_RND_MODE);
+  mpfr_set_ui (sum, 1, GFC_RND_MODE);
+  mpz_set_ui (num, 1);
+
+  mpfr_set (u, x, GFC_RND_MODE);
+  mpfr_sqr (u, u, GFC_RND_MODE);
+  mpfr_mul_ui (u, u, 2, GFC_RND_MODE);
+  mpfr_pow_si (u, u, -1, GFC_RND_MODE);
+
+  for (i = 1; i < MAX_ITER; i++)
+  {
+    mpfr_set (oldsum, sum, GFC_RND_MODE);
+
+    mpz_mul_ui (num, num, 2 * i - 1);
+    mpz_neg (num, num);
+
+    mpfr_set (w, u, GFC_RND_MODE);
+    mpfr_pow_ui (w, w, i, GFC_RND_MODE);
+
+    mpfr_set_z (v, num, GFC_RND_MODE);
+    mpfr_mul (v, v, w, GFC_RND_MODE);
+
+    mpfr_add (sum, sum, v, GFC_RND_MODE);
+
+    mpfr_set (sumtrunc, sum, GFC_RND_MODE);
+    if (mpfr_cmp (sumtrunc, oldsum) == 0)
+      break;
+  }
+
+  /* We should have converged by now; otherwise, ARG_LIMIT is probably
+     set too low.  */
+  gcc_assert (i < MAX_ITER);
+
+  /* Divide by x * sqrt(Pi).  */
+  mpfr_const_pi (u, GFC_RND_MODE);
+  mpfr_sqrt (u, u, GFC_RND_MODE);
+  mpfr_mul (u, u, x, GFC_RND_MODE);
+  mpfr_div (sum, sum, u, GFC_RND_MODE);
+
+  mpfr_set (res, sum, GFC_RND_MODE);
+  mpfr_set_default_prec (prec);
+
+  mpfr_clears (sum, x, u, v, w, oldsum, sumtrunc, NULL);
+  mpz_clear (num);
+}
+
+
+gfc_expr *
+gfc_simplify_erfc_scaled (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  if (mpfr_cmp_d (x->value.real, ARG_LIMIT) >= 0)
+    asympt_erfc_scaled (result->value.real, x->value.real);
+  else
+    fullprec_erfc_scaled (result->value.real, x->value.real);
+
+  return range_check (result, "ERFC_SCALED");
+}
+
+#undef MAX_ITER
+#undef ARG_LIMIT
+
+
 gfc_expr *
 gfc_simplify_epsilon (gfc_expr *e)
 {

gcc/testsuite/ChangeLog

 2009-05-16  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
 
+	PR fortran/33197
+	* gfortran.dg/erf_2.F90: New test.
+	* gfortran.dg/erfc_scaled_2.f90: New test.
+
+2009-05-16  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
+
 	PR fortran/31243
 	* gcc/testsuite/gfortran.dg/string_1.f90: New test.
 	* gcc/testsuite/gfortran.dg/string_2.f90: New test.

gcc/testsuite/gfortran.dg/erf_2.F90

+! { dg-do run }
+! { dg-options "-fno-range-check -ffree-line-length-none " }
+!
+! Check that simplification functions and runtime library agree on ERF,
+! ERFC and ERFC_SCALED.
+
+program test
+  implicit none
+
+  interface check
+    procedure check_r4
+    procedure check_r8
+  end interface check
+
+  real(kind=4) :: x4
+  real(kind=8) :: x8
+
+#define CHECK(a) \
+  x8 = a ; x4 = a ; \
+  call check(erf(real(a,kind=8)), erf(x8)) ; \
+  call check(erf(real(a,kind=4)), erf(x4)) ; \
+  call check(erfc(real(a,kind=8)), erfc(x8)) ; \
+  call check(erfc(real(a,kind=4)), erfc(x4)) ; \
+  call check(erfc_scaled(real(a,kind=8)), erfc_scaled(x8)) ; \
+  call check(erfc_scaled(real(a,kind=4)), erfc_scaled(x4)) ;
+
+  CHECK(0.0)
+  CHECK(0.9)
+  CHECK(1.9)
+  CHECK(19.)
+  CHECK(190.)
+
+  CHECK(-0.0)
+  CHECK(-0.9)
+  CHECK(-1.9)
+  CHECK(-19.)
+  CHECK(-190.)
+
+contains
+
+  subroutine check_r4 (a, b)
+    real(kind=4), intent(in) :: a, b
+    if (abs(a - b) > 10 * spacing(a)) call abort
+  end subroutine
+
+  subroutine check_r8 (a, b)
+    real(kind=8), intent(in) :: a, b
+    if (abs(a - b) > 10 * spacing(a)) call abort
+  end subroutine
+
+end program test

gcc/testsuite/gfortran.dg/erfc_scaled_2.f90

+! { dg-do compile }
+!
+! Check that ERFC_SCALED can be used in initialization expressions
+  real, parameter :: r = 100*erfc_scaled(12.7)
+  integer(kind=int(r)) :: i
+
+  real(kind=8), parameter :: r8 = 100*erfc_scaled(6.77)
+  integer(kind=int(r8)) :: j
+
+  i = 12
+  j = 8
+  print *, i, j
+
+  end