aboutsummaryrefslogtreecommitdiff
path: root/gcc
diff options
context:
space:
mode:
authorFrancois-Xavier Coudert <fxcoudert@gcc.gnu.org>2009-05-16 17:33:23 +0000
committerFrançois-Xavier Coudert <fxcoudert@gcc.gnu.org>2009-05-16 17:33:23 +0000
commit9b33a6a1db9e3e21eac8df9909d21d89b57fb507 (patch)
tree4fec6f5d2ae7971100be9fe6ac743f7eb89a125b /gcc
parentb0c068160f502c9d37ec02c8a514546937544eb8 (diff)
downloadgcc-9b33a6a1db9e3e21eac8df9909d21d89b57fb507.zip
gcc-9b33a6a1db9e3e21eac8df9909d21d89b57fb507.tar.gz
gcc-9b33a6a1db9e3e21eac8df9909d21d89b57fb507.tar.bz2
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
Diffstat (limited to 'gcc')
-rw-r--r--gcc/fortran/ChangeLog8
-rw-r--r--gcc/fortran/intrinsic.c5
-rw-r--r--gcc/fortran/intrinsic.h1
-rw-r--r--gcc/fortran/simplify.c137
-rw-r--r--gcc/testsuite/ChangeLog6
-rw-r--r--gcc/testsuite/gfortran.dg/erf_2.F9051
-rw-r--r--gcc/testsuite/gfortran.dg/erfc_scaled_2.f9014
7 files changed, 220 insertions, 2 deletions
diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 0b81464..199e7cc 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,5 +1,13 @@
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.
diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c
index ca125a3..612026e 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 83c5207..7e8bc73 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 68ebb56..01b252c 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 478ba1f..41a1338 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,5 +1,11 @@
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.
diff --git a/gcc/testsuite/gfortran.dg/erf_2.F90 b/gcc/testsuite/gfortran.dg/erf_2.F90
new file mode 100644
index 0000000..d9d3299
--- /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 0000000..97fa91f
--- /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