diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 0b8146430a425d2d34a61c4b461c91e5850ccd01..199e7ccb1fe30a8a2c7790f556a2b02e1a281b46 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,11 @@
+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
diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c
index ca125a363357af0170cb5021d56459c710af68f7..612026edadbc42f6ced1b5c19da1814fc6e8125c 100644
--- a/gcc/fortran/intrinsic.c
+++ b/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);
diff --git a/gcc/fortran/intrinsic.h b/gcc/fortran/intrinsic.h
index 83c5207785b7ec52e55a5feb2281b5fde93956c9..7e8bc73ec6fc4b01b3c0e8929323df7454ac99b8 100644
--- a/gcc/fortran/intrinsic.h
+++ b/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 *);
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index 68ebb56d5f227794f4161d4469738cfbee70d562..01b252cf2ad5bdd8ac45480ebc936bec795901bd 100644
--- a/gcc/fortran/simplify.c
+++ b/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)
{
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index 478ba1f0e892cf31164f57dca6e3c78e43a7b04d..41a1338ee7967ec89b96c76456ac4101b3de03d1 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,9 @@
+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
diff --git a/gcc/testsuite/gfortran.dg/erf_2.F90 b/gcc/testsuite/gfortran.dg/erf_2.F90
new file mode 100644
index 0000000000000000000000000000000000000000..d9d3299314e500c462acd3e82dd73e63e53cd0fc
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/erf_2.F90
@@ -0,0 +1,51 @@
+! { 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
diff --git a/gcc/testsuite/gfortran.dg/erfc_scaled_2.f90 b/gcc/testsuite/gfortran.dg/erfc_scaled_2.f90
new file mode 100644
index 0000000000000000000000000000000000000000..97fa91fb915af8282ac7df544c24a5b9bd8bd42c
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/erfc_scaled_2.f90
@@ -0,0 +1,14 @@
+! { 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