diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/digamma.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/digamma.c | 61 |
1 files changed, 35 insertions, 26 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/digamma.c b/Build/source/libs/mpfr/mpfr-src/src/digamma.c index 2cae66b4cd7..e9e8389b115 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/digamma.c +++ b/Build/source/libs/mpfr/mpfr-src/src/digamma.c @@ -1,6 +1,6 @@ /* mpfr_digamma -- digamma function of a floating-point number -Copyright 2009-2019 Free Software Foundation, Inc. +Copyright 2009-2020 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -22,6 +22,10 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #include "mpfr-impl.h" +/* FIXME: Check that MPFR_GET_EXP can only be called on regular values + (in r14025, this is not the case) and that there cannot be integer + overflows. */ + /* Put in s an approximation of digamma(x). Assumes x >= 2. Assumes s does not overlap with x. @@ -36,7 +40,7 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x) mpfr_exp_t e, exps, f, expu; unsigned long n; - MPFR_ASSERTN(MPFR_IS_POS(x) && (MPFR_EXP(x) >= 2)); + MPFR_ASSERTN (MPFR_IS_POS (x) && MPFR_GET_EXP (x) >= 2); mpfr_init2 (t, p); mpfr_init2 (u, p); @@ -44,15 +48,15 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x) mpfr_log (s, x, MPFR_RNDN); /* error <= 1/2 ulp */ mpfr_ui_div (t, 1, x, MPFR_RNDN); /* error <= 1/2 ulp */ - mpfr_div_2exp (t, t, 1, MPFR_RNDN); /* exact */ + mpfr_div_2ui (t, t, 1, MPFR_RNDN); /* exact */ mpfr_sub (s, s, t, MPFR_RNDN); /* error <= 1/2 + 1/2*2^(EXP(olds)-EXP(s)) + 1/2*2^(EXP(t)-EXP(s)). For x >= 2, log(x) >= 2*(1/(2x)), thus olds >= 2t, and olds - t >= olds/2, thus 0 <= EXP(olds)-EXP(s) <= 1, and EXP(t)-EXP(s) <= 0, thus error <= 1/2 + 1/2*2 + 1/2 <= 2 ulps. */ e = 2; /* initial error */ - mpfr_mul (invxx, x, x, MPFR_RNDZ); /* invxx = x^2 * (1 + theta) - for |theta| <= 2^(-p) */ + mpfr_sqr (invxx, x, MPFR_RNDZ); /* invxx = x^2 * (1 + theta) + for |theta| <= 2^(-p) */ mpfr_ui_div (invxx, 1, invxx, MPFR_RNDU); /* invxx = 1/x^2 * (1 + theta)^2 */ /* in the following we note err=xxx when the ratio between the approximation @@ -74,13 +78,13 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x) /* if the terms 'u' are decreasing by a factor two at least, then the error coming from those is bounded by sum((10n+4)/2^n, n=1..infinity) = 24 */ - exps = mpfr_get_exp (s); - expu = mpfr_get_exp (u); + exps = MPFR_GET_EXP (s); + expu = MPFR_GET_EXP (u); if (expu < exps - (mpfr_exp_t) p) break; mpfr_sub (s, s, u, MPFR_RNDN); /* error <= 24 + n/2 */ - if (mpfr_get_exp (s) < exps) - e <<= exps - mpfr_get_exp (s); + if (MPFR_GET_EXP (s) < exps) + e <<= exps - MPFR_GET_EXP (s); e ++; /* error in mpfr_sub */ f = 10 * n + 4; while (expu < exps) @@ -111,9 +115,9 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x) static int mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) { - mpfr_prec_t p = MPFR_PREC(y) + 10, q; + mpfr_prec_t p = MPFR_PREC(y) + 10; mpfr_t t, u, v; - mpfr_exp_t e1, expv; + mpfr_exp_t e1, expv, expx, q; int inex; MPFR_ZIV_DECL (loop); @@ -125,20 +129,22 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) q = PREC(x)-EXP(x) is ok, otherwise if -1 <= x < 0, q = PREC(x)-EXP(x) is ok, otherwise for x < -1, PREC(x) is ok if EXP(x) <= PREC(x), otherwise we need EXP(x) */ - if (MPFR_EXP(x) < 0) - q = MPFR_PREC(x) + 1 - MPFR_EXP(x); - else if (MPFR_EXP(x) <= MPFR_PREC(x)) + expx = MPFR_GET_EXP (x); + if (expx < 0) + q = MPFR_PREC(x) + 1 - expx; + else if (expx <= MPFR_PREC(x)) q = MPFR_PREC(x) + 1; else - q = MPFR_EXP(x); + q = expx; + MPFR_ASSERTN (q <= MPFR_PREC_MAX); mpfr_init2 (u, q); MPFR_DBGRES(inex = mpfr_ui_sub (u, 1, x, MPFR_RNDN)); MPFR_ASSERTN(inex == 0); /* if x is half an integer, cot(Pi*x) = 0, thus Digamma(x) = Digamma(1-x) */ - mpfr_mul_2exp (u, u, 1, MPFR_RNDN); + mpfr_mul_2ui (u, u, 1, MPFR_RNDN); inex = mpfr_integer_p (u); - mpfr_div_2exp (u, u, 1, MPFR_RNDN); + mpfr_div_2ui (u, u, 1, MPFR_RNDN); if (inex) { inex = mpfr_digamma (y, u, rnd_mode); @@ -153,10 +159,10 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) { mpfr_const_pi (v, MPFR_RNDN); /* v = Pi*(1+theta) for |theta|<=2^(-p) */ mpfr_mul (t, v, x, MPFR_RNDN); /* (1+theta)^2 */ - e1 = MPFR_EXP(t) - (mpfr_exp_t) p + 1; /* bound for t: err(t) <= 2^e1 */ + e1 = MPFR_GET_EXP(t) - (mpfr_exp_t) p + 1; /* bound for t: err(t) <= 2^e1 */ mpfr_cot (t, t, MPFR_RNDN); /* cot(t * (1+h)) = cot(t) - theta * (1 + cot(t)^2) with |theta|<=t*h */ - if (MPFR_EXP(t) > 0) + if (MPFR_GET_EXP(t) > 0) e1 = e1 + 2 * MPFR_EXP(t) + 1; else e1 = e1 + 1; @@ -165,9 +171,9 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_mul (t, t, v, MPFR_RNDN); e1 ++; mpfr_digamma (v, u, MPFR_RNDN); /* error <= 1/2 ulp */ - expv = MPFR_EXP(v); + expv = MPFR_GET_EXP (v); mpfr_sub (v, v, t, MPFR_RNDN); - if (MPFR_EXP(v) < MPFR_EXP(t)) + if (MPFR_GET_EXP (v) < MPFR_GET_EXP (t)) e1 += MPFR_EXP(t) - MPFR_EXP(v); /* scale error for t wrt new v */ /* now take into account the 1/2 ulp error for v */ if (expv - MPFR_EXP(v) - 1 > e1) @@ -209,7 +215,7 @@ mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex)); /* compute a precision q such that x+1 is exact */ - if (MPFR_PREC(x) < MPFR_EXP(x)) + if (MPFR_PREC(x) < MPFR_GET_EXP(x)) q = MPFR_EXP(x); else q = MPFR_PREC(x) + 1; @@ -265,11 +271,14 @@ mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) } for (erru = 0; j > 1; erru++, j = (j + 1) / 2); errt = mpfr_digamma_approx (t, x_plus_j); - expt = MPFR_EXP(t); + expt = MPFR_GET_EXP (t); mpfr_sub (t, t, u, MPFR_RNDN); - if (MPFR_EXP(t) < expt) + if (MPFR_GET_EXP (t) < expt) errt += expt - MPFR_EXP(t); - if (MPFR_EXP(t) < MPFR_EXP(u)) + /* Warning: if u is zero (which happens when x_plus_j >= min at the + beginning of the while loop above), EXP(u) is not defined. + In this case we have no error from u. */ + if (MPFR_NOTZERO(u) && MPFR_GET_EXP (t) < MPFR_GET_EXP (u)) erru += MPFR_EXP(u) - MPFR_EXP(t); if (errt > erru) errt = errt + 1; @@ -352,7 +361,7 @@ mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) |y - Digamma(x)| >= 2^(-2n-1)ufp(y), and rounding -1/x gives the correct result. If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1). A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)). */ - if (MPFR_EXP(x) < -2) + if (MPFR_GET_EXP (x) < -2) { if (MPFR_EXP(x) <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y))) { |