| /* mpfr_erf -- error function of a floating-point number |
| |
| Copyright 2001, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. |
| Contributed by Ludovic Meunier and Paul Zimmermann. |
| |
| This file is part of the GNU MPFR Library. |
| |
| The GNU MPFR Library is free software; you can redistribute it and/or modify |
| it under the terms of the GNU Lesser General Public License as published by |
| the Free Software Foundation; either version 2.1 of the License, or (at your |
| option) any later version. |
| |
| The GNU MPFR Library is distributed in the hope that it will be useful, but |
| WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
| or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
| License for more details. |
| |
| You should have received a copy of the GNU Lesser General Public License |
| along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to |
| the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
| MA 02110-1301, USA. */ |
| |
| #define MPFR_NEED_LONGLONG_H |
| #include "mpfr-impl.h" |
| |
| #define EXP1 2.71828182845904523536 /* exp(1) */ |
| |
| static int mpfr_erf_0 (mpfr_ptr, mpfr_srcptr, double, mp_rnd_t); |
| |
| int |
| mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) |
| { |
| mpfr_t xf; |
| int inex, large; |
| MPFR_SAVE_EXPO_DECL (expo); |
| |
| MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode), |
| ("y[%#R]=%R inexact=%d", y, y, inex)); |
| |
| if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) |
| { |
| if (MPFR_IS_NAN (x)) |
| { |
| MPFR_SET_NAN (y); |
| MPFR_RET_NAN; |
| } |
| else if (MPFR_IS_INF (x)) /* erf(+inf) = +1, erf(-inf) = -1 */ |
| return mpfr_set_si (y, MPFR_INT_SIGN (x), GMP_RNDN); |
| else /* erf(+0) = +0, erf(-0) = -0 */ |
| { |
| MPFR_ASSERTD (MPFR_IS_ZERO (x)); |
| return mpfr_set (y, x, GMP_RNDN); /* should keep the sign of x */ |
| } |
| } |
| |
| /* now x is neither NaN, Inf nor 0 */ |
| |
| /* first try expansion at x=0 when x is small, or asymptotic expansion |
| where x is large */ |
| |
| MPFR_SAVE_EXPO_MARK (expo); |
| |
| /* around x=0, we have erf(x) = 2x/sqrt(Pi) (1 - x^2/3 + ...), |
| with 1 - x^2/3 <= sqrt(Pi)*erf(x)/2/x <= 1 for x >= 0. This means that |
| if x^2/3 < 2^(-PREC(y)-1) we can decide of the correct rounding, |
| unless we have a worst-case for 2x/sqrt(Pi). */ |
| if (MPFR_EXP(x) < - (mp_exp_t) (MPFR_PREC(y) / 2)) |
| { |
| /* we use 2x/sqrt(Pi) (1 - x^2/3) <= erf(x) <= 2x/sqrt(Pi) for x > 0 |
| and 2x/sqrt(Pi) <= erf(x) <= 2x/sqrt(Pi) (1 - x^2/3) for x < 0. |
| In both cases |2x/sqrt(Pi) (1 - x^2/3)| <= |erf(x)| <= |2x/sqrt(Pi)|. |
| We will compute l and h such that l <= |2x/sqrt(Pi) (1 - x^2/3)| |
| and |2x/sqrt(Pi)| <= h. If l and h round to the same value to |
| precision PREC(y) and rounding rnd_mode, then we are done. */ |
| mpfr_t l, h; /* lower and upper bounds for erf(x) */ |
| int ok, inex2; |
| |
| mpfr_init2 (l, MPFR_PREC(y) + 17); |
| mpfr_init2 (h, MPFR_PREC(y) + 17); |
| /* first compute l */ |
| mpfr_mul (l, x, x, GMP_RNDU); |
| mpfr_div_ui (l, l, 3, GMP_RNDU); /* upper bound on x^2/3 */ |
| mpfr_ui_sub (l, 1, l, GMP_RNDZ); /* lower bound on 1 - x^2/3 */ |
| mpfr_const_pi (h, GMP_RNDU); /* upper bound of Pi */ |
| mpfr_sqrt (h, h, GMP_RNDU); /* upper bound on sqrt(Pi) */ |
| mpfr_div (l, l, h, GMP_RNDZ); /* lower bound on 1/sqrt(Pi) (1 - x^2/3) */ |
| mpfr_mul_2ui (l, l, 1, GMP_RNDZ); /* 2/sqrt(Pi) (1 - x^2/3) */ |
| mpfr_mul (l, l, x, GMP_RNDZ); /* |l| is a lower bound on |
| |2x/sqrt(Pi) (1 - x^2/3)| */ |
| /* now compute h */ |
| mpfr_const_pi (h, GMP_RNDD); /* lower bound on Pi */ |
| mpfr_sqrt (h, h, GMP_RNDD); /* lower bound on sqrt(Pi) */ |
| mpfr_div_2ui (h, h, 1, GMP_RNDD); /* lower bound on sqrt(Pi)/2 */ |
| /* since sqrt(Pi)/2 < 1, the following should not underflow */ |
| mpfr_div (h, x, h, MPFR_IS_POS(x) ? GMP_RNDU : GMP_RNDD); |
| /* round l and h to precision PREC(y) */ |
| inex = mpfr_prec_round (l, MPFR_PREC(y), rnd_mode); |
| inex2 = mpfr_prec_round (h, MPFR_PREC(y), rnd_mode); |
| /* Caution: we also need inex=inex2 (inex might be 0). */ |
| ok = SAME_SIGN (inex, inex2) && mpfr_cmp (l, h) == 0; |
| if (ok) |
| mpfr_set (y, h, rnd_mode); |
| mpfr_clear (l); |
| mpfr_clear (h); |
| if (ok) |
| goto end; |
| /* this test can still fail for small precision, for example |
| for x=-0.100E-2 with a target precision of 3 bits, since |
| the error term x^2/3 is not that small. */ |
| } |
| |
| mpfr_init2 (xf, 53); |
| mpfr_const_log2 (xf, GMP_RNDU); |
| mpfr_div (xf, x, xf, GMP_RNDZ); /* round to zero ensures we get a lower |
| bound of |x/log(2)| */ |
| mpfr_mul (xf, xf, x, GMP_RNDZ); |
| large = mpfr_cmp_ui (xf, MPFR_PREC (y) + 1) > 0; |
| mpfr_clear (xf); |
| |
| /* when x goes to infinity, we have erf(x) = 1 - 1/sqrt(Pi)/exp(x^2)/x + ... |
| and |erf(x) - 1| <= exp(-x^2) is true for any x >= 0, thus if |
| exp(-x^2) < 2^(-PREC(y)-1) the result is 1 or 1-epsilon. |
| This rewrites as x^2/log(2) > p+1. */ |
| if (MPFR_UNLIKELY (large)) |
| /* |erf x| = 1 or 1- */ |
| { |
| mp_rnd_t rnd2 = MPFR_IS_POS (x) ? rnd_mode : MPFR_INVERT_RND(rnd_mode); |
| if (rnd2 == GMP_RNDN || rnd2 == GMP_RNDU) |
| { |
| inex = MPFR_INT_SIGN (x); |
| mpfr_set_si (y, inex, rnd2); |
| } |
| else /* round to zero */ |
| { |
| inex = -MPFR_INT_SIGN (x); |
| mpfr_setmax (y, 0); /* warning: setmax keeps the old sign of y */ |
| MPFR_SET_SAME_SIGN (y, x); |
| } |
| } |
| else /* use Taylor */ |
| { |
| double xf2; |
| |
| /* FIXME: get rid of doubles/mpfr_get_d here */ |
| xf2 = mpfr_get_d (x, GMP_RNDN); |
| xf2 = xf2 * xf2; /* xf2 ~ x^2 */ |
| inex = mpfr_erf_0 (y, x, xf2, rnd_mode); |
| } |
| |
| end: |
| MPFR_SAVE_EXPO_FREE (expo); |
| return mpfr_check_range (y, inex, rnd_mode); |
| } |
| |
| /* return x*2^e */ |
| static double |
| mul_2exp (double x, mp_exp_t e) |
| { |
| if (e > 0) |
| { |
| while (e--) |
| x *= 2.0; |
| } |
| else |
| { |
| while (e++) |
| x /= 2.0; |
| } |
| |
| return x; |
| } |
| |
| /* evaluates erf(x) using the expansion at x=0: |
| |
| erf(x) = 2/sqrt(Pi) * sum((-1)^k*x^(2k+1)/k!/(2k+1), k=0..infinity) |
| |
| Assumes x is neither NaN nor infinite nor zero. |
| Assumes also that e*x^2 <= n (target precision). |
| */ |
| static int |
| mpfr_erf_0 (mpfr_ptr res, mpfr_srcptr x, double xf2, mp_rnd_t rnd_mode) |
| { |
| mp_prec_t n, m; |
| mp_exp_t nuk, sigmak; |
| double tauk; |
| mpfr_t y, s, t, u; |
| unsigned int k; |
| int log2tauk; |
| int inex; |
| MPFR_ZIV_DECL (loop); |
| |
| n = MPFR_PREC (res); /* target precision */ |
| |
| /* initial working precision */ |
| m = n + (mp_prec_t) (xf2 / LOG2) + 8 + MPFR_INT_CEIL_LOG2 (n); |
| |
| mpfr_init2 (y, m); |
| mpfr_init2 (s, m); |
| mpfr_init2 (t, m); |
| mpfr_init2 (u, m); |
| |
| MPFR_ZIV_INIT (loop, m); |
| for (;;) |
| { |
| mpfr_mul (y, x, x, GMP_RNDU); /* err <= 1 ulp */ |
| mpfr_set_ui (s, 1, GMP_RNDN); |
| mpfr_set_ui (t, 1, GMP_RNDN); |
| tauk = 0.0; |
| |
| for (k = 1; ; k++) |
| { |
| mpfr_mul (t, y, t, GMP_RNDU); |
| mpfr_div_ui (t, t, k, GMP_RNDU); |
| mpfr_div_ui (u, t, 2 * k + 1, GMP_RNDU); |
| sigmak = MPFR_GET_EXP (s); |
| if (k % 2) |
| mpfr_sub (s, s, u, GMP_RNDN); |
| else |
| mpfr_add (s, s, u, GMP_RNDN); |
| sigmak -= MPFR_GET_EXP(s); |
| nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s); |
| |
| if ((nuk < - (mp_exp_t) m) && ((double) k >= xf2)) |
| break; |
| |
| /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */ |
| tauk = 0.5 + mul_2exp (tauk, sigmak) |
| + mul_2exp (1.0 + 8.0 * (double) k, nuk); |
| } |
| |
| mpfr_mul (s, x, s, GMP_RNDU); |
| MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1); |
| |
| mpfr_const_pi (t, GMP_RNDZ); |
| mpfr_sqrt (t, t, GMP_RNDZ); |
| mpfr_div (s, s, t, GMP_RNDN); |
| tauk = 4.0 * tauk + 11.0; /* final ulp-error on s */ |
| log2tauk = __gmpfr_ceil_log2 (tauk); |
| |
| if (MPFR_LIKELY (MPFR_CAN_ROUND (s, m - log2tauk, n, rnd_mode))) |
| break; |
| |
| /* Actualisation of the precision */ |
| MPFR_ZIV_NEXT (loop, m); |
| mpfr_set_prec (y, m); |
| mpfr_set_prec (s, m); |
| mpfr_set_prec (t, m); |
| mpfr_set_prec (u, m); |
| |
| } |
| MPFR_ZIV_FREE (loop); |
| |
| inex = mpfr_set (res, s, rnd_mode); |
| |
| mpfr_clear (y); |
| mpfr_clear (t); |
| mpfr_clear (u); |
| mpfr_clear (s); |
| |
| return inex; |
| } |