diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2014-02-04 12:20:47 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2014-02-04 12:20:47 +0000 |
commit | 9330fa00c20dc23e24663096c3028cb9c189ba71 (patch) | |
tree | d72315cf48a4154a08a96c1caa0df9169f4ed4aa /src/digamma.c | |
parent | 514f5dc86b4c9f83cd60b81a5100926985ce1ca3 (diff) | |
download | mpfr-9330fa00c20dc23e24663096c3028cb9c189ba71.tar.gz |
use asymptotic expansion for large positive argument
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@8915 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src/digamma.c')
-rw-r--r-- | src/digamma.c | 20 |
1 files changed, 16 insertions, 4 deletions
diff --git a/src/digamma.c b/src/digamma.c index d3e4ac2ef..9b1c9207c 100644 --- a/src/digamma.c +++ b/src/digamma.c @@ -217,15 +217,27 @@ mpfr_digamma_positive (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)); - /* FIXME: for x large, use the asymptotic expansion: - https://en.wikipedia.org/wiki/Digamma_function#Computation_and_approximation - */ - /* 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); |