diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/erf.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/erf.c | 265 |
1 files changed, 0 insertions, 265 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/erf.c b/Build/source/libs/mpfr/mpfr-src/src/erf.c deleted file mode 100644 index 34884f4346d..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/erf.c +++ /dev/null @@ -1,265 +0,0 @@ -/* mpfr_erf -- error function of a floating-point number - -Copyright 2001, 2003-2020 Free Software Foundation, Inc. -Contributed by the AriC and Caramba projects, INRIA. - -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 3 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.LESSER. If not, see -https://www.gnu.org/licenses/ or 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" - -static int mpfr_erf_0 (mpfr_ptr, mpfr_srcptr, double, mpfr_rnd_t); - -int -mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) -{ - mpfr_t xf; - mp_limb_t xf_limb[(53 - 1) / GMP_NUMB_BITS + 1]; - int inex, large; - MPFR_SAVE_EXPO_DECL (expo); - - MPFR_LOG_FUNC - (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), - ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, 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), MPFR_RNDN); - else /* erf(+0) = +0, erf(-0) = -0 */ - { - MPFR_ASSERTD (MPFR_IS_ZERO (x)); - return mpfr_set (y, x, MPFR_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) < - (mpfr_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_sqr (l, x, MPFR_RNDU); - mpfr_div_ui (l, l, 3, MPFR_RNDU); /* upper bound on x^2/3 */ - mpfr_ui_sub (l, 1, l, MPFR_RNDZ); /* lower bound on 1 - x^2/3 */ - mpfr_const_pi (h, MPFR_RNDU); /* upper bound of Pi */ - mpfr_sqrt (h, h, MPFR_RNDU); /* upper bound on sqrt(Pi) */ - mpfr_div (l, l, h, MPFR_RNDZ); /* lower bound on 1/sqrt(Pi) (1 - x^2/3) */ - mpfr_mul_2ui (l, l, 1, MPFR_RNDZ); /* 2/sqrt(Pi) (1 - x^2/3) */ - mpfr_mul (l, l, x, MPFR_RNDZ); /* |l| is a lower bound on - |2x/sqrt(Pi) (1 - x^2/3)| */ - /* now compute h */ - mpfr_const_pi (h, MPFR_RNDD); /* lower bound on Pi */ - mpfr_sqrt (h, h, MPFR_RNDD); /* lower bound on sqrt(Pi) */ - mpfr_div_2ui (h, h, 1, MPFR_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) ? MPFR_RNDU : MPFR_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_TMP_INIT1(xf_limb, xf, 53); - mpfr_div (xf, x, __gmpfr_const_log2_RNDU, MPFR_RNDZ); /* round to zero - ensures we get a lower bound of |x/log(2)| */ - mpfr_mul (xf, xf, x, MPFR_RNDZ); - large = mpfr_cmp_ui (xf, MPFR_PREC (y) + 1) > 0; - - /* 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- */ - { - mpfr_rnd_t rnd2 = MPFR_IS_POS (x) ? rnd_mode : MPFR_INVERT_RND(rnd_mode); - if (rnd2 == MPFR_RNDN || rnd2 == MPFR_RNDU || rnd2 == MPFR_RNDA) - { - 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, MPFR_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, mpfr_exp_t e) -{ - /* Most of the times, the argument is negative */ - if (MPFR_UNLIKELY (e > 0)) - { - while (e--) - x *= 2.0; - } - else - { - while (e <= -16) - { - x *= (1.0 / 65536.0); - e += 16; - } - while (e++) - x *= 0.5; - } - - 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, mpfr_rnd_t rnd_mode) -{ - mpfr_prec_t n, m; - mpfr_exp_t nuk, sigmak; - mpfr_t y, s, t, u; - unsigned int k; - int inex; - MPFR_GROUP_DECL (group); - MPFR_ZIV_DECL (loop); - - n = MPFR_PREC (res); /* target precision */ - - /* initial working precision */ - m = n + (mpfr_prec_t) (xf2 / LOG2) + 8 + MPFR_INT_CEIL_LOG2 (n); - - MPFR_GROUP_INIT_4(group, m, y, s, t, u); - - MPFR_ZIV_INIT (loop, m); - for (;;) - { - mpfr_t tauk; - mpfr_exp_t log2tauk; - - mpfr_sqr (y, x, MPFR_RNDU); /* err <= 1 ulp */ - mpfr_set_ui (s, 1, MPFR_RNDN); - mpfr_set_ui (t, 1, MPFR_RNDN); - mpfr_init2 (tauk, 53); - mpfr_set_ui (tauk, 0, MPFR_RNDU); - - for (k = 1; ; k++) - { - mpfr_mul (t, y, t, MPFR_RNDU); - mpfr_div_ui (t, t, k, MPFR_RNDU); - mpfr_div_ui (u, t, 2 * k + 1, MPFR_RNDU); - sigmak = MPFR_GET_EXP (s); - if (k % 2) - mpfr_sub (s, s, u, MPFR_RNDN); - else - mpfr_add (s, s, u, MPFR_RNDN); - sigmak -= MPFR_GET_EXP(s); - nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s); - - if ((nuk < - (mpfr_exp_t) m) && (k >= xf2)) - break; - - /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */ - mpfr_mul_2si (tauk, tauk, sigmak, MPFR_RNDU); - mpfr_add_d (tauk, tauk, 0.5 + mul_2exp (1.0 + 8.0 * (double) k, nuk), - MPFR_RNDU); - } - - mpfr_mul (s, x, s, MPFR_RNDU); - MPFR_SET_EXP (s, MPFR_GET_EXP (s) + 1); - - mpfr_const_pi (t, MPFR_RNDZ); - mpfr_sqrt (t, t, MPFR_RNDZ); - mpfr_div (s, s, t, MPFR_RNDN); - /* tauk = 4 * tauk + 11: final ulp-error on s */ - mpfr_mul_2ui (tauk, tauk, 2, MPFR_RNDU); - mpfr_add_ui (tauk, tauk, 11, MPFR_RNDU); - /* In practice, one should not get an infinity, as it would require - a huge precision and a lot of time, but just in case... */ - MPFR_ASSERTN (!MPFR_IS_INF (tauk)); - log2tauk = MPFR_GET_EXP (tauk); - mpfr_clear (tauk); - - if (MPFR_LIKELY (MPFR_CAN_ROUND (s, m - log2tauk, n, rnd_mode))) - break; - - /* Actualisation of the precision */ - MPFR_ZIV_NEXT (loop, m); - MPFR_GROUP_REPREC_4 (group, m, y, s, t, u); - } - MPFR_ZIV_FREE (loop); - - inex = mpfr_set (res, s, rnd_mode); - - MPFR_GROUP_CLEAR (group); - - return inex; -} |