summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/digamma.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/digamma.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/digamma.c45
1 files changed, 30 insertions, 15 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/digamma.c b/Build/source/libs/mpfr/mpfr-src/src/digamma.c
index 319eaf22e8b..91e66c54469 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/digamma.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/digamma.c
@@ -34,8 +34,7 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x)
mpfr_prec_t p = MPFR_PREC (s);
mpfr_t t, u, invxx;
mpfr_exp_t e, exps, f, expu;
- mpz_t *INITIALIZED(B); /* variable B declared as initialized */
- unsigned long n0, n; /* number of allocated B[] */
+ unsigned long n;
MPFR_ASSERTN(MPFR_IS_POS(x) && (MPFR_EXP(x) >= 2));
@@ -59,12 +58,9 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x)
/* in the following we note err=xxx when the ratio between the approximation
and the exact result can be written (1 + theta)^xxx for |theta| <= 2^(-p),
following Higham's method */
- B = mpfr_bernoulli_internal ((mpz_t *) 0, 0);
mpfr_set_ui (t, 1, MPFR_RNDN); /* err = 0 */
for (n = 1;; n++)
{
- /* compute next Bernoulli number */
- B = mpfr_bernoulli_internal (B, n);
/* The main term is Bernoulli[2n]/(2n)/x^(2n) = B[n]/(2n+1)!(2n)/x^(2n)
= B[n]*t[n]/(2n) where t[n]/t[n-1] = 1/(2n)/(2n+1)/x^2. */
mpfr_mul (t, t, invxx, MPFR_RNDU); /* err = err + 3 */
@@ -72,7 +68,7 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x)
mpfr_div_ui (t, t, 2 * n + 1, MPFR_RNDU); /* err = err + 1 */
/* we thus have err = 5n here */
mpfr_div_ui (u, t, 2 * n, MPFR_RNDU); /* err = 5n+1 */
- mpfr_mul_z (u, u, B[n], MPFR_RNDU); /* err = 5n+2, and the
+ mpfr_mul_z (u, u, mpfr_bernoulli_cache(n), MPFR_RNDU);/* err = 5n+2, and the
absolute error is bounded
by 10n+4 ulp(u) [Rule 11] */
/* if the terms 'u' are decreasing by a factor two at least,
@@ -92,14 +88,9 @@ mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x)
f = (1 + f) / 2;
expu ++;
}
- e += f; /* total rouding error coming from 'u' term */
+ e += f; /* total rounding error coming from 'u' term */
}
- n0 = ++n;
- while (n--)
- mpz_clear (B[n]);
- (*__gmp_free_func) (B, n0 * sizeof (mpz_t));
-
mpfr_clear (t);
mpfr_clear (u);
mpfr_clear (invxx);
@@ -126,6 +117,10 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
int inex;
MPFR_ZIV_DECL (loop);
+ 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));
+
/* we want that 1-x is exact with precision q: if 0 < x < 1/2, then
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),
@@ -137,7 +132,8 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
else
q = MPFR_EXP(x);
mpfr_init2 (u, q);
- MPFR_ASSERTN(mpfr_ui_sub (u, 1, x, MPFR_RNDN) == 0);
+ 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);
@@ -208,11 +204,31 @@ mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
unsigned long j = 0, min;
MPFR_ZIV_DECL (loop);
+ 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));
+
/* compute a precision q such that x+1 is exact */
if (MPFR_PREC(x) < MPFR_EXP(x))
q = MPFR_EXP(x);
else
q = MPFR_PREC(x) + 1;
+
+ /* for very large x, use |digamma(x) - log(x)| < 1/x < 2^(1-EXP(x)) */
+ if (MPFR_PREC(y) + 10 < MPFR_EXP(x))
+ {
+ /* this ensures EXP(x) >= 3, thus x >= 4, thus log(x) > 1 */
+ mpfr_init2 (t, MPFR_PREC(y) + 10);
+ mpfr_log (t, x, MPFR_RNDZ);
+ if (MPFR_CAN_ROUND (t, MPFR_PREC(y) + 10, MPFR_PREC(y), rnd_mode))
+ {
+ inex = mpfr_set (y, t, rnd_mode);
+ mpfr_clear (t);
+ return inex;
+ }
+ mpfr_clear (t);
+ }
+
mpfr_init2 (x_plus_j, q);
mpfr_init2 (t, p);
@@ -285,7 +301,6 @@ mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
(("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))
@@ -312,7 +327,7 @@ mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
/* the following works also in case of overlap */
MPFR_SET_INF(y);
MPFR_SET_OPPOSITE_SIGN(y, x);
- mpfr_set_divby0 ();
+ MPFR_SET_DIVBY0 ();
MPFR_RET(0);
}
}