diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-10 23:23:15 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2021-02-10 23:23:15 +0000 |
commit | 32a9b16a85bc52289be725edb06ae6b95cafa9c8 (patch) | |
tree | e64659b3d5b8812126d5d26353afbf4dbb54ecc2 | |
parent | b781448aa0a03d332cbc3c2fa87dbd3a82763684 (diff) | |
download | mpfr-32a9b16a85bc52289be725edb06ae6b95cafa9c8.tar.gz |
[src/digamma.c] Fixed bug in mpfr_digamma_positive when t-u is zero.
[tests/tdigamma.c] Testcase.
(merged changesets r14399,14401-14404 from the trunk)
git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/branches/4.1@14425 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | src/digamma.c | 35 | ||||
-rw-r--r-- | tests/tdigamma.c | 21 |
2 files changed, 41 insertions, 15 deletions
diff --git a/src/digamma.c b/src/digamma.c index 6dba748e6..b971930fa 100644 --- a/src/digamma.c +++ b/src/digamma.c @@ -296,21 +296,26 @@ mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) errt = mpfr_digamma_approx (t, x_plus_j); expt = MPFR_GET_EXP (t); mpfr_sub (t, t, u, MPFR_RNDN); - if (MPFR_GET_EXP (t) < expt) - errt += expt - MPFR_EXP(t); - /* 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; - else if (errt == erru) - errt = errt + 2; - else - errt = erru + 1; - if (MPFR_CAN_ROUND (t, p - errt, MPFR_PREC(y), rnd_mode)) - break; + /* Warning! t may be zero (more likely in small precision). Note + that in this case, this is an exact zero, not an underflow. */ + if (MPFR_NOTZERO(t)) + { + if (MPFR_GET_EXP (t) < expt) + errt += expt - MPFR_EXP(t); + /* 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; + else if (errt == erru) + errt = errt + 2; + else + errt = erru + 1; + if (MPFR_CAN_ROUND (t, p - errt, MPFR_PREC(y), rnd_mode)) + break; + } MPFR_ZIV_NEXT (loop, p); mpfr_set_prec (t, p); mpfr_set_prec (u, p); diff --git a/tests/tdigamma.c b/tests/tdigamma.c index 6841220b0..612e5e098 100644 --- a/tests/tdigamma.c +++ b/tests/tdigamma.c @@ -90,6 +90,26 @@ bug20210206 (void) set_emax (emax); } +/* another test that fails with GMP_CHECK_RANDOMIZE=1612741376857003 + on revision 14398 */ +static void +bug20210208 (void) +{ + mpfr_t x, y; + int inex; + + mpfr_init2 (x, 73); + mpfr_init2 (y, 1); + mpfr_set_str (x, "1.4613470547060071827450", 10, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_digamma (y, x, MPFR_RNDU); + MPFR_ASSERTN (mpfr_cmp_si_2exp (y, -1, -12) == 0); + MPFR_ASSERTN (inex > 0); + MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_INEXACT); + mpfr_clear (x); + mpfr_clear (y); +} + int main (int argc, char *argv[]) { @@ -97,6 +117,7 @@ main (int argc, char *argv[]) special (); bug20210206 (); + bug20210208 (); test_generic (MPFR_PREC_MIN, 200, 20); |