diff options
Diffstat (limited to 'gcc')
| -rw-r--r-- | gcc/fortran/ChangeLog | 8 | ||||
| -rw-r--r-- | gcc/fortran/intrinsic.c | 5 | ||||
| -rw-r--r-- | gcc/fortran/intrinsic.h | 1 | ||||
| -rw-r--r-- | gcc/fortran/simplify.c | 137 | ||||
| -rw-r--r-- | gcc/testsuite/ChangeLog | 6 | ||||
| -rw-r--r-- | gcc/testsuite/gfortran.dg/erf_2.F90 | 51 | ||||
| -rw-r--r-- | gcc/testsuite/gfortran.dg/erfc_scaled_2.f90 | 14 |
7 files changed, 220 insertions, 2 deletions
diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index 0b8146430a4..199e7ccb1fe 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 ca125a36335..612026edadb 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 83c5207785b..7e8bc73ec6f 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 68ebb56d5f2..01b252cf2ad 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 478ba1f0e89..41a1338ee79 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 00000000000..d9d3299314e --- /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 00000000000..97fa91fb915 --- /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 |