diff options
| -rw-r--r-- | src/digamma.c | 27 | ||||
| -rw-r--r-- | tests/tdigamma.c | 21 |
2 files changed, 35 insertions, 13 deletions
diff --git a/src/digamma.c b/src/digamma.c index b971930fa..473093d35 100644 --- a/src/digamma.c +++ b/src/digamma.c @@ -173,16 +173,19 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) mpfr_digamma (v, u, MPFR_RNDN); /* error <= 1/2 ulp */ expv = MPFR_GET_EXP (v); mpfr_sub (v, v, t, MPFR_RNDN); - 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) - e1 = expv - MPFR_EXP(v) - 1; - else - e1 ++; - e1 ++; /* rounding error for mpfr_sub */ - if (MPFR_CAN_ROUND (v, p - e1, MPFR_PREC(y), rnd_mode)) - break; + if (!MPFR_IS_ZERO(v)) + { + 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) + e1 = expv - MPFR_EXP(v) - 1; + else + e1 ++; + e1 ++; /* rounding error for mpfr_sub */ + if (MPFR_CAN_ROUND (v, p - e1, MPFR_PREC(y), rnd_mode)) + break; + } MPFR_ZIV_NEXT (loop, p); mpfr_set_prec (t, p); mpfr_set_prec (v, p); @@ -416,10 +419,8 @@ mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) } } - if (MPFR_IS_NEG(x)) - inex = mpfr_digamma_reflection (y, x, rnd_mode); /* if x < 1/2 we use the reflection formula */ - else if (MPFR_EXP(x) < 0) + if (MPFR_IS_NEG(x) || MPFR_EXP(x) < 0) inex = mpfr_digamma_reflection (y, x, rnd_mode); else inex = mpfr_digamma_positive (y, x, rnd_mode); diff --git a/tests/tdigamma.c b/tests/tdigamma.c index 612e5e098..c24537ee4 100644 --- a/tests/tdigamma.c +++ b/tests/tdigamma.c @@ -110,6 +110,26 @@ bug20210208 (void) mpfr_clear (y); } +/* another test that fails with GMP_CHECK_RANDOMIZE=1613197421465830 + on revision 14429 */ +static void +bug20210215 (void) +{ + mpfr_t x, y; + int inex; + + mpfr_init2 (x, 510); + mpfr_init2 (y, 4); + mpfr_set_str (x, "-8.2923051438433494998166335341807999322052669984208422481227138906096000469898717007386115912802685588348601663465077353194268894939972221117314512518182580e+35", 10, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_digamma (y, x, MPFR_RNDU); + MPFR_ASSERTN (mpfr_cmp_ui (y, 88) == 0); + MPFR_ASSERTN (inex > 0); + MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_INEXACT); + mpfr_clear (x); + mpfr_clear (y); +} + int main (int argc, char *argv[]) { @@ -118,6 +138,7 @@ main (int argc, char *argv[]) special (); bug20210206 (); bug20210208 (); + bug20210215 (); test_generic (MPFR_PREC_MIN, 200, 20); |