summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/erf.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/erf.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/erf.c265
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;
-}