diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2019-06-02 17:05:36.121226916 +0000 +++ mpfr-4.0.2-b/PATCHES 2019-06-02 17:05:36.157226621 +0000 @@ -0,0 +1 @@ +include-float diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2019-01-31 20:29:48.000000000 +0000 +++ mpfr-4.0.2-b/VERSION 2019-06-02 17:05:36.157226621 +0000 @@ -1 +1 @@ -4.0.2 +4.0.2-p1 diff -Naurd mpfr-4.0.2-a/src/mpfr-impl.h mpfr-4.0.2-b/src/mpfr-impl.h --- mpfr-4.0.2-a/src/mpfr-impl.h 2019-01-27 18:30:16.000000000 +0000 +++ mpfr-4.0.2-b/src/mpfr-impl.h 2019-06-02 17:05:36.145226719 +0000 @@ -57,6 +57,7 @@ #include #include +#include /* for FLT_RADIX, etc., tested below */ /****************************************************** diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2019-01-31 20:29:48.000000000 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2019-06-02 17:05:36.153226653 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2" +#define MPFR_VERSION_STRING "4.0.2-p1" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2019-01-31 20:29:48.000000000 +0000 +++ mpfr-4.0.2-b/src/version.c 2019-06-02 17:05:36.153226653 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2"; + return "4.0.2-p1"; } diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-03-30 13:02:08.103248669 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-03-30 13:02:08.143248189 +0000 @@ -0,0 +1 @@ +int-overflow diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2019-06-02 17:05:36.157226621 +0000 +++ mpfr-4.0.2-b/VERSION 2020-03-30 13:02:08.143248189 +0000 @@ -1 +1 @@ -4.0.2-p1 +4.0.2-p2 diff -Naurd mpfr-4.0.2-a/src/agm.c mpfr-4.0.2-b/src/agm.c --- mpfr-4.0.2-a/src/agm.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/agm.c 2020-03-30 13:02:08.131248333 +0000 @@ -220,7 +220,7 @@ mpfr_add (vf, u, v, MPFR_RNDN); /* No overflow? */ mpfr_div_2ui (vf, vf, 1, MPFR_RNDN); /* See proof in algorithms.tex */ - if (4*eq > p) + if (eq > p / 4) { mpfr_t w; MPFR_BLOCK_DECL (flags3); diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2019-06-02 17:05:36.153226653 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-03-30 13:02:08.143248189 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p1" +#define MPFR_VERSION_STRING "4.0.2-p2" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/pow.c mpfr-4.0.2-b/src/pow.c --- mpfr-4.0.2-a/src/pow.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/pow.c 2020-03-30 13:02:08.131248333 +0000 @@ -34,8 +34,7 @@ mpfr_rnd_t rnd_mode, int *inexact) { mpz_t a, c; - mpfr_exp_t d, b; - unsigned long i; + mpfr_exp_t d, b, i; int res; MPFR_ASSERTD (!MPFR_IS_SINGULAR (y)); @@ -48,7 +47,9 @@ if (MPFR_IS_NEG (y)) return 0; /* x is not a power of two => x^-y is not exact */ - /* compute d such that y = c*2^d with c odd integer */ + /* Compute d such that y = c*2^d with c odd integer. + Since c comes from a regular MPFR number, due to the constraints on the + exponent and the precision, there can be no integer overflow below. */ mpz_init (c); d = mpfr_get_z_2exp (c, y); i = mpz_scan1 (c, 0); @@ -58,7 +59,9 @@ /* Since y is not an integer, d is necessarily < 0 */ MPFR_ASSERTD (d < 0); - /* Compute a,b such that x=a*2^b */ + /* Compute a,b such that x=a*2^b. + Since a comes from a regular MPFR number, due to the constrainst on the + exponent and the precision, there can be no integer overflow below. */ mpz_init (a); b = mpfr_get_z_2exp (a, x); i = mpz_scan1 (a, 0); diff -Naurd mpfr-4.0.2-a/src/rem1.c mpfr-4.0.2-b/src/rem1.c --- mpfr-4.0.2-a/src/rem1.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/rem1.c 2020-03-30 13:02:08.131248333 +0000 @@ -100,9 +100,11 @@ mpz_abs (my, my); q_is_odd = 0; - /* divide my by 2^k if possible to make operations mod my easier */ + /* Divide my by 2^k if possible to make operations mod my easier. + Since my comes from a regular MPFR number, due to the constraints on the + exponent and the precision, there can be no integer overflow below. */ { - unsigned long k = mpz_scan1 (my, 0); + mpfr_exp_t k = mpz_scan1 (my, 0); ey += k; mpz_fdiv_q_2exp (my, my, k); } diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2019-06-02 17:05:36.153226653 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-03-30 13:02:08.143248189 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p1"; + return "4.0.2-p2"; } diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-03-30 13:05:16.224990169 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-03-30 13:05:16.264989689 +0000 @@ -0,0 +1 @@ +set-int diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-03-30 13:02:08.143248189 +0000 +++ mpfr-4.0.2-b/VERSION 2020-03-30 13:05:16.264989689 +0000 @@ -1 +1 @@ -4.0.2-p2 +4.0.2-p3 diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-03-30 13:02:08.143248189 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-03-30 13:05:16.260989737 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p2" +#define MPFR_VERSION_STRING "4.0.2-p3" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/set_si_2exp.c mpfr-4.0.2-b/src/set_si_2exp.c --- mpfr-4.0.2-a/src/set_si_2exp.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/set_si_2exp.c 2020-03-30 13:05:16.248989881 +0000 @@ -40,6 +40,15 @@ mp_limb_t ai, *xp; int inex = 0; + /* Early underflow/overflow checking is necessary to avoid + integer overflow or errors due to special exponent values. */ + if (MPFR_UNLIKELY (e < __gmpfr_emin - (mpfr_exp_t) + (sizeof (unsigned long) * CHAR_BIT + 1))) + return mpfr_underflow (x, rnd_mode == MPFR_RNDN ? + MPFR_RNDZ : rnd_mode, i < 0 ? -1 : 1); + if (MPFR_UNLIKELY (e >= __gmpfr_emax)) + return mpfr_overflow (x, rnd_mode, i < 0 ? -1 : 1); + /* FIXME: support int limbs (e.g. 16-bit limbs on 16-bit proc) */ ai = SAFE_ABS (unsigned long, i); MPFR_ASSERTN (SAFE_ABS (unsigned long, i) == ai); diff -Naurd mpfr-4.0.2-a/src/set_ui_2exp.c mpfr-4.0.2-b/src/set_ui_2exp.c --- mpfr-4.0.2-a/src/set_ui_2exp.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/set_ui_2exp.c 2020-03-30 13:05:16.248989881 +0000 @@ -41,6 +41,15 @@ mp_limb_t *xp; int inex = 0; + /* Early underflow/overflow checking is necessary to avoid + integer overflow or errors due to special exponent values. */ + if (MPFR_UNLIKELY (e < __gmpfr_emin - (mpfr_exp_t) + (sizeof (unsigned long) * CHAR_BIT + 1))) + return mpfr_underflow (x, rnd_mode == MPFR_RNDN ? + MPFR_RNDZ : rnd_mode, i < 0 ? -1 : 1); + if (MPFR_UNLIKELY (e >= __gmpfr_emax)) + return mpfr_overflow (x, rnd_mode, i < 0 ? -1 : 1); + /* FIXME: support int limbs (e.g. 16-bit limbs on 16-bit proc) */ MPFR_ASSERTD (i == (mp_limb_t) i); diff -Naurd mpfr-4.0.2-a/src/set_uj.c mpfr-4.0.2-b/src/set_uj.c --- mpfr-4.0.2-a/src/set_uj.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/set_uj.c 2020-03-30 13:05:16.248989881 +0000 @@ -41,7 +41,7 @@ int mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd) { - int cnt; + int cnt, inex; mp_size_t i, k; mp_limb_t limb; mp_limb_t yp[uintmaxpml]; @@ -57,6 +57,10 @@ MPFR_RET(0); } + /* early overflow detection to avoid a possible integer overflow below */ + if (MPFR_UNLIKELY(e >= __gmpfr_emax)) + return mpfr_overflow (x, rnd, MPFR_SIGN_POS); + MPFR_ASSERTN (sizeof(uintmax_t) % sizeof(mp_limb_t) == 0); /* Create an auxiliary var */ @@ -107,7 +111,9 @@ e += k * GMP_NUMB_BITS - cnt; /* Update Expo */ MPFR_ASSERTD (MPFR_LIMB_MSB(yp[numberof (yp) - 1]) != 0); - /* Check expo underflow / overflow (can't use mpfr_check_range) */ + MPFR_RNDRAW (inex, x, yp, uintmax_bit_size, rnd, MPFR_SIGN_POS, e++); + + /* Check expo underflow / overflow */ if (MPFR_UNLIKELY(e < __gmpfr_emin)) { /* The following test is necessary because in the rounding to the @@ -116,16 +122,18 @@ * _ |x| < 2^(emin-2), or * _ |x| = 2^(emin-2) and the absolute value of the exact * result is <= 2^(emin-2). */ - if (rnd == MPFR_RNDN && (e+1 < __gmpfr_emin || mpfr_powerof2_raw(y))) + if (rnd == MPFR_RNDN && + (e + 1 < __gmpfr_emin || + (mpfr_powerof2_raw (x) && inex >= 0))) rnd = MPFR_RNDZ; return mpfr_underflow (x, rnd, MPFR_SIGN_POS); } if (MPFR_UNLIKELY(e > __gmpfr_emax)) return mpfr_overflow (x, rnd, MPFR_SIGN_POS); - MPFR_SET_EXP (y, e); - /* Final: set x to y (rounding if necessary) */ - return mpfr_set (x, y, rnd); + MPFR_SET_SIGN (x, MPFR_SIGN_POS); + MPFR_SET_EXP (x, e); + MPFR_RET (inex); } #endif diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-03-30 13:02:08.143248189 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-03-30 13:05:16.264989689 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p2"; + return "4.0.2-p3"; } diff -Naurd mpfr-4.0.2-a/tests/tset_si.c mpfr-4.0.2-b/tests/tset_si.c --- mpfr-4.0.2-a/tests/tset_si.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tset_si.c 2020-03-30 13:05:16.248989881 +0000 @@ -90,6 +90,135 @@ mpfr_clear (x); } +#define REXP 1024 + +static void +test_2exp_extreme_aux (void) +{ + mpfr_t x1, x2, y; + mpfr_exp_t e, ep[1 + 8 * 5], eb[] = + { MPFR_EMIN_MIN, -REXP, REXP, MPFR_EMAX_MAX, MPFR_EXP_MAX }; + mpfr_flags_t flags1, flags2; + int i, j, rnd, inex1, inex2; + char s; + + ep[0] = MPFR_EXP_MIN; + for (i = 0; i < numberof(eb); i++) + for (j = 0; j < 8; j++) + ep[1 + 8 * i + j] = eb[i] - j; + + mpfr_inits2 (3, x1, x2, (mpfr_ptr) 0); + mpfr_init2 (y, 32); + + for (i = 0; i < numberof(ep); i++) + for (j = -31; j <= 31; j++) + RND_LOOP_NO_RNDF (rnd) + { + int sign = j < 0 ? -1 : 1; + + /* Compute the expected value, inex and flags */ + inex1 = mpfr_set_si (y, j, MPFR_RNDN); + MPFR_ASSERTN (inex1 == 0); + inex1 = mpfr_set (x1, y, (mpfr_rnd_t) rnd); + /* x1 is the rounded value and inex1 the ternary value, + assuming that the exponent argument is 0 (this is the + rounded significand of the final result, assuming an + unbounded exponent range). The multiplication by a + power of 2 is exact, unless underflow/overflow occurs. + The tests on the exponent below avoid integer overflows + (ep[i] may take extreme values). */ + e = mpfr_get_exp (x1); + mpfr_clear_flags (); + if (j != 0 && ep[i] < __gmpfr_emin - e) /* underflow */ + { + mpfr_rnd_t r = + (rnd == MPFR_RNDN && + (ep[i] < __gmpfr_emin - mpfr_get_exp (y) - 1 || + IS_POW2 (sign * j))) ? + MPFR_RNDZ : (mpfr_rnd_t) rnd; + inex1 = mpfr_underflow (x1, r, sign); + flags1 = __gmpfr_flags; + } + else if (j != 0 && ep[i] > __gmpfr_emax - e) /* overflow */ + { + inex1 = mpfr_overflow (x1, (mpfr_rnd_t) rnd, sign); + flags1 = __gmpfr_flags; + } + else + { + if (j != 0) + mpfr_set_exp (x1, ep[i] + e); + flags1 = inex1 != 0 ? MPFR_FLAGS_INEXACT : 0; + } + + /* Test mpfr_set_si_2exp */ + mpfr_clear_flags (); + inex2 = mpfr_set_si_2exp (x2, j, ep[i], (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (x1, x2))) + { + s = 's'; + goto err_extreme; + } + + if (j < 0) + continue; + + /* Test mpfr_set_ui_2exp */ + mpfr_clear_flags (); + inex2 = mpfr_set_ui_2exp (x2, j, ep[i], (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (x1, x2))) + { + s = 'u'; + err_extreme: + printf ("Error in extreme mpfr_set_%ci_2exp for i=%d j=%d %s\n", + s, i, j, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("emin=%" MPFR_EXP_FSPEC "d " + "emax=%" MPFR_EXP_FSPEC "d\n", + (mpfr_eexp_t) __gmpfr_emin, + (mpfr_eexp_t) __gmpfr_emax); + printf ("ep[%d] = %" MPFR_EXP_FSPEC "d\n", + i, (mpfr_eexp_t) ep[i]); + printf ("Expected "); + mpfr_dump (x1); + printf ("with inex = %d and flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (x2); + printf ("with inex = %d and flags =", inex2); + flags_out (flags2); + exit (1); + } + } + + mpfr_clears (x1, x2, y, (mpfr_ptr) 0); +} + +static void +test_2exp_extreme (void) +{ + mpfr_exp_t emin, emax; + + emin = mpfr_get_emin (); + emax = mpfr_get_emax (); + + set_emin (MPFR_EMIN_MIN); + set_emax (MPFR_EMAX_MAX); + test_2exp_extreme_aux (); + + set_emin (-REXP); + set_emax (REXP); + test_2exp_extreme_aux (); + + set_emin (emin); + set_emax (emax); +} + static void test_macros (void) { @@ -639,6 +768,7 @@ mpfr_clear (x); test_2exp (); + test_2exp_extreme (); test_macros (); test_macros_keyword (); test_get_ui_smallneg (); diff -Naurd mpfr-4.0.2-a/tests/tset_sj.c mpfr-4.0.2-b/tests/tset_sj.c --- mpfr-4.0.2-a/tests/tset_sj.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tset_sj.c 2020-03-30 13:05:16.248989881 +0000 @@ -176,6 +176,154 @@ mpfr_clear (x); } +#define REXP 1024 + +static void +test_2exp_extreme_aux (void) +{ + mpfr_t x1, x2, y; + mpfr_exp_t e, ep[1 + 8 * 5], eb[] = + { MPFR_EMIN_MIN, -REXP, REXP, MPFR_EMAX_MAX, MPFR_EXP_MAX }; + mpfr_flags_t flags1, flags2; + int i, j, rnd, inex1, inex2; + char s; + + ep[0] = MPFR_EXP_MIN; + for (i = 0; i < numberof(eb); i++) + for (j = 0; j < 8; j++) + ep[1 + 8 * i + j] = eb[i] - j; + + mpfr_inits2 (3, x1, x2, (mpfr_ptr) 0); + mpfr_init2 (y, 32); + + /* The cast to int is needed if numberof(ep) is of unsigned type, in + which case the condition without the cast would be always false. */ + for (i = -2; i < (int) numberof(ep); i++) + for (j = -31; j <= 31; j++) + RND_LOOP_NO_RNDF (rnd) + { + int sign = j < 0 ? -1 : 1; + intmax_t em; + + if (i < 0) + { + em = i == -2 ? INTMAX_MIN : INTMAX_MAX; + mpfr_clear_flags (); + inex1 = j == 0 ? + (mpfr_set_zero (x1, +1), 0) : i == -2 ? + mpfr_underflow (x1, rnd == MPFR_RNDN ? + MPFR_RNDZ : (mpfr_rnd_t) rnd, sign) : + mpfr_overflow (x1, (mpfr_rnd_t) rnd, sign); + flags1 = __gmpfr_flags; + } + else + { + em = ep[i]; + /* Compute the expected value, inex and flags */ + inex1 = mpfr_set_si (y, j, MPFR_RNDN); + MPFR_ASSERTN (inex1 == 0); + inex1 = mpfr_set (x1, y, (mpfr_rnd_t) rnd); + /* x1 is the rounded value and inex1 the ternary value, + assuming that the exponent argument is 0 (this is the + rounded significand of the final result, assuming an + unbounded exponent range). The multiplication by a + power of 2 is exact, unless underflow/overflow occurs. + The tests on the exponent below avoid integer overflows + (ep[i] may take extreme values). */ + e = mpfr_get_exp (x1); + mpfr_clear_flags (); + if (j != 0 && ep[i] < __gmpfr_emin - e) /* underflow */ + { + mpfr_rnd_t r = + (rnd == MPFR_RNDN && + (ep[i] < __gmpfr_emin - mpfr_get_exp (y) - 1 || + IS_POW2 (sign * j))) ? + MPFR_RNDZ : (mpfr_rnd_t) rnd; + inex1 = mpfr_underflow (x1, r, sign); + flags1 = __gmpfr_flags; + } + else if (j != 0 && ep[i] > __gmpfr_emax - e) /* overflow */ + { + inex1 = mpfr_overflow (x1, (mpfr_rnd_t) rnd, sign); + flags1 = __gmpfr_flags; + } + else + { + if (j != 0) + mpfr_set_exp (x1, ep[i] + e); + flags1 = inex1 != 0 ? MPFR_FLAGS_INEXACT : 0; + } + } + + /* Test mpfr_set_sj_2exp */ + mpfr_clear_flags (); + inex2 = mpfr_set_sj_2exp (x2, j, em, (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (x1, x2))) + { + s = 's'; + goto err_extreme; + } + + if (j < 0) + continue; + + /* Test mpfr_set_uj_2exp */ + mpfr_clear_flags (); + inex2 = mpfr_set_uj_2exp (x2, j, em, (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) && + mpfr_equal_p (x1, x2))) + { + s = 'u'; + err_extreme: + printf ("Error in extreme mpfr_set_%cj_2exp for i=%d j=%d %s\n", + s, i, j, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("emin=%" MPFR_EXP_FSPEC "d " + "emax=%" MPFR_EXP_FSPEC "d\n", + (mpfr_eexp_t) __gmpfr_emin, + (mpfr_eexp_t) __gmpfr_emax); +#ifdef MPFR_PRINTF_MAXLM + printf ("e = %" MPFR_PRINTF_MAXLM "d\n", em); +#endif + printf ("Expected "); + mpfr_dump (x1); + printf ("with inex = %d and flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (x2); + printf ("with inex = %d and flags =", inex2); + flags_out (flags2); + exit (1); + } + } + + mpfr_clears (x1, x2, y, (mpfr_ptr) 0); +} + +static void +test_2exp_extreme (void) +{ + mpfr_exp_t emin, emax; + + emin = mpfr_get_emin (); + emax = mpfr_get_emax (); + + set_emin (MPFR_EMIN_MIN); + set_emax (MPFR_EMAX_MAX); + test_2exp_extreme_aux (); + + set_emin (-REXP); + set_emax (REXP); + test_2exp_extreme_aux (); + + set_emin (emin); + set_emax (emax); +} + int main (int argc, char *argv[]) { @@ -185,6 +333,7 @@ check_set_uj_2exp (); check_set_sj (); check_set_sj_2exp (); + test_2exp_extreme (); tests_end_mpfr (); return 0; diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-03-30 13:09:17.466071979 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-03-30 13:09:17.502071540 +0000 @@ -0,0 +1 @@ +sub1-ubf diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-03-30 13:05:16.264989689 +0000 +++ mpfr-4.0.2-b/VERSION 2020-03-30 13:09:17.502071540 +0000 @@ -1 +1 @@ -4.0.2-p3 +4.0.2-p4 diff -Naurd mpfr-4.0.2-a/src/add1.c mpfr-4.0.2-b/src/add1.c --- mpfr-4.0.2-a/src/add1.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/add1.c 2020-03-30 13:09:17.490071686 +0000 @@ -41,7 +41,7 @@ if (MPFR_UNLIKELY (MPFR_IS_UBF (b))) { - exp = mpfr_ubf_zexp2exp (MPFR_ZEXP (b)); + exp = MPFR_UBF_GET_EXP (b); if (exp > __gmpfr_emax) return mpfr_overflow (a, rnd_mode, MPFR_SIGN (b));; } diff -Naurd mpfr-4.0.2-a/src/mpfr-impl.h mpfr-4.0.2-b/src/mpfr-impl.h --- mpfr-4.0.2-a/src/mpfr-impl.h 2019-06-02 17:05:36.145226719 +0000 +++ mpfr-4.0.2-b/src/mpfr-impl.h 2020-03-30 13:09:17.490071686 +0000 @@ -2406,11 +2406,28 @@ } #endif -#define MPFR_ZEXP(x) \ - ((void) (x)->_mpfr_exp /* to check that x has a correct type */, \ +/* Get the _mpfr_zexp field (pointer to a mpz_t) of a UBF object. + For practical reasons, the type of the argument x can be either + mpfr_ubf_ptr or mpfr_ptr, since the latter is used in functions + that accept both MPFR numbers and UBF's; this is checked by the + code "(x)->_mpfr_exp" (the "sizeof" prevents an access, which + could be invalid when MPFR_ZEXP(x) is used for an assignment, + and also avoids breaking the aliasing rules if they are dealt + with in the future). + This macro can be used when building a UBF. So we do not check + that the _mpfr_exp field has the value MPFR_EXP_UBF. */ +#define MPFR_ZEXP(x) \ + ((void) sizeof ((x)->_mpfr_exp), \ ((mpfr_ubf_ptr) (x))->_mpfr_zexp) +/* If x is a UBF, clear its mpz_t exponent. */ #define MPFR_UBF_CLEAR_EXP(x) \ ((void) (MPFR_IS_UBF (x) && (mpz_clear (MPFR_ZEXP (x)), 0))) +/* Like MPFR_GET_EXP, but accepts UBF (with exponent saturated to + the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]). */ +#define MPFR_UBF_GET_EXP(x) \ + (MPFR_IS_UBF (x) ? mpfr_ubf_zexp2exp (MPFR_ZEXP (x)) : \ + MPFR_GET_EXP ((mpfr_ptr) (x))) + #endif /* __MPFR_IMPL_H__ */ diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-03-30 13:05:16.260989737 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-03-30 13:09:17.498071589 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p3" +#define MPFR_VERSION_STRING "4.0.2-p4" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/sub1.c mpfr-4.0.2-b/src/sub1.c --- mpfr-4.0.2-a/src/sub1.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/sub1.c 2020-03-30 13:09:17.490071686 +0000 @@ -91,9 +91,20 @@ if (MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c))) { - exp_b = MPFR_IS_UBF (b) ? - mpfr_ubf_zexp2exp (MPFR_ZEXP (b)) : MPFR_GET_EXP (b); + exp_b = MPFR_UBF_GET_EXP (b); + /* Early underflow detection. Rare, but a test is needed anyway + since in the "MAX (aq, bq) + 2 <= diff_exp" branch, the exponent + may decrease and MPFR_EXP_MIN would yield an integer overflow. */ + if (MPFR_UNLIKELY (exp_b < __gmpfr_emin - 1)) + { + if (rnd_mode == MPFR_RNDN) + rnd_mode = MPFR_RNDZ; + return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); + } diff_exp = mpfr_ubf_diff_exp (b, c); + /* mpfr_set4 below used with MPFR_RNDF does not support UBF. */ + if (rnd_mode == MPFR_RNDF) + rnd_mode = MPFR_RNDN; } else { @@ -123,11 +134,11 @@ if (rnd_mode == MPFR_RNDF) return mpfr_set4 (a, b, MPFR_RNDZ, MPFR_SIGN (a)); - MPFR_EXP (a) = exp_b; /* may be up to MPFR_EXP_MAX */ + exp_a = exp_b; /* may be any out-of-range value due to UBF */ MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq, rnd_mode, MPFR_SIGN (a), - if (MPFR_EXP (a) != MPFR_EXP_MAX) - ++ MPFR_EXP (a)); + if (exp_a != MPFR_EXP_MAX) + exp_a ++); MPFR_LOG_MSG (("inexact=%d\n", inexact)); if (inexact == 0) { @@ -139,7 +150,7 @@ if (! MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a))) { inexact = MPFR_INT_SIGN (a); - goto check_overflow; + goto end_of_c_small; } } else /* inexact != 0 */ @@ -164,7 +175,7 @@ which means we get a wrong rounded result if x == 1, i.e. inexact == MPFR_EVEN_INEX (for positive numbers). */ if (MPFR_LIKELY (inexact != MPFR_EVEN_INEX * MPFR_INT_SIGN (a))) - goto check_overflow; + goto end_of_c_small; } /* We need to take the value preceding |a|. We can't use mpfr_nexttozero due to a possible out-of-range exponent. @@ -174,16 +185,26 @@ mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh); if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0)) { - MPFR_EXP (a) --; + exp_a --; /* The following is valid whether an = 1 or an > 1. */ ap[an-1] |= MPFR_LIMB_HIGHBIT; } inexact = - MPFR_INT_SIGN (a); - check_overflow: - if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax)) + end_of_c_small: + /* The underflow case is possible only with UBF. The overflow case + is also possible with normal FP due to rounding. */ + if (MPFR_UNLIKELY (exp_a > __gmpfr_emax)) return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)); - else - MPFR_RET (inexact); + if (MPFR_UNLIKELY (exp_a < __gmpfr_emin)) + { + if (rnd_mode == MPFR_RNDN && + (exp_a < __gmpfr_emin - 1 || + (inexact * MPFR_INT_SIGN (a) >= 0 && mpfr_powerof2_raw (a)))) + rnd_mode = MPFR_RNDZ; + return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); + } + MPFR_SET_EXP (a, exp_a); + MPFR_RET (inexact); } /* reserve a space to store b aligned with the result, i.e. shifted by @@ -656,6 +677,15 @@ if (MPFR_LIKELY(cancel)) { cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */ + MPFR_ASSERTD (cancel >= 0); + /* Detect an underflow case to avoid a possible integer overflow + with UBF in the computation of exp_a. */ + if (MPFR_UNLIKELY (exp_b < __gmpfr_emin - 1)) + { + if (rnd_mode == MPFR_RNDN) + rnd_mode = MPFR_RNDZ; + return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); + } exp_a = exp_b - cancel; /* The following assertion corresponds to a limitation of the MPFR implementation. It may fail with a 32-bit ABI and huge precisions, diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-03-30 13:05:16.264989689 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-03-30 13:09:17.502071540 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p3"; + return "4.0.2-p4"; } diff -Naurd mpfr-4.0.2-a/tests/tfma.c mpfr-4.0.2-b/tests/tfma.c --- mpfr-4.0.2-a/tests/tfma.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tfma.c 2020-03-30 13:09:17.490071686 +0000 @@ -47,7 +47,7 @@ mpfr_add (r1, r1, c, (mpfr_rnd_t) rnd)) { if (rnd == MPFR_RNDF) - break; + continue; printf ("test_exact internal error for (%d,%d,%d,%d,%s)\n", i, j, k, rnd, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); exit (1); @@ -664,7 +664,7 @@ Note: The purpose of the s * 2^(emin-7) term is to yield double rounding when scaling for k = 4, s != 0, MPFR_RNDN. */ - RND_LOOP (rnd) + RND_LOOP_NO_RNDF (rnd) { mpfr_clear_flags (); inex1 = mpfr_set_si_2exp (t1, sign * (8*k+s-64), e-6, diff -Naurd mpfr-4.0.2-a/tests/tfmma.c mpfr-4.0.2-b/tests/tfmma.c --- mpfr-4.0.2-a/tests/tfmma.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tfmma.c 2020-03-30 13:09:17.490071686 +0000 @@ -480,6 +480,123 @@ mpfr_clears (x, y, z, (mpfr_ptr) 0); } +static void +underflow_tests (void) +{ + mpfr_t x, y, z; + mpfr_prec_t p; + mpfr_exp_t emin; + + emin = mpfr_get_emin (); + mpfr_set_emin (-17); + for (p = MPFR_PREC_MIN; p <= 1024; p++) + { + mpfr_inits2 (p, x, y, (mpfr_ptr) 0); + mpfr_init2 (z, p + 1); + mpfr_set_str_binary (x, "1e-10"); + mpfr_set_str_binary (y, "1e-11"); + mpfr_clear_underflow (); + mpfr_fmms (z, x, x, y, y, MPFR_RNDN); + /* the exact result is 2^-20-2^-22, and should underflow */ + MPFR_ASSERTN(mpfr_underflow_p ()); + MPFR_ASSERTN(mpfr_zero_p (z)); + MPFR_ASSERTN(mpfr_signbit (z) == 0); + mpfr_clears (x, y, z, (mpfr_ptr) 0); + } + mpfr_set_emin (emin); +} + +static void +bug20180604 (void) +{ + mpfr_t x, y, yneg, z; + mpfr_exp_t emin; + int ret; + + emin = mpfr_get_emin (); + mpfr_set_emin (-1073741821); + mpfr_inits2 (564, x, y, yneg, (mpfr_ptr) 0); + mpfr_init2 (z, 2256); + mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100E-737194993"); + mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010E-737194903"); + + mpfr_clear_underflow (); + ret = mpfr_fmms (z, x, x, y, y, MPFR_RNDN); + MPFR_ASSERTN(mpfr_underflow_p ()); + MPFR_ASSERTN(mpfr_zero_p (z)); + MPFR_ASSERTN(mpfr_signbit (z) == 1); + MPFR_ASSERTN(ret > 0); + + mpfr_clear_underflow (); + ret = mpfr_fmms (z, y, y, x, x, MPFR_RNDN); + MPFR_ASSERTN(mpfr_underflow_p ()); + MPFR_ASSERTN(mpfr_zero_p (z)); + MPFR_ASSERTN(mpfr_signbit (z) == 0); + MPFR_ASSERTN(ret < 0); + + mpfr_neg (yneg, y, MPFR_RNDN); + mpfr_clear_underflow (); + ret = mpfr_fmms (z, x, x, y, yneg, MPFR_RNDN); + MPFR_ASSERTN(mpfr_underflow_p ()); + MPFR_ASSERTN(mpfr_zero_p (z)); + MPFR_ASSERTN(mpfr_signbit (z) == 0); + MPFR_ASSERTN(ret < 0); + + mpfr_clear_underflow (); + ret = mpfr_fmms (z, y, yneg, x, x, MPFR_RNDN); + MPFR_ASSERTN(mpfr_underflow_p ()); + MPFR_ASSERTN(mpfr_zero_p (z)); + MPFR_ASSERTN(mpfr_signbit (z) == 1); + MPFR_ASSERTN(ret > 0); + + mpfr_clears (x, y, yneg, z, (mpfr_ptr) 0); + mpfr_set_emin (emin); +} + +/* this bug was discovered from mpc_mul */ +static void +bug20200206 (void) +{ + mpfr_exp_t emin = mpfr_get_emin (); + mpfr_t xre, xim, yre, yim, zre; + + mpfr_set_emin (-1073); + mpfr_inits2 (53, xre, xim, yre, yim, zre, (mpfr_ptr) 0); + mpfr_set_str (xre, "-0x8.294611b331c8p-904", 16, MPFR_RNDN); + mpfr_set_str (xim, "-0x1.859278c2992fap-676", 16, MPFR_RNDN); + mpfr_set_str (yre, "0x9.ac54802a95f8p-820", 16, MPFR_RNDN); + mpfr_set_str (yim, "0x3.17e59e7612aap-412", 16, MPFR_RNDN); + mpfr_fmms (zre, xre, yre, xim, yim, MPFR_RNDN); + if (mpfr_regular_p (zre) && mpfr_get_exp (zre) < -1073) + { + printf ("Error, mpfr_fmms returns an out-of-range exponent:\n"); + mpfr_dump (zre); + exit (1); + } + mpfr_clears (xre, xim, yre, yim, zre, (mpfr_ptr) 0); + mpfr_set_emin (emin); +} + +/* check for integer overflow (see bug fixed in r13798) */ +static void +extreme_underflow (void) +{ + mpfr_t x, y, z; + mpfr_exp_t emin = mpfr_get_emin (); + + set_emin (MPFR_EMIN_MIN); + mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0); + mpfr_set_ui_2exp (x, 1, MPFR_EMIN_MIN - 1, MPFR_RNDN); + mpfr_set (y, x, MPFR_RNDN); + mpfr_nextabove (x); + mpfr_clear_flags (); + mpfr_fmms (z, x, x, y, y, MPFR_RNDN); + MPFR_ASSERTN (__gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT)); + MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z)); + mpfr_clears (x, y, z, (mpfr_ptr) 0); + set_emin (emin); +} + /* Test double-rounding cases in mpfr_set_1_2, which is called when all the precisions are the same. With set.c before r13347, this triggers errors for neg=0. */ @@ -548,6 +665,9 @@ { tests_start_mpfr (); + bug20200206 (); + bug20180604 (); + underflow_tests (); random_tests (); zero_tests (); max_tests (); @@ -555,6 +675,7 @@ half_plus_half (); bug20170405 (); double_rounding (); + extreme_underflow (); tests_end_mpfr (); return 0; diff -Naurd mpfr-4.0.2-a/tests/tsub.c mpfr-4.0.2-b/tests/tsub.c --- mpfr-4.0.2-a/tests/tsub.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tsub.c 2020-03-30 13:09:17.490071686 +0000 @@ -1159,6 +1159,358 @@ } } +/* Tests on UBF. + * + * Note: mpfr_sub1sp will never be used as it does not support UBF. + * Thus there is no need to generate tests for both mpfr_sub1 and + * mpfr_sub1sp. + * + * Note that mpfr_sub1 has a special branch "c small", where the second + * argument c is sufficiently smaller than the ulp of the first argument + * and the ulp of the result: MAX (aq, bq) + 2 <= diff_exp. + * Tests should be done for both the main branch and this special branch + * when this makes sense. + */ +#define REXP 1024 + +static void test_ubf_aux (void) +{ + mpfr_ubf_t x[11]; + mpfr_ptr p[11]; + int ex[11]; + mpfr_t ee, y, z, w; + int i, j, k, neg, inexact, rnd; + const int kn = 2; + mpfr_exp_t e[] = + { MPFR_EXP_MIN, MPFR_EMIN_MIN, -REXP, 0, + REXP, MPFR_EMAX_MAX, MPFR_EXP_MAX }; + + mpfr_init2 (ee, sizeof (mpfr_exp_t) * CHAR_BIT); + mpfr_inits2 (64, y, z, (mpfr_ptr) 0); + mpfr_init2 (w, 2); + + for (i = 0; i < 11; i++) + p[i] = (mpfr_ptr) x[i]; + + /* exact zero result, with small and large exponents */ + for (i = 0; i < 2; i++) + { + mpfr_init2 (p[i], 5 + (randlimb () % 128)); + mpfr_set_ui (p[i], 17, MPFR_RNDN); + mpz_init (MPFR_ZEXP (p[i])); + MPFR_SET_UBF (p[i]); + } + for (j = 0; j < numberof (e); j++) + { + inexact = mpfr_set_exp_t (ee, e[j], MPFR_RNDN); + MPFR_ASSERTD (inexact == 0); + inexact = mpfr_get_z (MPFR_ZEXP (p[0]), ee, MPFR_RNDN); + MPFR_ASSERTD (inexact == 0); + mpz_sub_ui (MPFR_ZEXP (p[0]), MPFR_ZEXP (p[0]), kn); + + for (k = -kn; k <= kn; k++) + { + /* exponent: e[j] + k, with |k| <= kn */ + mpz_set (MPFR_ZEXP (p[1]), MPFR_ZEXP (p[0])); + + for (neg = 0; neg <= 1; neg++) + { + RND_LOOP (rnd) + { + /* Note: x[0] and x[1] are equal MPFR numbers, but do not + test mpfr_sub with arg2 == arg3 as pointers in order to + skip potentially optimized mpfr_sub code. */ + inexact = mpfr_sub (z, p[0], p[1], (mpfr_rnd_t) rnd); + if (inexact != 0 || MPFR_NOTZERO (z) || + (rnd != MPFR_RNDD ? MPFR_IS_NEG (z) : MPFR_IS_POS (z))) + { + printf ("Error 1 in test_ubf for exact zero result: " + "j=%d k=%d neg=%d, rnd=%s\nGot ", j, k, neg, + mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + mpfr_dump (z); + printf ("inexact = %d\n", inexact); + exit (1); + } + } + + for (i = 0; i < 2; i++) + MPFR_CHANGE_SIGN (p[i]); + } + + mpz_add_ui (MPFR_ZEXP (p[0]), MPFR_ZEXP (p[0]), 1); + } + } + for (i = 0; i < 2; i++) + { + MPFR_UBF_CLEAR_EXP (p[i]); + mpfr_clear (p[i]); + } + + /* Up to a given exponent (for the result) and sign, test: + * (t + .11010) - (t + .00001) = .11001 + * (t + 8) - (t + 111.00111) = .11001 + * where t = 0 or a power of 2, e.g. 2^200. Test various exponents + * (including those near the underflow/overflow boundaries) so that + * the subtraction yields a normal number, an overflow or an underflow. + * In MPFR_RNDA, also test with a 2-bit precision target, as this + * yields an exponent change. + * + * Also test the "MAX (aq, bq) + 2 <= diff_exp" branch of sub1.c with + * .1 - epsilon (possible decrease of the exponent) and .111 - epsilon + * in precision 2 (possible increase of the exponent). The first test + * triggers a possible decrease of the exponent (see bug fixed in r13806). + * The second test triggers a possible increase of the exponent (see the + * "exp_a != MPFR_EXP_MAX" test to avoid an integer overflow). + */ + for (i = 0; i < 8; i++) + { + static int v[4] = { 26, 1, 256, 231 }; + + mpfr_init2 (p[i], i < 4 ? 5 + (randlimb () % 128) : 256); + if (i < 4) + mpfr_set_si_2exp (p[i], v[i], -5, MPFR_RNDN); + else + { + mpfr_set_si_2exp (p[i], 1, 200, MPFR_RNDN); + mpfr_add (p[i], p[i], p[i-4], MPFR_RNDN); + } + ex[i] = mpfr_get_exp (p[i]) + 5; + MPFR_ASSERTD (ex[i] >= 0); + } + mpfr_inits2 (3, p[8], p[9], p[10], (mpfr_ptr) 0); + mpfr_set_si_2exp (p[8], 1, 0, MPFR_RNDN); + ex[8] = 5; + mpfr_set_si_2exp (p[9], 1, 0, MPFR_RNDN); /* will be epsilon */ + ex[9] = 0; + mpfr_set_si_2exp (p[10], 7, 0, MPFR_RNDN); + ex[10] = 5; + + for (i = 0; i < 11; i++) + { + mpz_init (MPFR_ZEXP (p[i])); + MPFR_SET_UBF (p[i]); + } + + for (j = 0; j < numberof (e); j++) + { + inexact = mpfr_set_exp_t (ee, e[j], MPFR_RNDN); + MPFR_ASSERTD (inexact == 0); + inexact = mpfr_get_z (MPFR_ZEXP (p[0]), ee, MPFR_RNDN); + MPFR_ASSERTD (inexact == 0); + for (i = 1; i < 11; i++) + mpz_set (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[0])); + for (i = 0; i < 11; i++) + { + mpz_add_ui (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[i]), ex[i]); + mpz_sub_ui (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[i]), 5 + kn); + } + mpz_sub_ui (MPFR_ZEXP (p[9]), MPFR_ZEXP (p[9]), 256); + for (k = -kn; k <= kn; k++) + { + for (neg = 0; neg <= 1; neg++) + { + int sign = neg ? -1 : 1; + + RND_LOOP (rnd) + for (i = 0; i <= 10; i += 2) + { + mpfr_exp_t e0; + mpfr_flags_t flags, flags_y; + int inex_y; + + if (i >= 8) + { + int d; + + e0 = MPFR_UBF_GET_EXP (p[i]); + if (e0 < MPFR_EXP_MIN + 3) + e0 += 3; + + if (rnd == MPFR_RNDN) + d = i == 8 ? (e0 == __gmpfr_emin - 1 ? 3 : 4) : 6; + else if (MPFR_IS_LIKE_RNDZ (rnd, neg)) + d = i == 8 ? 3 : 6; + else + d = i == 8 ? 4 : 8; + + mpfr_clear_flags (); + inex_y = mpfr_set_si_2exp (w, sign * d, e0 - 3, + (mpfr_rnd_t) rnd); + flags_y = __gmpfr_flags | MPFR_FLAGS_INEXACT; + if (inex_y == 0) + inex_y = rnd == MPFR_RNDN ? + sign * (i == 8 ? 1 : -1) : + MPFR_IS_LIKE_RNDD ((mpfr_rnd_t) rnd, sign) ? + -1 : 1; + mpfr_set (y, w, MPFR_RNDN); + + mpfr_clear_flags (); + inexact = mpfr_sub (w, p[i], p[9], (mpfr_rnd_t) rnd); + flags = __gmpfr_flags; + + /* For MPFR_RNDF, only do a basic test. */ + MPFR_ASSERTN (mpfr_check (w)); + if (rnd == MPFR_RNDF) + continue; + + goto testw; + } + + mpfr_clear_flags (); + inexact = mpfr_sub (z, p[i], p[i+1], (mpfr_rnd_t) rnd); + flags = __gmpfr_flags; + + /* For MPFR_RNDF, only do a basic test. */ + MPFR_ASSERTN (mpfr_check (z)); + if (rnd == MPFR_RNDF) + continue; + + e0 = MPFR_UBF_GET_EXP (p[0]); + + if (e0 < __gmpfr_emin) + { + mpfr_rnd_t r = + rnd == MPFR_RNDN && e0 < __gmpfr_emin - 1 ? + MPFR_RNDZ : (mpfr_rnd_t) rnd; + flags_y = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; + inex_y = mpfr_underflow (y, r, sign); + } + else if (e0 > __gmpfr_emax) + { + flags_y = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; + inex_y = mpfr_overflow (y, (mpfr_rnd_t) rnd, sign); + } + else + { + mpfr_set_si_2exp (y, sign * 25, e0 - 5, MPFR_RNDN); + flags_y = 0; + inex_y = 0; + } + + if (flags != flags_y || + ! SAME_SIGN (inexact, inex_y) || + ! mpfr_equal_p (y, z)) + { + printf ("Error 2 in test_ubf with " + "j=%d k=%d neg=%d i=%d rnd=%s\n", + j, k, neg, i, + mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("emin=%" MPFR_EXP_FSPEC "d " + "emax=%" MPFR_EXP_FSPEC "d\n", + (mpfr_eexp_t) __gmpfr_emin, + (mpfr_eexp_t) __gmpfr_emax); + printf ("b = "); + mpfr_dump (p[i]); + printf ("c = "); + mpfr_dump (p[i+1]); + printf ("Expected "); + mpfr_dump (y); + printf ("with inex = %d and flags =", inex_y); + flags_out (flags_y); + printf ("Got "); + mpfr_dump (z); + printf ("with inex = %d and flags =", inexact); + flags_out (flags); + exit (1); + } + + /* Do the following 2-bit precision test only in RNDA. */ + if (rnd != MPFR_RNDA) + continue; + + mpfr_clear_flags (); + inexact = mpfr_sub (w, p[i], p[i+1], MPFR_RNDA); + flags = __gmpfr_flags; + if (e0 < MPFR_EXP_MAX) + e0++; + + if (e0 < __gmpfr_emin) + { + flags_y = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT; + inex_y = mpfr_underflow (y, MPFR_RNDA, sign); + } + else if (e0 > __gmpfr_emax) + { + flags_y = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT; + inex_y = mpfr_overflow (y, MPFR_RNDA, sign); + } + else + { + mpfr_set_si_2exp (y, sign, e0 - 1, MPFR_RNDN); + flags_y = MPFR_FLAGS_INEXACT; + inex_y = sign; + } + + testw: + if (flags != flags_y || + ! SAME_SIGN (inexact, inex_y) || + ! mpfr_equal_p (y, w)) + { + printf ("Error 3 in test_ubf with " + "j=%d k=%d neg=%d i=%d rnd=%s\n", + j, k, neg, i, + mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("emin=%" MPFR_EXP_FSPEC "d " + "emax=%" MPFR_EXP_FSPEC "d\n", + (mpfr_eexp_t) __gmpfr_emin, + (mpfr_eexp_t) __gmpfr_emax); + printf ("b = "); + mpfr_dump (p[i]); + printf ("c = "); + mpfr_dump (p[i <= 8 ? i+1 : 9]); + printf ("Expected "); + /* Set y to a 2-bit precision just for the output. + Since we exit, this will have no other effect. */ + mpfr_prec_round (y, 2, MPFR_RNDA); + mpfr_dump (y); + printf ("with inex = %d and flags =", inex_y); + flags_out (flags_y); + printf ("Got "); + mpfr_dump (w); + printf ("with inex = %d and flags =", inexact); + flags_out (flags); + exit (1); + } + } + + for (i = 0; i < 11; i++) + MPFR_CHANGE_SIGN (p[i]); + } + + for (i = 0; i < 11; i++) + mpz_add_ui (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[i]), 1); + } + } + for (i = 0; i < 11; i++) + { + MPFR_UBF_CLEAR_EXP (p[i]); + mpfr_clear (p[i]); + } + + mpfr_clears (ee, y, z, w, (mpfr_ptr) 0); +} + +/* Run the tests on UBF with the maximum exponent range and with a + reduced exponent range. */ +static void test_ubf (void) +{ + mpfr_exp_t emin, emax; + + emin = mpfr_get_emin (); + emax = mpfr_get_emax (); + + set_emin (MPFR_EMIN_MIN); + set_emax (MPFR_EMAX_MAX); + test_ubf_aux (); + + set_emin (-REXP); + set_emax (REXP); + test_ubf_aux (); + + set_emin (emin); + set_emax (emax); +} + #define TEST_FUNCTION test_sub #define TWO_ARGS #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) @@ -1188,6 +1540,7 @@ for (i=0; i<50; i++) check_two_sum (p); test_generic (MPFR_PREC_MIN, 800, 100); + test_ubf (); tests_end_mpfr (); return 0; diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-03-30 16:31:45.813415121 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-03-30 16:33:54.623893528 +0000 @@ -0,0 +1 @@ +const diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-03-30 13:09:17.502071540 +0000 +++ mpfr-4.0.2-b/VERSION 2020-03-30 15:17:31.535330224 +0000 @@ -1 +1 @@ -4.0.2-p4 +4.0.2-p5 diff -Naurd mpfr-4.0.2-a/acinclude.m4 mpfr-4.0.2-b/acinclude.m4 --- mpfr-4.0.2-a/acinclude.m4 2019-01-27 18:30:16.000000000 +0000 +++ mpfr-4.0.2-b/acinclude.m4 2019-01-27 18:30:16.000000000 +0000 @@ -438,12 +438,14 @@ AC_MSG_WARN([platform and/or document the behavior.]) fi -dnl Check if the chars '0' to '9' are consecutive values +dnl Check if the chars '0' to '9', 'a' to 'z', and 'A' to 'Z' are +dnl consecutive values. +dnl The const is necessary with GCC's "-Wwrite-strings -Werror". AC_MSG_CHECKING([if charset has consecutive values]) AC_RUN_IFELSE([AC_LANG_PROGRAM([[ -char *number = "0123456789"; -char *lower = "abcdefghijklmnopqrstuvwxyz"; -char *upper = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"; +const char *number = "0123456789"; +const char *lower = "abcdefghijklmnopqrstuvwxyz"; +const char *upper = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"; ]],[[ int i; unsigned char *p; diff -Naurd mpfr-4.0.2-a/configure mpfr-4.0.2-b/configure --- mpfr-4.0.2-a/configure 2019-01-31 20:43:19.000000000 +0000 +++ mpfr-4.0.2-b/configure 2020-03-30 16:33:29.340190160 +0000 @@ -15868,9 +15868,9 @@ cat confdefs.h - <<_ACEOF >conftest.$ac_ext /* end confdefs.h. */ -char *number = "0123456789"; -char *lower = "abcdefghijklmnopqrstuvwxyz"; -char *upper = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"; +const char *number = "0123456789"; +const char *lower = "abcdefghijklmnopqrstuvwxyz"; +const char *upper = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"; int main (void) diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-03-30 13:09:17.498071589 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-03-30 15:17:31.535330224 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p4" +#define MPFR_VERSION_STRING "4.0.2-p5" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-03-30 13:09:17.502071540 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-03-30 15:17:31.535330224 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p4"; + return "4.0.2-p5"; } diff -Naurd mpfr-4.0.2-a/tests/tdiv.c mpfr-4.0.2-b/tests/tdiv.c --- mpfr-4.0.2-a/tests/tdiv.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tdiv.c 2020-03-30 15:17:31.523330366 +0000 @@ -23,7 +23,7 @@ #include "mpfr-test.h" static void -check_equal (mpfr_srcptr a, mpfr_srcptr a2, char *s, +check_equal (mpfr_srcptr a, mpfr_srcptr a2, const char *s, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t r) { if (SAME_VAL (a, a2)) diff -Naurd mpfr-4.0.2-a/tests/tfpif.c mpfr-4.0.2-b/tests/tfpif.c --- mpfr-4.0.2-a/tests/tfpif.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tfpif.c 2020-03-30 15:17:31.523330366 +0000 @@ -30,8 +30,8 @@ static void doit (int argc, char *argv[], mpfr_prec_t p1, mpfr_prec_t p2) { - char *filenameCompressed = FILE_NAME_RW; - char *data = FILE_NAME_R; + const char *filenameCompressed = FILE_NAME_RW; + const char *data = FILE_NAME_R; int status; FILE *fh; mpfr_t x[9]; @@ -198,7 +198,7 @@ static void check_bad (void) { - char *filenameCompressed = FILE_NAME_RW; + const char *filenameCompressed = FILE_NAME_RW; int status; FILE *fh; mpfr_t x; @@ -340,7 +340,7 @@ static void extra (void) { - char *data = FILE_NAME_R2; + const char *data = FILE_NAME_R2; mpfr_t x; FILE *fp; int ret; diff -Naurd mpfr-4.0.2-a/tests/tmul_2exp.c mpfr-4.0.2-b/tests/tmul_2exp.c --- mpfr-4.0.2-a/tests/tmul_2exp.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tmul_2exp.c 2020-03-30 15:17:31.523330366 +0000 @@ -258,7 +258,7 @@ mpfr_exp_t old_emax; mpfr_t x, y1, y2; int neg, r, op; - static char *sop[4] = { "mul_2ui", "mul_2si", "div_2ui", "div_2si" }; + static const char *sop[4] = { "mul_2ui", "mul_2si", "div_2ui", "div_2si" }; old_emax = mpfr_get_emax (); set_emax (emax); diff -Naurd mpfr-4.0.2-a/tests/trandom.c mpfr-4.0.2-b/tests/trandom.c --- mpfr-4.0.2-a/tests/trandom.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/trandom.c 2020-03-30 15:17:31.523330366 +0000 @@ -185,7 +185,7 @@ { gmp_randstate_t s; mpfr_t x; - char *str = "0.1010111100000000000000000000000000000000E-32"; + const char *str = "0.1010111100000000000000000000000000000000E-32"; int k; gmp_randinit_default (s); diff -Naurd mpfr-4.0.2-a/tests/trint.c mpfr-4.0.2-b/tests/trint.c --- mpfr-4.0.2-a/tests/trint.c 2019-01-25 02:00:06.000000000 +0000 +++ mpfr-4.0.2-b/tests/trint.c 2020-03-30 15:17:31.523330366 +0000 @@ -308,7 +308,7 @@ #if __MPFR_STDC (199901L) static void -test_fct (double (*f)(double), int (*g)(), char *s, mpfr_rnd_t r) +test_fct (double (*f)(double), int (*g)(), const char *s, mpfr_rnd_t r) { double d, y; mpfr_t dd, yy; diff -Naurd mpfr-4.0.2-a/tests/tsum.c mpfr-4.0.2-b/tests/tsum.c --- mpfr-4.0.2-a/tests/tsum.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tsum.c 2020-03-30 15:17:31.523330366 +0000 @@ -285,7 +285,7 @@ static void check_more_special (void) { - char *str[NC] = { "NaN", "+Inf", "-Inf", "+0", "-0", "+1", "-1" }; + const char *str[NC] = { "NaN", "+Inf", "-Inf", "+0", "-0", "+1", "-1" }; int i, r, k[NS]; mpfr_t c[NC], s[NS], sum; mpfr_ptr p[NS]; @@ -724,7 +724,7 @@ { mpfr_t sum, t[4]; mpfr_ptr p[4]; - char *s[4] = { + const char *s[4] = { "0x1p1000", "-0x0.fffffffffffff80000000000000001p1000", "-0x1p947", @@ -771,7 +771,11 @@ { mpfr_t sum1, sum2, t[3]; mpfr_ptr p[3]; - char *s[3] = { "0.10000111110101000010101011100001", "1E-100", "0.1E95" }; + const char *s[3] = { + "0.10000111110101000010101011100001", + "1E-100", + "0.1E95" + }; int i, r; mpfr_inits2 (58, sum1, sum2, (mpfr_ptr) 0); @@ -1268,7 +1272,7 @@ { mpfr_t r, t[4]; mpfr_ptr p[4]; - char *s[4] = { "0.10E20", "0", "0.11E382", "0.10E826" }; + const char *s[4] = { "0.10E20", "0", "0.11E382", "0.10E826" }; int i; mpfr_init2 (r, 2); diff -Naurd mpfr-4.0.2-a/tests/turandom.c mpfr-4.0.2-b/tests/turandom.c --- mpfr-4.0.2-a/tests/turandom.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/turandom.c 2020-03-30 15:17:31.523330366 +0000 @@ -526,7 +526,7 @@ #define N 6 /* Run this program with the MPFR_REPROD_ABI_OUTPUT environment variable set to get the array of strings. */ - char *t[5 * N] = { + const char *t[5 * N] = { "1.0@-1", "3.0@-1", "7.0@-1", diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-03-30 16:50:17.064231191 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-03-30 16:50:17.096230810 +0000 @@ -0,0 +1 @@ +array-length diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-03-30 15:17:31.535330224 +0000 +++ mpfr-4.0.2-b/VERSION 2020-03-30 16:50:17.096230810 +0000 @@ -1 +1 @@ -4.0.2-p5 +4.0.2-p6 diff -Naurd mpfr-4.0.2-a/src/mpfr-impl.h mpfr-4.0.2-b/src/mpfr-impl.h --- mpfr-4.0.2-a/src/mpfr-impl.h 2020-03-30 13:09:17.490071686 +0000 +++ mpfr-4.0.2-b/src/mpfr-impl.h 2020-03-30 16:50:17.084230953 +0000 @@ -2026,7 +2026,21 @@ struct mpfr_group_t { size_t alloc; mp_limb_t *mant; +#if MPFR_GROUP_STATIC_SIZE != 0 mp_limb_t tab[MPFR_GROUP_STATIC_SIZE]; +#else + /* In order to detect memory leaks when testing, MPFR_GROUP_STATIC_SIZE + can be set to 0, in which case tab will not be used. ISO C does not + support zero-length arrays[*], thus let's use a flexible array member + (which will be equivalent here). Note: this is new in C99, but this + is just used for testing. + [*] Zero-length arrays are a GNU extension: + https://gcc.gnu.org/onlinedocs/gcc/Zero-Length.html + and as such an extension is forbidden in ISO C, it triggers an + error with -Werror=pedantic. + */ + mp_limb_t tab[]; +#endif }; #define MPFR_GROUP_DECL(g) struct mpfr_group_t g diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-03-30 15:17:31.535330224 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-03-30 16:50:17.096230810 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p5" +#define MPFR_VERSION_STRING "4.0.2-p6" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-03-30 15:17:31.535330224 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-03-30 16:50:17.096230810 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p5"; + return "4.0.2-p6"; } diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-04-03 13:54:03.807946879 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-04-03 13:54:03.891945830 +0000 @@ -0,0 +1 @@ +sub1-ubftest diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-03-30 16:50:17.096230810 +0000 +++ mpfr-4.0.2-b/VERSION 2020-04-03 13:54:03.891945830 +0000 @@ -1 +1 @@ -4.0.2-p6 +4.0.2-p7 diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-03-30 16:50:17.096230810 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-04-03 13:54:03.891945830 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p6" +#define MPFR_VERSION_STRING "4.0.2-p7" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-03-30 16:50:17.096230810 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-04-03 13:54:03.891945830 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p6"; + return "4.0.2-p7"; } diff -Naurd mpfr-4.0.2-a/tests/tsub.c mpfr-4.0.2-b/tests/tsub.c --- mpfr-4.0.2-a/tests/tsub.c 2020-03-30 13:09:17.490071686 +0000 +++ mpfr-4.0.2-b/tests/tsub.c 2020-04-03 13:54:03.827946629 +0000 @@ -1203,9 +1203,9 @@ for (j = 0; j < numberof (e); j++) { inexact = mpfr_set_exp_t (ee, e[j], MPFR_RNDN); - MPFR_ASSERTD (inexact == 0); + MPFR_ASSERTN (inexact == 0); inexact = mpfr_get_z (MPFR_ZEXP (p[0]), ee, MPFR_RNDN); - MPFR_ASSERTD (inexact == 0); + MPFR_ASSERTN (inexact == 0); mpz_sub_ui (MPFR_ZEXP (p[0]), MPFR_ZEXP (p[0]), kn); for (k = -kn; k <= kn; k++) @@ -1266,23 +1266,31 @@ { static int v[4] = { 26, 1, 256, 231 }; - mpfr_init2 (p[i], i < 4 ? 5 + (randlimb () % 128) : 256); + mpfr_init2 (p[i], i < 4 ? 8 + (randlimb () % 128) : 256); if (i < 4) - mpfr_set_si_2exp (p[i], v[i], -5, MPFR_RNDN); + { + inexact = mpfr_set_si_2exp (p[i], v[i], -5, MPFR_RNDN); + MPFR_ASSERTN (inexact == 0); + } else { - mpfr_set_si_2exp (p[i], 1, 200, MPFR_RNDN); - mpfr_add (p[i], p[i], p[i-4], MPFR_RNDN); + inexact = mpfr_set_si_2exp (p[i], 1, 200, MPFR_RNDN); + MPFR_ASSERTN (inexact == 0); + inexact = mpfr_add (p[i], p[i], p[i-4], MPFR_RNDN); + MPFR_ASSERTN (inexact == 0); } ex[i] = mpfr_get_exp (p[i]) + 5; - MPFR_ASSERTD (ex[i] >= 0); + MPFR_ASSERTN (ex[i] >= 0); } mpfr_inits2 (3, p[8], p[9], p[10], (mpfr_ptr) 0); - mpfr_set_si_2exp (p[8], 1, 0, MPFR_RNDN); + inexact = mpfr_set_si_2exp (p[8], 1, 0, MPFR_RNDN); + MPFR_ASSERTN (inexact == 0); ex[8] = 5; - mpfr_set_si_2exp (p[9], 1, 0, MPFR_RNDN); /* will be epsilon */ + inexact = mpfr_set_si_2exp (p[9], 1, 0, MPFR_RNDN); /* will be epsilon */ + MPFR_ASSERTN (inexact == 0); ex[9] = 0; - mpfr_set_si_2exp (p[10], 7, 0, MPFR_RNDN); + inexact = mpfr_set_si_2exp (p[10], 7, 0, MPFR_RNDN); + MPFR_ASSERTN (inexact == 0); ex[10] = 5; for (i = 0; i < 11; i++) @@ -1294,9 +1302,9 @@ for (j = 0; j < numberof (e); j++) { inexact = mpfr_set_exp_t (ee, e[j], MPFR_RNDN); - MPFR_ASSERTD (inexact == 0); + MPFR_ASSERTN (inexact == 0); inexact = mpfr_get_z (MPFR_ZEXP (p[0]), ee, MPFR_RNDN); - MPFR_ASSERTD (inexact == 0); + MPFR_ASSERTN (inexact == 0); for (i = 1; i < 11; i++) mpz_set (MPFR_ZEXP (p[i]), MPFR_ZEXP (p[0])); for (i = 0; i < 11; i++) diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-06-02 11:26:59.031114881 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-06-02 11:26:59.175114040 +0000 @@ -0,0 +1 @@ +bernoulli-ziv diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-04-03 13:54:03.891945830 +0000 +++ mpfr-4.0.2-b/VERSION 2020-06-02 11:26:59.175114040 +0000 @@ -1 +1 @@ -4.0.2-p7 +4.0.2-p8 diff -Naurd mpfr-4.0.2-a/doc/README.dev mpfr-4.0.2-b/doc/README.dev --- mpfr-4.0.2-a/doc/README.dev 2019-01-07 14:06:05.000000000 +0000 +++ mpfr-4.0.2-b/doc/README.dev 2020-06-02 11:26:59.055114741 +0000 @@ -540,7 +540,9 @@ By default, a fixed seed is used. Only developers and testers should change the seed. -+ MPFR_CHECK_LARGEMEM: Define to enable expensive tests. ++ MPFR_CHECK_LARGEMEM: Define to enable tests that take a lot of memory. + ++ MPFR_CHECK_EXPENSIVE: Define to enable tests that take a lot of time. + MPFR_CHECK_LIBC_PRINTF: Define to enable comparisons with the printf diff -Naurd mpfr-4.0.2-a/src/bernoulli.c mpfr-4.0.2-b/src/bernoulli.c --- mpfr-4.0.2-a/src/bernoulli.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/bernoulli.c 2020-06-02 11:26:59.055114741 +0000 @@ -39,7 +39,7 @@ using Von Staudt–Clausen theorem, which says that the denominator of B[n] divides the product of all primes p such that p-1 divides n. Since B[n] = zeta(n) * 2*n!/(2pi)^n, we compute an approximation of - d * zeta(n) * 2*n!/(2pi)^n and round it to the nearest integer. */ + (2n+1)! * zeta(n) * 2*n!/(2pi)^n and round it to the nearest integer. */ static void mpfr_bernoulli_internal (mpz_t *b, unsigned long n) { @@ -86,8 +86,9 @@ mpfr_mul_ui (z, z, n, MPFR_RNDU); p = mpfr_get_ui (z, MPFR_RNDU); /* (n/e/2/pi)^n <= 2^p */ mpfr_clear (z); - /* the +14 term ensures no rounding failure up to n=10000 */ - prec += p + mpz_sizeinbase (den, 2) + 14; + prec += p + mpz_sizeinbase (den, 2); + /* the +2 term ensures no rounding failure up to n=10000 */ + prec += __gmpfr_ceil_log2 (prec) + 2; } try_again: @@ -184,7 +185,6 @@ mpz_mul_ui (t, t, n + 1); mpz_divexact (t, t, den); /* t was still n! */ mpz_mul (num, num, t); - mpz_set_ui (den, 1); mpfr_clear (y); mpfr_clear (z); diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-04-03 13:54:03.891945830 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-06-02 11:26:59.175114040 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p7" +#define MPFR_VERSION_STRING "4.0.2-p8" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-04-03 13:54:03.891945830 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-06-02 11:26:59.175114040 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p7"; + return "4.0.2-p8"; } diff -Naurd mpfr-4.0.2-a/tests/tgamma.c mpfr-4.0.2-b/tests/tgamma.c --- mpfr-4.0.2-a/tests/tgamma.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/tests/tgamma.c 2020-06-02 11:26:59.055114741 +0000 @@ -1055,6 +1055,48 @@ set_emax (emax); } +/* Bug reported by Frithjof Blomquist on May 19, 2020. + For the record, this bug was present since r8981 + (in mpfr_bernoulli_internal, den was wrongly reset to 1 in case + of failure in Ziv's loop). The bug only occurred up from r8986 + where the initial precision was reduced, but was potentially + present in any case of failure of Ziv's loop. */ +static void +bug20200519 (void) +{ + mpfr_prec_t prec = 25093; + mpfr_t x, y, z, d; + double dd; + size_t min_memory_limit, old_memory_limit; + + old_memory_limit = tests_memory_limit; + min_memory_limit = 24000000; + if (tests_memory_limit > 0 && tests_memory_limit < min_memory_limit) + tests_memory_limit = min_memory_limit; + + mpfr_init2 (x, prec); + mpfr_init2 (y, prec); + mpfr_init2 (z, prec + 100); + mpfr_init2 (d, 24); + mpfr_set_d (x, 2.5, MPFR_RNDN); + mpfr_gamma (y, x, MPFR_RNDN); + mpfr_gamma (z, x, MPFR_RNDN); + mpfr_sub (d, y, z, MPFR_RNDN); + mpfr_mul_2si (d, d, prec - mpfr_get_exp (y), MPFR_RNDN); + dd = mpfr_get_d (d, MPFR_RNDN); + if (dd < -0.5 || 0.5 < dd) + { + printf ("Error in bug20200519: dd=%f\n", dd); + exit (1); + } + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); + mpfr_clear (d); + + tests_memory_limit = old_memory_limit; +} + int main (int argc, char *argv[]) { @@ -1086,6 +1128,11 @@ data_check ("data/gamma", mpfr_gamma, "mpfr_gamma"); + /* this test takes about one minute */ + if (getenv ("MPFR_CHECK_EXPENSIVE") != NULL && + getenv ("MPFR_CHECK_LARGEMEM") != NULL) + bug20200519 (); + tests_end_mpfr (); return 0; } diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES --- mpfr-4.0.2-a/PATCHES 2020-06-02 12:39:10.963722927 +0000 +++ mpfr-4.0.2-b/PATCHES 2020-06-02 12:39:10.999722675 +0000 @@ -0,0 +1 @@ +fpif-dead-code diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION --- mpfr-4.0.2-a/VERSION 2020-06-02 11:26:59.175114040 +0000 +++ mpfr-4.0.2-b/VERSION 2020-06-02 12:39:10.999722675 +0000 @@ -1 +1 @@ -4.0.2-p8 +4.0.2-p9 diff -Naurd mpfr-4.0.2-a/src/fpif.c mpfr-4.0.2-b/src/fpif.c --- mpfr-4.0.2-a/src/fpif.c 2019-01-07 13:53:20.000000000 +0000 +++ mpfr-4.0.2-b/src/fpif.c 2020-06-02 12:39:10.987722759 +0000 @@ -90,11 +90,7 @@ { \ (buffer) = (unsigned char *) mpfr_reallocate_func \ ((buffer), *(buffer_size), (wanted_size)); \ - if ((buffer) == NULL) \ - { \ - *(buffer_size) = 0; \ - return NULL; \ - } \ + MPFR_ASSERTN((buffer) != 0); \ } \ *(buffer_size) = (wanted_size); \ } \ diff -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h --- mpfr-4.0.2-a/src/mpfr.h 2020-06-02 11:26:59.175114040 +0000 +++ mpfr-4.0.2-b/src/mpfr.h 2020-06-02 12:39:10.995722702 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 2 -#define MPFR_VERSION_STRING "4.0.2-p8" +#define MPFR_VERSION_STRING "4.0.2-p9" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c --- mpfr-4.0.2-a/src/version.c 2020-06-02 11:26:59.175114040 +0000 +++ mpfr-4.0.2-b/src/version.c 2020-06-02 12:39:10.999722675 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.2-p8"; + return "4.0.2-p9"; }