diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-04-27 12:40:46.838268769 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-04-27 12:40:46.870268475 +0000 @@ -0,0 +1 @@ +sub1sp1n-reuse diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-02-07 12:50:31.000000000 +0000 +++ mpfr-4.0.1-b/VERSION 2018-04-27 12:40:46.870268475 +0000 @@ -1 +1 @@ -4.0.1 +4.0.1-p1 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-02-07 12:50:31.000000000 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-04-27 12:40:46.870268475 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1" +#define MPFR_VERSION_STRING "4.0.1-p1" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/sub1sp.c mpfr-4.0.1-b/src/sub1sp.c --- mpfr-4.0.1-a/src/sub1sp.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/sub1sp.c 2018-04-27 12:40:46.858268585 +0000 @@ -375,13 +375,15 @@ } else /* cases (a), (c), (d) and (e) */ { - ap[0] = -MPFR_LIMB_ONE; /* rb=1 in case (e) and case (c) */ rb = d > GMP_NUMB_BITS + 1 || (d == GMP_NUMB_BITS + 1 && cp[0] == MPFR_LIMB_HIGHBIT); /* sb = 1 in case (d) and (e) */ sb = d > GMP_NUMB_BITS + 1 || (d == GMP_NUMB_BITS + 1 && cp[0] > MPFR_LIMB_HIGHBIT); + /* Warning: only set ap[0] last, otherwise in case ap=cp, + the above comparisons involving cp[0] would be wrong */ + ap[0] = -MPFR_LIMB_ONE; } } } diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-02-07 12:50:31.000000000 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-04-27 12:40:46.870268475 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1"; + return "4.0.1-p1"; } diff -Naurd mpfr-4.0.1-a/tests/tsub.c mpfr-4.0.1-b/tests/tsub.c --- mpfr-4.0.1-a/tests/tsub.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tsub.c 2018-04-27 12:40:46.858268585 +0000 @@ -22,12 +22,13 @@ #include "mpfr-test.h" -#ifdef CHECK_EXTERNAL static int test_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) { +#ifdef CHECK_EXTERNAL int res; int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c); + if (ok) { mpfr_print_raw (b); @@ -42,10 +43,69 @@ printf ("\n"); } return res; -} -#else -#define test_sub mpfr_sub +#else /* reuse test */ + int inex; + + inex = mpfr_sub (a, b, c, rnd_mode); + + if (a != b && a != c && ! MPFR_IS_NAN (a)) + { + mpfr_t t; + int reuse_b, reuse_c, inex_r; + + reuse_b = MPFR_PREC (a) == MPFR_PREC (b); + reuse_c = MPFR_PREC (a) == MPFR_PREC (c); + + if (reuse_b || reuse_c) + mpfr_init2 (t, MPFR_PREC (a)); + + if (reuse_b) + { + mpfr_set (t, b, MPFR_RNDN); + inex_r = mpfr_sub (t, t, c, rnd_mode); + if (!(mpfr_equal_p (t, a) && SAME_SIGN (inex_r, inex))) + { + printf ("reuse of b error in b - c in %s for\n", + mpfr_print_rnd_mode (rnd_mode)); + printf ("b = "); + mpfr_dump (b); + printf ("c = "); + mpfr_dump (c); + printf ("Expected "); mpfr_dump (a); + printf (" with inex = %d\n", inex); + printf ("Got "); mpfr_dump (t); + printf (" with inex = %d\n", inex_r); + exit (1); + } + } + + if (reuse_c) + { + mpfr_set (t, c, MPFR_RNDN); + inex_r = mpfr_sub (t, b, t, rnd_mode); + if (!(mpfr_equal_p (t, a) && SAME_SIGN (inex_r, inex))) + { + printf ("reuse of c error in b - c in %s for\n", + mpfr_print_rnd_mode (rnd_mode)); + printf ("b = "); + mpfr_dump (b); + printf ("c = "); + mpfr_dump (c); + printf ("Expected "); mpfr_dump (a); + printf (" with inex = %d\n", inex); + printf ("Got "); mpfr_dump (t); + printf (" with inex = %d\n", inex_r); + exit (1); + } + } + + if (reuse_b || reuse_c) + mpfr_clear (t); + } + + return inex; #endif +} static void check_diverse (void) @@ -940,6 +1000,80 @@ } } +/* Fails with r12281: "reuse of c error in b - c in MPFR_RNDN". */ +static void +bug20180217 (void) +{ + mpfr_t x, y, z1, z2; + int r, p, d, i, inex1, inex2; + + for (p = 3; p <= 3 + 4 * GMP_NUMB_BITS; p++) + { + mpfr_inits2 (p, x, y, z1, z2, (mpfr_ptr) 0); + for (d = p; d <= p+4; d++) + { + mpfr_set_ui (x, 1, MPFR_RNDN); + mpfr_set_ui_2exp (y, 1, -d, MPFR_RNDN); + for (i = 0; i < 3; i++) + { + RND_LOOP_NO_RNDF (r) + { + mpfr_set (z1, x, MPFR_RNDN); + if (d == p) + { + mpfr_nextbelow (z1); + if (i == 0) + inex1 = 0; + else if (r == MPFR_RNDD || r == MPFR_RNDZ || + (r == MPFR_RNDN && i > 1)) + { + mpfr_nextbelow (z1); + inex1 = -1; + } + else + inex1 = 1; + } + else if (r == MPFR_RNDD || r == MPFR_RNDZ || + (r == MPFR_RNDN && d == p+1 && i > 0)) + { + mpfr_nextbelow (z1); + inex1 = -1; + } + else + inex1 = 1; + inex2 = test_sub (z2, x, y, (mpfr_rnd_t) r); + if (!(mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2))) + { + printf ("Error in bug20180217 with " + "p=%d, d=%d, i=%d, %s\n", p, d, i, + mpfr_print_rnd_mode ((mpfr_rnd_t) r)); + printf ("x = "); + mpfr_dump (x); + printf ("y = "); + mpfr_dump (y); + printf ("Expected "); mpfr_dump (z1); + printf (" with inex = %d\n", inex1); + printf ("Got "); mpfr_dump (z2); + printf (" with inex = %d\n", inex2); + exit (1); + } + } + if (i == 0) + mpfr_nextabove (y); + else + { + if (p < 6) + break; + mpfr_nextbelow (y); + mpfr_mul_ui (y, y, 25, MPFR_RNDD); + mpfr_div_2ui (y, y, 4, MPFR_RNDN); + } + } + } + mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0); + } +} + #define TEST_FUNCTION test_sub #define TWO_ARGS #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS) @@ -962,6 +1096,7 @@ check_inexact (); check_max_almosteven (); bug_ddefour (); + bug20180217 (); for (p=2; p<200; p++) for (i=0; i<50; i++) check_two_sum (p); diff -Naurd mpfr-4.0.1-a/tests/tsub1sp.c mpfr-4.0.1-b/tests/tsub1sp.c --- mpfr-4.0.1-a/tests/tsub1sp.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tsub1sp.c 2018-04-27 12:40:46.858268585 +0000 @@ -284,6 +284,91 @@ } } +static void +coverage (void) +{ + mpfr_t a, b, c, d, u; + int inex; + + /* coverage test in mpfr_sub1sp: case d=1, limb > MPFR_LIMB_HIGHBIT, RNDF + and also RNDZ */ + mpfr_init2 (a, 3 * GMP_NUMB_BITS); + mpfr_init2 (b, 3 * GMP_NUMB_BITS); + mpfr_init2 (c, 3 * GMP_NUMB_BITS); + mpfr_init2 (d, 3 * GMP_NUMB_BITS); + mpfr_init2 (u, 3 * GMP_NUMB_BITS); + mpfr_set_ui (b, 1, MPFR_RNDN); + mpfr_nextbelow (b); /* b = 1 - 2^(-p) */ + mpfr_set_prec (c, GMP_NUMB_BITS); + mpfr_set_ui_2exp (c, 1, -1, MPFR_RNDN); + mpfr_nextbelow (c); + mpfr_nextbelow (c); /* c = 1/2 - 2*2^(-GMP_NUMB_BITS-1) */ + mpfr_prec_round (c, 3 * GMP_NUMB_BITS, MPFR_RNDN); + mpfr_nextbelow (c); /* c = 1/2 - 2*2^(-GMP_NUMB_BITS-1) - 2^(-p-1) */ + /* b-c = c */ + mpfr_sub (a, b, c, MPFR_RNDF); + mpfr_sub (d, b, c, MPFR_RNDD); + mpfr_sub (u, b, c, MPFR_RNDU); + /* check a = d or u */ + MPFR_ASSERTN(mpfr_equal_p (a, d) || mpfr_equal_p (a, u)); + + /* coverage test in mpfr_sub1sp: case d=p, RNDN, sb = 0, significand of b + is even but b<>2^e, (case 1e) */ + mpfr_set_prec (a, 3 * GMP_NUMB_BITS); + mpfr_set_prec (b, 3 * GMP_NUMB_BITS); + mpfr_set_prec (c, 3 * GMP_NUMB_BITS); + mpfr_set_ui (b, 1, MPFR_RNDN); + mpfr_nextabove (b); + mpfr_nextabove (b); + mpfr_set_ui_2exp (c, 1, -3 * GMP_NUMB_BITS, MPFR_RNDN); + inex = mpfr_sub (a, b, c, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_equal_p (a, b)); + + mpfr_clear (a); + mpfr_clear (b); + mpfr_clear (c); + mpfr_clear (d); + mpfr_clear (u); +} + +/* bug in mpfr_sub1sp1n, made generic */ +static void +bug20180217 (mpfr_prec_t pmax) +{ + mpfr_t a, b, c; + int inex; + mpfr_prec_t p; + + for (p = MPFR_PREC_MIN; p <= pmax; p++) + { + mpfr_init2 (a, p); + mpfr_init2 (b, p); + mpfr_init2 (c, p); + mpfr_set_ui (b, 1, MPFR_RNDN); /* b = 1 */ + mpfr_set_ui_2exp (c, 1, -p-1, MPFR_RNDN); /* c = 2^(-p-1) */ + /* a - b = 1 - 2^(-p-1) and should be rounded to 1 (case 2f of + mpfr_sub1sp) */ + inex = mpfr_sub (a, b, c, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); + /* check also when a=b */ + mpfr_set_ui (a, 1, MPFR_RNDN); + inex = mpfr_sub (a, a, c, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); + /* and when a=c */ + mpfr_set_ui (b, 1, MPFR_RNDN); /* b = 1 */ + mpfr_set_ui_2exp (a, 1, -p-1, MPFR_RNDN); + inex = mpfr_sub (a, b, a, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0); + mpfr_clear (a); + mpfr_clear (b); + mpfr_clear (c); + } +} + int main (void) { @@ -291,6 +376,8 @@ tests_start_mpfr (); + bug20180217 (1024); + coverage (); compare_sub_sub1sp (); test20170208 (); bug20170109 (); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-04-27 12:45:53.139452673 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-04-27 12:45:53.171452379 +0000 @@ -0,0 +1 @@ +fma diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-04-27 12:40:46.870268475 +0000 +++ mpfr-4.0.1-b/VERSION 2018-04-27 12:45:53.171452379 +0000 @@ -1 +1 @@ -4.0.1-p1 +4.0.1-p2 diff -Naurd mpfr-4.0.1-a/src/fma.c mpfr-4.0.1-b/src/fma.c --- mpfr-4.0.1-a/src/fma.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/fma.c 2018-04-27 12:45:53.159452489 +0000 @@ -225,194 +225,73 @@ if (MPFR_IS_INF (u)) /* overflow */ { + int sign_u = MPFR_SIGN (u); + MPFR_LOG_MSG (("Overflow on x*y\n", 0)); + MPFR_GROUP_CLEAR (group); /* we no longer need u */ /* Let's eliminate the obvious case where x*y and z have the same sign. No possible cancellation -> real overflow. Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3, - then |x*y| >= 2^(emax+1), and |x*y + z| >= 2^emax. This case + then |x*y| >= 2^(emax+1), and |x*y + z| > 2^emax. This case is also an overflow. */ - if (MPFR_SIGN (u) == MPFR_SIGN (z) || e >= __gmpfr_emax + 3) + if (sign_u == MPFR_SIGN (z) || e >= __gmpfr_emax + 3) { - MPFR_GROUP_CLEAR (group); MPFR_SAVE_EXPO_FREE (expo); - return mpfr_overflow (s, rnd_mode, MPFR_SIGN (z)); + return mpfr_overflow (s, rnd_mode, sign_u); } - - /* E(x) + E(y) <= emax+2, therefore |x*y| < 2^(emax+2), and - (x/4)*y does not overflow (let's recall that the result - is exact with an unbounded exponent range). It does not - underflow either, because x*y overflows and the exponent - range is large enough. */ - inexact = mpfr_div_2ui (u, x, 2, MPFR_RNDN); - MPFR_ASSERTN (inexact == 0); - inexact = mpfr_mul (u, u, y, MPFR_RNDN); - MPFR_ASSERTN (inexact == 0); - - /* Now, we need to add z/4... But it may underflow! */ - { - mpfr_t zo4; - mpfr_srcptr zz; - MPFR_BLOCK_DECL (flags); - - if (MPFR_GET_EXP (u) > MPFR_GET_EXP (z) && - MPFR_GET_EXP (u) - MPFR_GET_EXP (z) > MPFR_PREC (u)) - { - /* |z| < ulp(u)/2, therefore one can use z instead of z/4. */ - zz = z; - } - else - { - mpfr_init2 (zo4, MPFR_PREC (z)); - if (mpfr_div_2ui (zo4, z, 2, MPFR_RNDZ)) - { - /* The division by 4 underflowed! */ - MPFR_ASSERTN (0); /* TODO... */ - } - zz = zo4; - } - - /* Let's recall that u = x*y/4 and zz = z/4 (or z if the - following addition would give the same result). */ - MPFR_BLOCK (flags, inexact = mpfr_add (s, u, zz, rnd_mode)); - /* u and zz have different signs, so that an overflow - is not possible. But an underflow is theoretically - possible! */ - if (MPFR_UNDERFLOW (flags)) - { - MPFR_ASSERTN (zz != z); - MPFR_ASSERTN (0); /* TODO... */ - mpfr_clears (zo4, u, (mpfr_ptr) 0); - } - else - { - int inex2; - - if (zz != z) - mpfr_clear (zo4); - MPFR_GROUP_CLEAR (group); - MPFR_ASSERTN (! MPFR_OVERFLOW (flags)); - inex2 = mpfr_mul_2ui (s, s, 2, rnd_mode); - if (inex2) /* overflow */ - { - inexact = inex2; - MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); - } - goto end; - } - } } - else /* underflow: one has |xy| < 2^(emin-1). */ + else /* underflow: one has |x*y| < 2^(emin-1). */ { - unsigned long scale = 0; - mpfr_t scaled_z; - mpfr_srcptr new_z; - mpfr_exp_t diffexp; - mpfr_prec_t pzs; - int xy_underflows; - MPFR_LOG_MSG (("Underflow on x*y\n", 0)); - /* Let's scale z so that ulp(z) > 2^emin and ulp(s) > 2^emin - (the + 1 on MPFR_PREC (s) is necessary because the exponent - of the result can be EXP(z) - 1). */ - diffexp = MPFR_GET_EXP (z) - __gmpfr_emin; - pzs = MAX (MPFR_PREC (z), MPFR_PREC (s) + 1); - MPFR_LOG_MSG (("diffexp=%" MPFR_EXP_FSPEC "d pzs=%Pd\n", - diffexp, pzs)); - if (diffexp <= pzs) - { - mpfr_uexp_t uscale; - mpfr_t scaled_v; - MPFR_BLOCK_DECL (flags); - - uscale = (mpfr_uexp_t) pzs - diffexp + 1; - MPFR_ASSERTN (uscale > 0); - MPFR_ASSERTN (uscale <= ULONG_MAX); - scale = uscale; - mpfr_init2 (scaled_z, MPFR_PREC (z)); - inexact = mpfr_mul_2ui (scaled_z, z, scale, MPFR_RNDN); - MPFR_ASSERTN (inexact == 0); /* TODO: overflow case */ - new_z = scaled_z; - /* Now we need to recompute u = xy * 2^scale. */ - MPFR_BLOCK (flags, - if (MPFR_GET_EXP (x) < MPFR_GET_EXP (y)) - { - mpfr_init2 (scaled_v, precx); - mpfr_mul_2ui (scaled_v, x, scale, MPFR_RNDN); - mpfr_mul (u, scaled_v, y, MPFR_RNDN); - } - else - { - mpfr_init2 (scaled_v, precy); - mpfr_mul_2ui (scaled_v, y, scale, MPFR_RNDN); - mpfr_mul (u, x, scaled_v, MPFR_RNDN); - }); - mpfr_clear (scaled_v); - MPFR_ASSERTN (! MPFR_OVERFLOW (flags)); - xy_underflows = MPFR_UNDERFLOW (flags); - } - else - { - new_z = z; - xy_underflows = 1; - } - - MPFR_LOG_MSG (("scale=%lu xy_underflows=%d\n", - scale, xy_underflows)); - - if (xy_underflows) + /* Easy cases: when 2^(emin-1) <= 1/2 * min(ulp(z),ulp(s)), + one can replace x*y by sign(x*y) * 2^(emin-1). Note that + this is even true in case of equality for MPFR_RNDN thanks + to the even-rounding rule. + The + 1 on MPFR_PREC (s) is necessary because the exponent + of the result can be EXP(z) - 1. */ + if (MPFR_GET_EXP (z) - __gmpfr_emin >= + MAX (MPFR_PREC (z), MPFR_PREC (s) + 1)) { - /* Let's replace xy by sign(xy) * 2^(emin-1). */ MPFR_PREC (u) = MPFR_PREC_MIN; mpfr_setmin (u, __gmpfr_emin); MPFR_SET_SIGN (u, MPFR_MULT_SIGN (MPFR_SIGN (x), MPFR_SIGN (y))); + mpfr_clear_flags (); + goto add; } - { - MPFR_BLOCK_DECL (flags); + MPFR_GROUP_CLEAR (group); /* we no longer need u */ + } - MPFR_BLOCK (flags, inexact = mpfr_add (s, u, new_z, rnd_mode)); - MPFR_LOG_MSG (("inexact=%d\n", inexact)); - MPFR_GROUP_CLEAR (group); - if (scale != 0) - { - int inex2; + /* Let's use UBF to resolve the overflow/underflow issues. */ + { + mpfr_ubf_t uu; + mp_size_t un; + mpfr_limb_ptr up; + MPFR_TMP_DECL(marker); - mpfr_clear (scaled_z); - /* Here an overflow is theoretically possible, in which case - the result may be wrong, hence the assert. An underflow - is not possible, but let's check that anyway. */ - MPFR_ASSERTN (! MPFR_OVERFLOW (flags)); /* TODO... */ - MPFR_ASSERTN (! MPFR_UNDERFLOW (flags)); /* not possible */ - if (rnd_mode == MPFR_RNDN && - MPFR_GET_EXP (s) == __gmpfr_emin - 1 + scale && - mpfr_powerof2_raw (s)) - { - MPFR_LOG_MSG (("Double rounding\n", 0)); - rnd_mode = (MPFR_IS_NEG (s) ? inexact <= 0 : inexact >= 0) - ? MPFR_RNDZ : MPFR_RNDA; - } - inex2 = mpfr_div_2ui (s, s, scale, rnd_mode); - MPFR_LOG_MSG (("inex2=%d\n", inex2)); - if (inex2) /* underflow */ - inexact = inex2; - } - } + MPFR_LOG_MSG (("Use UBF\n", 0)); - /* FIXME/TODO: I'm not sure that the following is correct. - Check for possible spurious exceptions due to intermediate - computations. */ - MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); - goto end; - } + MPFR_TMP_MARK (marker); + un = MPFR_LIMB_SIZE (x) + MPFR_LIMB_SIZE (y); + MPFR_TMP_INIT (up, uu, (mpfr_prec_t) un * GMP_NUMB_BITS, un); + mpfr_ubf_mul_exact (uu, x, y); + mpfr_clear_flags (); + inexact = mpfr_add (s, (mpfr_srcptr) uu, z, rnd_mode); + MPFR_UBF_CLEAR_EXP (uu); + MPFR_TMP_FREE (marker); + } + } + else + { + add: + inexact = mpfr_add (s, u, z, rnd_mode); + MPFR_GROUP_CLEAR (group); } - inexact = mpfr_add (s, u, z, rnd_mode); - MPFR_GROUP_CLEAR (group); MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); - end: MPFR_SAVE_EXPO_FREE (expo); return mpfr_check_range (s, inexact, rnd_mode); } diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-04-27 12:40:46.870268475 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-04-27 12:45:53.171452379 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p1" +#define MPFR_VERSION_STRING "4.0.1-p2" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-04-27 12:40:46.870268475 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-04-27 12:45:53.171452379 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p1"; + return "4.0.1-p2"; } diff -Naurd mpfr-4.0.1-a/tests/tfma.c mpfr-4.0.1-b/tests/tfma.c --- mpfr-4.0.1-a/tests/tfma.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tfma.c 2018-04-27 12:45:53.163452452 +0000 @@ -196,6 +196,238 @@ } static void +test_overflow3 (void) +{ + mpfr_t x, y, z, r; + int inex; + mpfr_prec_t p = 8; + mpfr_flags_t ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT, flags; + int i, j, k; + unsigned int neg; + + mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0); + for (i = 0; i < 2; i++) + { + mpfr_init2 (r, 2 * p + i); + mpfr_set_ui_2exp (x, 1, mpfr_get_emax () - 1, MPFR_RNDN); + mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */ + for (j = 1; j <= 2; j++) + for (k = 0; k <= 1; k++) + { + mpfr_set_si_2exp (z, -1, mpfr_get_emax () - mpfr_get_prec (r) - j, + MPFR_RNDN); + if (k) + mpfr_nextabove (z); + for (neg = 0; neg <= 3; neg++) + { + mpfr_clear_flags (); + /* (The following applies for neg = 0 or 2, all the signs + need to be reversed for neg = 1 or 3.) + We have x*y = 2^emax and + z = - 2^(emax-2p-i-j) * (1-k*2^(-p)), thus + x*y+z = 2^emax - 2^(emax-2p-i-j) + k*2^(emax-3p-i-j) + should overflow. Indeed it is >= the midpoint of + 2^emax - 2^(emax-2p-i) and 2^emax, the midpoint + being obtained for j = 1 and k = 0. */ + inex = mpfr_fma (r, x, y, z, MPFR_RNDN); + flags = __gmpfr_flags; + if (! mpfr_inf_p (r) || flags != ex_flags || + ((neg & 1) == 0 ? + (inex <= 0 || MPFR_IS_NEG (r)) : + (inex >= 0 || MPFR_IS_POS (r)))) + { + printf ("Error in test_overflow3 for " + "i=%d j=%d k=%d neg=%u\n", i, j, k, neg); + printf ("Expected %c@Inf@\n with inex %c 0 and flags:", + (neg & 1) == 0 ? '+' : '-', + (neg & 1) == 0 ? '>' : '<'); + flags_out (ex_flags); + printf ("Got "); + mpfr_dump (r); + printf (" with inex = %d and flags:", inex); + flags_out (flags); + exit (1); + } + if (neg == 0 || neg == 2) + mpfr_neg (x, x, MPFR_RNDN); + if (neg == 1 || neg == 3) + mpfr_neg (y, y, MPFR_RNDN); + mpfr_neg (z, z, MPFR_RNDN); + } /* neg */ + } /* k */ + mpfr_clear (r); + } /* i */ + mpfr_clears (x, y, z, (mpfr_ptr) 0); +} + +static void +test_overflow4 (void) +{ + mpfr_t x, y, z, r1, r2; + mpfr_exp_t emax, e; + mpfr_prec_t px; + mpfr_flags_t flags1, flags2; + int inex1, inex2; + int ei, i, j; + int below; + unsigned int neg; + + emax = mpfr_get_emax (); + + mpfr_init2 (y, MPFR_PREC_MIN); + mpfr_set_ui (y, 2, MPFR_RNDN); /* y = 2 */ + + mpfr_init2 (z, 8); + + for (px = 17; px < 256; px *= 2) + { + mpfr_init2 (x, px); + mpfr_inits2 (px - 8, r1, r2, (mpfr_ptr) 0); + for (ei = 0; ei <= 1; ei++) + { + e = ei ? emax : 0; + mpfr_set_ui_2exp (x, 1, e - 1, MPFR_RNDN); + mpfr_nextabove (x); /* x = 2^(e - 1) + 2^(e - px) */ + /* x*y = 2^e + 2^(e - px + 1), which internally overflows + when e = emax. */ + for (i = -4; i <= 4; i++) + for (j = 2; j <= 3; j++) + { + mpfr_set_si_2exp (z, -j, e - px + i, MPFR_RNDN); + /* If |z| <= 2^(e - px + 1), then x*y + z >= 2^e and + RZ(x*y + z) = 2^e with an unbounded exponent range. + If |z| > 2^(e - px + 1), then RZ(x*y + z) is the + predecessor of 2^e (since |z| < ulp(r)/2); this + occurs when i > 0 and when i = 0 and j > 2 */ + mpfr_set_ui_2exp (r1, 1, e - 1, MPFR_RNDN); + below = i > 0 || (i == 0 && j > 2); + if (below) + mpfr_nextbelow (r1); + mpfr_clear_flags (); + inex1 = mpfr_mul_2ui (r1, r1, 1, MPFR_RNDZ); + if (below || e < emax) + { + inex1 = i == 0 && j == 2 ? 0 : -1; + flags1 = inex1 ? MPFR_FLAGS_INEXACT : 0; + } + else + { + MPFR_ASSERTN (inex1 < 0); + flags1 = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW; + MPFR_ASSERTN (flags1 == __gmpfr_flags); + } + for (neg = 0; neg <= 3; neg++) + { + mpfr_clear_flags (); + inex2 = mpfr_fma (r2, x, y, z, MPFR_RNDZ); + flags2 = __gmpfr_flags; + if (! (mpfr_equal_p (r1, r2) && + SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in test_overflow4 for " + "px=%d ei=%d i=%d j=%d neg=%u\n", + (int) px, ei, i, j, neg); + printf ("Expected "); + mpfr_dump (r1); + printf ("with inex = %d and flags:", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (r2); + printf ("with inex = %d and flags:", inex2); + flags_out (flags2); + exit (1); + } + if (neg == 0 || neg == 2) + mpfr_neg (x, x, MPFR_RNDN); + if (neg == 1 || neg == 3) + mpfr_neg (y, y, MPFR_RNDN); + mpfr_neg (z, z, MPFR_RNDN); + mpfr_neg (r1, r1, MPFR_RNDN); + inex1 = - inex1; + } + } + } + mpfr_clears (x, r1, r2, (mpfr_ptr) 0); + } + + mpfr_clears (y, z, (mpfr_ptr) 0); +} + +static void +test_overflow5 (void) +{ + mpfr_t x, y, z, r1, r2; + mpfr_exp_t emax; + int inex1, inex2; + int i, rnd; + unsigned int neg, negr; + + emax = mpfr_get_emax (); + + mpfr_init2 (x, 123); + mpfr_init2 (y, 45); + mpfr_init2 (z, 67); + mpfr_inits2 (89, r1, r2, (mpfr_ptr) 0); + + mpfr_set_ui_2exp (x, 1, emax - 1, MPFR_RNDN); + + for (i = 3; i <= 17; i++) + { + mpfr_set_ui (y, i, MPFR_RNDN); + mpfr_set_ui_2exp (z, 1, emax - 1, MPFR_RNDN); + for (neg = 0; neg < 8; neg++) + { + mpfr_setsign (x, x, neg & 1, MPFR_RNDN); + mpfr_setsign (y, y, neg & 2, MPFR_RNDN); + mpfr_setsign (z, z, neg & 4, MPFR_RNDN); + + /* |x*y + z| = (i +/- 1) * 2^(emax - 1) >= 2^emax (overflow) + and x*y + z has the same sign as x*y. */ + negr = (neg ^ (neg >> 1)) & 1; + + RND_LOOP (rnd) + { + mpfr_set_inf (r1, 1); + if (MPFR_IS_LIKE_RNDZ ((mpfr_rnd_t) rnd, negr)) + { + mpfr_nextbelow (r1); + inex1 = -1; + } + else + inex1 = 1; + + if (negr) + { + mpfr_neg (r1, r1, MPFR_RNDN); + inex1 = - inex1; + } + + mpfr_clear_flags (); + inex2 = mpfr_fma (r2, x, y, z, (mpfr_rnd_t) rnd); + MPFR_ASSERTN (__gmpfr_flags == + (MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW)); + + if (! (mpfr_equal_p (r1, r2) && SAME_SIGN (inex1, inex2))) + { + printf ("Error in test_overflow5 for i=%d neg=%u %s\n", + i, neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("Expected "); + mpfr_dump (r1); + printf ("with inex = %d\n", inex1); + printf ("Got "); + mpfr_dump (r2); + printf ("with inex = %d\n", inex2); + exit (1); + } + } /* rnd */ + } /* neg */ + } /* i */ + + mpfr_clears (x, y, z, r1, r2, (mpfr_ptr) 0); +} + +static void test_underflow1 (void) { mpfr_t x, y, z, r; @@ -281,59 +513,128 @@ static void test_underflow2 (void) { - mpfr_t x, y, z, r; - int b, i, inex, same, err = 0; + mpfr_t x, y, z, r1, r2; + int e, b, i, prec = 32, pz, inex; + unsigned int neg; - mpfr_inits2 (32, x, y, z, r, (mpfr_ptr) 0); + mpfr_init2 (x, MPFR_PREC_MIN); + mpfr_inits2 (prec, y, z, r1, r2, (mpfr_ptr) 0); - mpfr_set_si_2exp (z, 1, mpfr_get_emin (), MPFR_RNDN); /* z = 2^emin */ - mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN); /* x = 2^emin */ + mpfr_set_si_2exp (x, 1, mpfr_get_emin () - 1, MPFR_RNDN); + /* x = 2^(emin-1) */ - for (b = 0; b <= 1; b++) + for (e = -1; e <= prec + 2; e++) { - for (i = 15; i <= 17; i++) + mpfr_set (z, x, MPFR_RNDN); + /* z = x = 2^(emin+e) */ + for (b = 0; b <= 1; b++) { - mpfr_set_si_2exp (y, i, -4 - MPFR_PREC (z), MPFR_RNDN); - /* z = 1.000...00b - * xy = 01111 - * or 10000 - * or 10001 - */ - mpfr_clear_flags (); - inex = mpfr_fma (r, x, y, z, MPFR_RNDN); -#define STRTU2 "Error in test_underflow2 (b = %d, i = %d)\n " - if (__gmpfr_flags != MPFR_FLAGS_INEXACT) - { - printf (STRTU2 "flags = %u instead of %u\n", b, i, - (unsigned int) __gmpfr_flags, - (unsigned int) MPFR_FLAGS_INEXACT); - err = 1; - } - same = i == 15 || (i == 16 && b == 0); - if (same ? (inex >= 0) : (inex <= 0)) - { - printf (STRTU2 "incorrect ternary value (%d instead of %c 0)\n", - b, i, inex, same ? '<' : '>'); - err = 1; - } - mpfr_set (y, z, MPFR_RNDN); - if (!same) - mpfr_nextabove (y); - if (! mpfr_equal_p (r, y)) + for (pz = prec - 4 * (b == 0); pz <= prec + 4; pz++) { - printf (STRTU2 "expected ", b, i); - mpfr_dump (y); - printf (" got "); - mpfr_dump (r); - err = 1; - } - } - mpfr_nextabove (z); - } + inex = mpfr_prec_round (z, pz, MPFR_RNDZ); + MPFR_ASSERTN (inex == 0); + for (i = 15; i <= 17; i++) + { + mpfr_flags_t flags1, flags2; + int inex1, inex2; - if (err) - exit (1); - mpfr_clears (x, y, z, r, (mpfr_ptr) 0); + mpfr_set_si_2exp (y, i, -4 - prec, MPFR_RNDN); + + /* <--- r ---> + * z = 1.000...00b with b = 0 or 1 + * xy = 01111 (i = 15) + * or 10000 (i = 16) + * or 10001 (i = 17) + * + * x, y, and z will be modified to test the different sign + * combinations. In the case b = 0 (i.e. |z| is a power of + * 2) and x*y has a different sign from z, then y will be + * divided by 2, so that i = 16 corresponds to a midpoint. + */ + + for (neg = 0; neg < 8; neg++) + { + int xyneg, prev_binade; + + mpfr_setsign (x, x, neg & 1, MPFR_RNDN); + mpfr_setsign (y, y, neg & 2, MPFR_RNDN); + mpfr_setsign (z, z, neg & 4, MPFR_RNDN); + + xyneg = (neg ^ (neg >> 1)) & 1; /* true iff x*y < 0 */ + + /* If a change of binade occurs by adding x*y to z + exactly, then take into account the fact that the + midpoint has a different exponent. */ + prev_binade = b == 0 && (xyneg ^ MPFR_IS_NEG (z)); + if (prev_binade) + mpfr_div_2ui (y, y, 1, MPFR_RNDN); + + mpfr_set (r1, z, MPFR_RNDN); + flags1 = MPFR_FLAGS_INEXACT; + + if (e == -1 && i == 17 && b == 0 && + (xyneg ^ (neg >> 2)) != 0) + { + /* Special underflow case. */ + flags1 |= MPFR_FLAGS_UNDERFLOW; + inex1 = xyneg ? 1 : -1; + } + else if (i == 15 || (i == 16 && b == 0)) + { + /* round toward z */ + inex1 = xyneg ? 1 : -1; + } + else if (xyneg) + { + /* round away from z, with x*y < 0 */ + mpfr_nextbelow (r1); + inex1 = -1; + } + else + { + /* round away from z, with x*y > 0 */ + mpfr_nextabove (r1); + inex1 = 1; + } + + mpfr_clear_flags (); + inex2 = mpfr_fma (r2, x, y, z, MPFR_RNDN); + flags2 = __gmpfr_flags; + + if (! (mpfr_equal_p (r1, r2) && + SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in test_underflow2 for " + "e=%d b=%d pz=%d i=%d neg=%u\n", + e, b, pz, i, neg); + printf ("Expected "); + mpfr_dump (r1); + printf ("with inex = %d and flags:", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (r2); + printf ("with inex = %d and flags:", inex2); + flags_out (flags2); + exit (1); + } + + /* Restore y. */ + if (prev_binade) + mpfr_mul_2ui (y, y, 1, MPFR_RNDN); + } /* neg */ + } /* i */ + } /* pz */ + + inex = mpfr_prec_round (z, prec, MPFR_RNDZ); + MPFR_ASSERTN (inex == 0); + MPFR_SET_POS (z); + mpfr_nextabove (z); + } /* b */ + mpfr_mul_2ui (x, x, 1, MPFR_RNDN); + } /* e */ + + mpfr_clears (x, y, z, r1, r2, (mpfr_ptr) 0); } static void @@ -397,6 +698,185 @@ mpfr_clears (x, y, z, t1, t2, (mpfr_ptr) 0); } +/* Test s = x*y + z with PREC(z) > PREC(s) + 1, x*y underflows, where + z + x*y and z + sign(x*y) * 2^(emin-1) do not give the same result. + x = 2^emin + y = 2^(-8) + z = 2^emin * (2^PREC(s) + k - 2^(-1)) + with k = 3 for MPFR_RNDN and k = 2 for the directed rounding modes. + Also test the opposite versions with neg != 0. +*/ +static void +test_underflow4 (void) +{ + mpfr_t x, y, z, s1, s2; + mpfr_prec_t ps = 32; + int inex, rnd; + + mpfr_inits2 (MPFR_PREC_MIN, x, y, (mpfr_ptr) 0); + mpfr_inits2 (ps, s1, s2, (mpfr_ptr) 0); + mpfr_init2 (z, ps + 2); + + inex = mpfr_set_si_2exp (x, 1, mpfr_get_emin (), MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_set_si_2exp (y, 1, -8, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + RND_LOOP_NO_RNDF (rnd) + { + mpfr_flags_t flags1, flags2; + int inex1, inex2; + unsigned int neg; + + inex = mpfr_set_si_2exp (z, 1 << 1, ps, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_div_2ui (z, z, 1, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_add_ui (z, z, rnd == MPFR_RNDN ? 3 : 2, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_mul (z, z, x, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + for (neg = 0; neg <= 3; neg++) + { + inex1 = mpfr_set (s1, z, (mpfr_rnd_t) rnd); + flags1 = MPFR_FLAGS_INEXACT; + + mpfr_clear_flags (); + inex2 = mpfr_fma (s2, x, y, z, (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (! (mpfr_equal_p (s1, s2) && + SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in test_underflow4 for neg=%u %s\n", + neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("Expected "); + mpfr_dump (s1); + printf (" with inex ~ %d, flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (s2); + printf (" with inex ~ %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + + if (neg == 0 || neg == 2) + mpfr_neg (x, x, MPFR_RNDN); + if (neg == 1 || neg == 3) + mpfr_neg (y, y, MPFR_RNDN); + mpfr_neg (z, z, MPFR_RNDN); + } + } + + mpfr_clears (x, y, z, s1, s2, (mpfr_ptr) 0); +} + +/* Test s = x*y + z on: + x = +/- 2^emin + y = +/- 2^(-3) + z = +/- 2^(emin + PREC(s)) and MPFR numbers close to this value. + with PREC(z) from PREC(s) - 2 to PREC(s) + 8. +*/ +static void +test_underflow5 (void) +{ + mpfr_t w, x, y, z, s1, s2, t; + mpfr_exp_t emin; + mpfr_prec_t ps = 32; + int inex, i, j, rnd; + unsigned int neg; + + mpfr_inits2 (MPFR_PREC_MIN, w, x, y, (mpfr_ptr) 0); + mpfr_inits2 (ps, s1, s2, (mpfr_ptr) 0); + mpfr_init2 (t, ps + 9); + + emin = mpfr_get_emin (); + + inex = mpfr_set_si_2exp (w, 1, emin, MPFR_RNDN); /* w = 2^emin */ + MPFR_ASSERTN (inex == 0); + inex = mpfr_set_si (x, 1, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_set_si_2exp (y, 1, -3, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + for (i = -2; i <= 8; i++) + { + mpfr_init2 (z, ps + i); + inex = mpfr_set_si_2exp (z, 1, ps, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + for (j = 1; j <= 32; j++) + mpfr_nextbelow (z); + + for (j = -32; j <= 32; j++) + { + for (neg = 0; neg < 8; neg++) + { + mpfr_setsign (x, x, neg & 1, MPFR_RNDN); + mpfr_setsign (y, y, neg & 2, MPFR_RNDN); + mpfr_setsign (z, z, neg & 4, MPFR_RNDN); + + inex = mpfr_fma (t, x, y, z, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + inex = mpfr_mul (x, x, w, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_mul (z, z, w, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + + RND_LOOP_NO_RNDF (rnd) + { + mpfr_flags_t flags1, flags2; + int inex1, inex2; + + mpfr_clear_flags (); + inex1 = mpfr_mul (s1, w, t, (mpfr_rnd_t) rnd); + flags1 = __gmpfr_flags; + + mpfr_clear_flags (); + inex2 = mpfr_fma (s2, x, y, z, (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (! (mpfr_equal_p (s1, s2) && + SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in test_underflow5 on " + "i=%d j=%d neg=%u %s\n", i, j, neg, + mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + printf ("Expected "); + mpfr_dump (s1); + printf (" with inex ~ %d, flags =", inex1); + flags_out (flags1); + printf ("Got "); + mpfr_dump (s2); + printf (" with inex ~ %d, flags =", inex2); + flags_out (flags2); + exit (1); + } + } /* rnd */ + + inex = mpfr_div (x, x, w, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + inex = mpfr_div (z, z, w, MPFR_RNDN); + MPFR_ASSERTN (inex == 0); + } /* neg */ + + MPFR_SET_POS (z); /* restore the value before the loop on neg */ + mpfr_nextabove (z); + } /* j */ + + mpfr_clear (z); + } /* i */ + + mpfr_clears (w, x, y, s1, s2, t, (mpfr_ptr) 0); +} + static void bug20101018 (void) { @@ -533,6 +1013,7 @@ { mpfr_t x, y, z, s; mpfr_exp_t emin, emax; + int i; tests_start_mpfr (); @@ -823,21 +1304,39 @@ test_exact (); - test_overflow1 (); - test_overflow2 (); - test_underflow1 (); - test_underflow2 (); - test_underflow3 (1); + for (i = 0; i <= 2; i++) + { + if (i == 0) + { + /* corresponds to the generic code + mpfr_check_range */ + set_emin (-1024); + set_emax (1024); + } + else if (i == 1) + { + set_emin (MPFR_EMIN_MIN); + set_emax (MPFR_EMAX_MAX); + } + else + { + MPFR_ASSERTN (i == 2); + if (emin == MPFR_EMIN_MIN && emax == MPFR_EMAX_MAX) + break; + set_emin (emin); + set_emax (emax); + } - set_emin (MPFR_EMIN_MIN); - set_emax (MPFR_EMAX_MAX); - test_overflow1 (); - test_overflow2 (); - test_underflow1 (); - test_underflow2 (); - test_underflow3 (2); - set_emin (emin); - set_emax (emax); + test_overflow1 (); + test_overflow2 (); + test_overflow3 (); + test_overflow4 (); + test_overflow5 (); + test_underflow1 (); + test_underflow2 (); + test_underflow3 (i); + test_underflow4 (); + test_underflow5 (); + } tests_end_mpfr (); return 0; diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-04-27 12:47:54.194308761 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-04-27 12:47:54.230308407 +0000 @@ -0,0 +1 @@ +sqr_1n-underflow diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-04-27 12:45:53.171452379 +0000 +++ mpfr-4.0.1-b/VERSION 2018-04-27 12:47:54.226308446 +0000 @@ -1 +1 @@ -4.0.1-p2 +4.0.1-p3 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-04-27 12:45:53.171452379 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-04-27 12:47:54.226308446 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p2" +#define MPFR_VERSION_STRING "4.0.1-p3" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/sqr.c mpfr-4.0.1-b/src/sqr.c --- mpfr-4.0.1-a/src/sqr.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/sqr.c 2018-04-27 12:47:54.218308525 +0000 @@ -161,15 +161,18 @@ if (MPFR_UNLIKELY(ax < __gmpfr_emin)) { /* As seen in mpfr_mul_1, we cannot have a0 = 111...111 here if there - was not exponent decrease (ax--) above. - In the case of an exponent decrease, it is not possible for - GMP_NUMB_BITS=32 since the largest b0 such that b0^2 < 2^(2*32-1) - is b0=3037000499, but its square has only 30 leading ones. - For GMP_NUMB_BITS=64 it is possible: the largest b0 is - 13043817825332782212, and its square has 64 leading ones. */ - if ((ax == __gmpfr_emin - 1) && (ap[0] == ~MPFR_LIMB_HIGHBIT) && - ((rnd_mode == MPFR_RNDN && rb) || - (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb)))) + was not an exponent decrease (ax--) above. + In the case of an exponent decrease: + - For GMP_NUMB_BITS=32, a0 = 111...111 is not possible since the + largest b0 such that b0^2 < 2^(2*32-1) is b0=3037000499, but + its square has only 30 leading ones. + - For GMP_NUMB_BITS=64, a0 = 111...111 is possible: the largest b0 + is 13043817825332782212, and its square has 64 leading ones; but + since the next bit is rb=0, for RNDN, we always have an underflow. + For the test below, note that a is positive. + */ + if (ax == __gmpfr_emin - 1 && ap[0] == MPFR_LIMB_MAX && + MPFR_IS_LIKE_RNDA (rnd_mode, 0)) goto rounding; /* no underflow */ /* For RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2) we have to change to RNDZ. This corresponds to: diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-04-27 12:45:53.171452379 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-04-27 12:47:54.226308446 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p2"; + return "4.0.1-p3"; } diff -Naurd mpfr-4.0.1-a/tests/tsqr.c mpfr-4.0.1-b/tests/tsqr.c --- mpfr-4.0.1-a/tests/tsqr.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tsqr.c 2018-04-27 12:47:54.218308525 +0000 @@ -188,6 +188,156 @@ #endif } +static void +coverage (mpfr_prec_t pmax) +{ + mpfr_prec_t p; + + for (p = MPFR_PREC_MIN; p <= pmax; p++) + { + mpfr_t a, b, c; + int inex; + mpfr_exp_t emin; + + mpfr_init2 (a, p); + mpfr_init2 (b, p); + mpfr_init2 (c, p); + + /* exercise carry in most significant bits of a, with overflow */ + mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ); + mpfr_sqrt (b, b, MPFR_RNDU); + mpfr_div_2exp (c, b, 1, MPFR_RNDN); + mpfr_sqr (c, c, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_sqr (a, b, MPFR_RNDN); + /* if EXP(c) > emax-2, there is overflow */ + if (mpfr_get_exp (c) > mpfr_get_emax () - 2) + { + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); + MPFR_ASSERTN(mpfr_overflow_p ()); + } + else /* no overflow */ + { + /* 2^p-1 is a square only for p=1 */ + MPFR_ASSERTN((p == 1 && inex == 0) || (p > 1 && inex < 0)); + MPFR_ASSERTN(!mpfr_overflow_p ()); + mpfr_set_ui_2exp (c, 1, mpfr_get_emax (), MPFR_RNDZ); + MPFR_ASSERTN(mpfr_equal_p (a, c)); + } + + /* same as above, with RNDU */ + mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ); + mpfr_sqrt (b, b, MPFR_RNDU); + mpfr_div_2exp (c, b, 1, MPFR_RNDN); + mpfr_sqr (c, c, MPFR_RNDU); + mpfr_clear_flags (); + inex = mpfr_sqr (a, b, MPFR_RNDU); + /* if EXP(c) > emax-2, there is overflow */ + if (mpfr_get_exp (c) > mpfr_get_emax () - 2) + { + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); + MPFR_ASSERTN(mpfr_overflow_p ()); + } + else /* no overflow */ + { + /* 2^p-1 is a square only for p=1 */ + MPFR_ASSERTN((p == 1 && inex == 0) || (p > 1 && inex < 0)); + MPFR_ASSERTN(!mpfr_overflow_p ()); + mpfr_set_ui_2exp (c, 1, mpfr_get_emax (), MPFR_RNDZ); + MPFR_ASSERTN(mpfr_equal_p (a, c)); + } + + /* exercise trivial overflow */ + mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ); + mpfr_sqrt (b, b, MPFR_RNDU); + mpfr_mul_2exp (b, b, 1, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_sqr (a, b, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); + MPFR_ASSERTN(mpfr_overflow_p ()); + + /* exercise trivial underflow */ + mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDZ); + mpfr_sqrt (b, b, MPFR_RNDU); + mpfr_div_2exp (b, b, 1, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_sqr (a, b, MPFR_RNDN); + MPFR_ASSERTN(inex < 0); + MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0); + MPFR_ASSERTN(mpfr_underflow_p ()); + + /* exercise square between 0.5*2^emin and its predecessor (emin even) */ + emin = mpfr_get_emin (); + mpfr_set_emin (emin + (emin & 1)); /* now emin is even */ + mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDN); + inex = mpfr_sqrt (b, b, MPFR_RNDZ); + MPFR_ASSERTN(inex != 0); /* sqrt(2) is not exact */ + mpfr_mul_2exp (c, b, 1, MPFR_RNDN); + mpfr_sqr (c, c, MPFR_RNDN); + mpfr_clear_flags (); + inex = mpfr_sqr (a, b, MPFR_RNDN); + if (mpfr_get_exp (c) < mpfr_get_emin () + 2) /* underflow */ + { + /* if c > 0.5*2^(emin+1), we should round to 0.5*2^emin */ + if (mpfr_cmp_ui_2exp (c, 1, mpfr_get_emin ()) > 0) + { + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0); + MPFR_ASSERTN(mpfr_underflow_p ()); + } + else /* we should round to 0 */ + { + MPFR_ASSERTN(inex < 0); + MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0); + MPFR_ASSERTN(mpfr_underflow_p ()); + } + } + else + { + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0); + MPFR_ASSERTN(!mpfr_underflow_p ()); + } + mpfr_set_emin (emin); + + /* exercise exact square root 2^(emin-2) for emin even */ + emin = mpfr_get_emin (); + mpfr_set_emin (emin + (emin & 1)); /* now emin is even */ + mpfr_set_ui_2exp (b, 1, (mpfr_get_emin () - 2) / 2, MPFR_RNDN); + inex = mpfr_sqr (a, b, MPFR_RNDN); + MPFR_ASSERTN(inex < 0); + MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0); + MPFR_ASSERTN(mpfr_underflow_p ()); + mpfr_set_emin (emin); + + /* same as above, for RNDU */ + emin = mpfr_get_emin (); + mpfr_set_emin (emin + (emin & 1)); /* now emin is even */ + mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDN); + inex = mpfr_sqrt (b, b, MPFR_RNDZ); + MPFR_ASSERTN(inex != 0); /* sqrt(2) is not exact */ + mpfr_mul_2exp (c, b, 1, MPFR_RNDN); + mpfr_sqr (c, c, MPFR_RNDU); + mpfr_clear_flags (); + inex = mpfr_sqr (a, b, MPFR_RNDU); + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0); + /* we have underflow if c < 2^(emin+1) */ + if (mpfr_cmp_ui_2exp (c, 1, mpfr_get_emin () + 1) < 0) + MPFR_ASSERTN(mpfr_underflow_p ()); + else + MPFR_ASSERTN(!mpfr_underflow_p ()); + mpfr_set_emin (emin); + + mpfr_clear (a); + mpfr_clear (b); + mpfr_clear (c); + } +} + int main (void) { @@ -195,6 +345,7 @@ tests_start_mpfr (); + coverage (1024); check_mpn_sqr (); check_special (); test_underflow (); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-04-27 12:50:10.592974822 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-04-27 12:50:10.624974512 +0000 @@ -0,0 +1 @@ +get_str diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-04-27 12:47:54.226308446 +0000 +++ mpfr-4.0.1-b/VERSION 2018-04-27 12:50:10.624974512 +0000 @@ -1 +1 @@ -4.0.1-p3 +4.0.1-p4 diff -Naurd mpfr-4.0.1-a/doc/mpfr.info mpfr-4.0.1-b/doc/mpfr.info --- mpfr-4.0.1-a/doc/mpfr.info 2018-02-07 12:58:03.000000000 +0000 +++ mpfr-4.0.1-b/doc/mpfr.info 2018-04-27 12:50:17.128911341 +0000 @@ -1321,10 +1321,9 @@ size_t N, mpfr_t OP, mpfr_rnd_t RND) Convert OP to a string of digits in base B, with rounding in the direction RND, where N is either zero (see below) or the number of - significant digits output in the string; in the latter case, N must - be greater or equal to 2. The base may vary from 2 to 62; - otherwise the function does nothing and immediately returns a null - pointer. + significant digits output in the string. The base may vary from 2 + to 62; otherwise the function does nothing and immediately returns + a null pointer. If the input is NaN, then the returned string is ‘@NaN@’ and the NaN flag is set. If the input is +Inf (resp. −Inf), then the @@ -4471,21 +4470,21 @@ * mpfr_expm1: Special Functions. (line 40) * mpfr_fac_ui: Special Functions. (line 136) * mpfr_fits_intmax_p: Conversion Functions. - (line 168) + (line 167) * mpfr_fits_sint_p: Conversion Functions. - (line 164) + (line 163) * mpfr_fits_slong_p: Conversion Functions. - (line 162) + (line 161) * mpfr_fits_sshort_p: Conversion Functions. - (line 166) + (line 165) * mpfr_fits_uintmax_p: Conversion Functions. - (line 167) + (line 166) * mpfr_fits_uint_p: Conversion Functions. - (line 163) + (line 162) * mpfr_fits_ulong_p: Conversion Functions. - (line 161) + (line 160) * mpfr_fits_ushort_p: Conversion Functions. - (line 165) + (line 164) * mpfr_flags_clear: Exception Related Functions. (line 190) * mpfr_flags_restore: Exception Related Functions. @@ -4520,7 +4519,7 @@ * mpfr_free_cache2: Special Functions. (line 295) * mpfr_free_pool: Special Functions. (line 309) * mpfr_free_str: Conversion Functions. - (line 156) + (line 155) * mpfr_frexp: Conversion Functions. (line 49) * mpfr_gamma: Special Functions. (line 155) @@ -4928,30 +4927,30 @@ Node: Assignment Functions47553 Node: Combined Initialization and Assignment Functions57499 Node: Conversion Functions58800 -Node: Basic Arithmetic Functions69381 -Node: Comparison Functions80277 -Node: Special Functions83765 -Node: Input and Output Functions101974 -Node: Formatted Output Functions106751 -Node: Integer and Remainder Related Functions116956 -Node: Rounding-Related Functions124484 -Node: Miscellaneous Functions131001 -Node: Exception Related Functions141493 -Node: Compatibility with MPF151733 -Node: Custom Interface154679 -Node: Internals159310 -Node: API Compatibility160854 -Node: Type and Macro Changes162802 -Node: Added Functions165985 -Node: Changed Functions170499 -Node: Removed Functions177095 -Node: Other Changes177825 -Node: MPFR and the IEEE 754 Standard179526 -Node: Contributors182143 -Node: References185200 -Node: GNU Free Documentation License187084 -Node: Concept Index209677 -Node: Function and Type Index216049 +Node: Basic Arithmetic Functions69323 +Node: Comparison Functions80219 +Node: Special Functions83707 +Node: Input and Output Functions101916 +Node: Formatted Output Functions106693 +Node: Integer and Remainder Related Functions116898 +Node: Rounding-Related Functions124426 +Node: Miscellaneous Functions130943 +Node: Exception Related Functions141435 +Node: Compatibility with MPF151675 +Node: Custom Interface154621 +Node: Internals159252 +Node: API Compatibility160796 +Node: Type and Macro Changes162744 +Node: Added Functions165927 +Node: Changed Functions170441 +Node: Removed Functions177037 +Node: Other Changes177767 +Node: MPFR and the IEEE 754 Standard179468 +Node: Contributors182085 +Node: References185142 +Node: GNU Free Documentation License187026 +Node: Concept Index209619 +Node: Function and Type Index215991  End Tag Table diff -Naurd mpfr-4.0.1-a/doc/mpfr.texi mpfr-4.0.1-b/doc/mpfr.texi --- mpfr-4.0.1-a/doc/mpfr.texi 2018-02-07 12:50:31.000000000 +0000 +++ mpfr-4.0.1-b/doc/mpfr.texi 2018-04-27 12:50:10.612974628 +0000 @@ -1655,8 +1655,8 @@ @deftypefun {char *} mpfr_get_str (char *@var{str}, mpfr_exp_t *@var{expptr}, int @var{b}, size_t @var{n}, mpfr_t @var{op}, mpfr_rnd_t @var{rnd}) Convert @var{op} to a string of digits in base @var{b}, with rounding in the direction @var{rnd}, where @var{n} is either zero (see below) or the -number of significant digits output in the string; in the latter case, -@var{n} must be greater or equal to 2. The base may vary from 2 to 62; +number of significant digits output in the string. +The base may vary from 2 to 62; otherwise the function does nothing and immediately returns a null pointer. If the input is NaN, then the returned string is @samp{@@NaN@@} and the @@ -1699,8 +1699,7 @@ but in some very rare cases, it might be @math{m+1} (the smallest case for bases up to 62 is when @var{p} equals 186564318007 for bases 7 and 49). -@c In the source src/get_str.c, this is due to the approximate mpfr_ceil_mul, -@c but also m = 1 is changed to 2. +@c In the source src/get_str.c, this is due to the approximate mpfr_ceil_mul. If @var{str} is a null pointer, space for the significand is allocated using the allocation function (@pxref{Memory Handling}) and a pointer to the string diff -Naurd mpfr-4.0.1-a/src/get_str.c mpfr-4.0.1-b/src/get_str.c --- mpfr-4.0.1-a/src/get_str.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/get_str.c 2018-04-27 12:50:10.612974628 +0000 @@ -2325,15 +2325,12 @@ */ m = 1 + mpfr_ceil_mul (IS_POW2(b) ? MPFR_PREC(x) - 1 : MPFR_PREC(x), b, 1); - if (m < 2) - m = 2; } MPFR_LOG_MSG (("m=%zu\n", m)); - /* The code below for non-power-of-two bases works for m=1; - this is important for the internal use of mpfr_get_str. */ - MPFR_ASSERTN (m >= 2 || (!IS_POW2(b) && m >= 1)); + /* The code below works for m=1, both for power-of-two and non-power-of-two + bases; this is important for the internal use of mpfr_get_str. */ /* x is a floating-point number */ @@ -2376,6 +2373,8 @@ /* the first digit will contain only r bits */ prec = (m - 1) * pow2 + r; /* total number of bits */ + /* if m=1 then 1 <= prec <= pow2, and since prec=1 is now valid in MPFR, + the power-of-two code also works for m=1 */ n = MPFR_PREC2LIMBS (prec); MPFR_TMP_MARK (marker); diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-04-27 12:47:54.226308446 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-04-27 12:50:10.620974551 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p3" +#define MPFR_VERSION_STRING "4.0.1-p4" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-04-27 12:47:54.226308446 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-04-27 12:50:10.624974512 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p3"; + return "4.0.1-p4"; } diff -Naurd mpfr-4.0.1-a/tests/tget_str.c mpfr-4.0.1-b/tests/tget_str.c --- mpfr-4.0.1-a/tests/tget_str.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tget_str.c 2018-04-27 12:50:10.612974628 +0000 @@ -1267,6 +1267,41 @@ #define ITER 1000 +static void +coverage (void) +{ + mpfr_t x; + char s[42]; + mpfr_exp_t e; + int b = 3; + size_t m = 40; + + mpfr_init2 (x, 128); + + /* exercise corner case in mpfr_get_str_aux: exact case (e < 0), where r + rounds to a power of 2, and f is a multiple of GMP_NUMB_BITS */ + mpfr_set_ui_2exp (x, 1, 64, MPFR_RNDU); + mpfr_nextbelow (x); + /* x = 2^64 - 2^(-64) */ + mpfr_get_str (s, &e, b, m, x, MPFR_RNDU); + /* s is the base-3 string for 6148914691236517206 (in base 10) */ + MPFR_ASSERTN(strcmp (s, "1111222002212212010121102012021021021200") == 0); + MPFR_ASSERTN(e == 41); + + /* exercise corner case in mpfr_get_str: input is m=0, then it is changed + to m=1 */ + mpfr_set_prec (x, 1); + mpfr_set_ui (x, 1, MPFR_RNDN); + mpfr_get_str (s, &e, 2, 0, x, MPFR_RNDN); + MPFR_ASSERTN(strcmp (s, "1") == 0); + MPFR_ASSERTN(e == 1); + mpfr_get_str (s, &e, 2, 1, x, MPFR_RNDN); + MPFR_ASSERTN(strcmp (s, "1") == 0); + MPFR_ASSERTN(e == 1); + + mpfr_clear (x); +} + int main (int argc, char *argv[]) { @@ -1281,6 +1316,7 @@ tests_start_mpfr (); + coverage (); check_small (); check_special (2, 2); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-04-27 12:52:13.875783093 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-04-27 12:52:13.911782747 +0000 @@ -0,0 +1 @@ +cmp_q-special diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-04-27 12:50:10.624974512 +0000 +++ mpfr-4.0.1-b/VERSION 2018-04-27 12:52:13.911782747 +0000 @@ -1 +1 @@ -4.0.1-p4 +4.0.1-p5 diff -Naurd mpfr-4.0.1-a/src/gmp_op.c mpfr-4.0.1-b/src/gmp_op.c --- mpfr-4.0.1-a/src/gmp_op.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/gmp_op.c 2018-04-27 12:52:13.899782862 +0000 @@ -452,11 +452,15 @@ mpfr_prec_t p; MPFR_SAVE_EXPO_DECL (expo); - if (MPFR_UNLIKELY (mpq_denref (q) == 0)) + if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (q)) == 0)) { /* q is an infinity or NaN */ - mpfr_init2 (t, 2); + mpfr_flags_t old_flags; + + mpfr_init2 (t, MPFR_PREC_MIN); + old_flags = __gmpfr_flags; mpfr_set_q (t, q, MPFR_RNDN); + __gmpfr_flags = old_flags; res = mpfr_cmp (x, t); mpfr_clear (t); return res; diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-04-27 12:50:10.620974551 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-04-27 12:52:13.907782785 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p4" +#define MPFR_VERSION_STRING "4.0.1-p5" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-04-27 12:50:10.624974512 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-04-27 12:52:13.911782747 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p4"; + return "4.0.1-p5"; } diff -Naurd mpfr-4.0.1-a/tests/tgmpop.c mpfr-4.0.1-b/tests/tgmpop.c --- mpfr-4.0.1-a/tests/tgmpop.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tgmpop.c 2018-04-27 12:52:13.899782862 +0000 @@ -307,16 +307,39 @@ mpfr_init2 (z, MPFR_PREC_MIN); mpq_init (y); - /* check the erange flag when x is NaN */ + /* Check the flags when x is NaN: the erange flags must be set, and + only this one. */ mpfr_set_nan (x); mpq_set_ui (y, 17, 1); - mpfr_clear_erangeflag (); + mpfr_clear_flags (); res1 = mpfr_cmp_q (x, y); - if (res1 != 0 || mpfr_erangeflag_p () == 0) + if (res1 != 0 || __gmpfr_flags != MPFR_FLAGS_ERANGE) { printf ("Error for mpfr_cmp_q (NaN, 17)\n"); printf ("Return value: expected 0, got %d\n", res1); - printf ("Erange flag: expected set, got %d\n", mpfr_erangeflag_p ()); + printf ("Expected flags:"); + flags_out (MPFR_FLAGS_ERANGE); + printf ("Got flags: "); + flags_out (__gmpfr_flags); + exit (1); + } + + /* Check the flags when y is NaN: the erange flags must be set, and + only this one. */ + mpfr_set_ui (x, 42, MPFR_RNDN); + /* A NaN rational is represented by 0/0 (MPFR extension). */ + mpz_set_ui (mpq_numref (y), 0); + mpz_set_ui (mpq_denref (y), 0); + mpfr_clear_flags (); + res1 = mpfr_cmp_q (x, y); + if (res1 != 0 || __gmpfr_flags != MPFR_FLAGS_ERANGE) + { + printf ("Error for mpfr_cmp_q (42, NaN)\n"); + printf ("Return value: expected 0, got %d\n", res1); + printf ("Expected flags:"); + flags_out (MPFR_FLAGS_ERANGE); + printf ("Got flags: "); + flags_out (__gmpfr_flags); exit (1); } @@ -341,6 +364,55 @@ } } } + + /* check for y = 1/0 */ + mpz_set_ui (mpq_numref (y), 1); + mpz_set_ui (mpq_denref (y), 0); + mpfr_set_ui (x, 1, MPFR_RNDN); + MPFR_ASSERTN(mpfr_cmp_q (x, y) < 0); + mpfr_set_inf (x, -1); + MPFR_ASSERTN(mpfr_cmp_q (x, y) < 0); + mpfr_set_inf (x, +1); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + mpfr_set_nan (x); + mpfr_clear_erangeflag (); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + MPFR_ASSERTN(mpfr_erangeflag_p ()); + + /* check for y = -1/0 */ + mpz_set_si (mpq_numref (y), -1); + mpz_set_ui (mpq_denref (y), 0); + mpfr_set_ui (x, 1, MPFR_RNDN); + MPFR_ASSERTN(mpfr_cmp_q (x, y) > 0); + mpfr_set_inf (x, -1); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + mpfr_set_inf (x, +1); + MPFR_ASSERTN(mpfr_cmp_q (x, y) > 0); + mpfr_set_nan (x); + mpfr_clear_erangeflag (); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + MPFR_ASSERTN(mpfr_erangeflag_p ()); + + /* check for y = 0/0 */ + mpz_set_ui (mpq_numref (y), 0); + mpz_set_ui (mpq_denref (y), 0); + mpfr_set_ui (x, 1, MPFR_RNDN); + mpfr_clear_erangeflag (); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + MPFR_ASSERTN(mpfr_erangeflag_p ()); + mpfr_set_inf (x, -1); + mpfr_clear_erangeflag (); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + MPFR_ASSERTN(mpfr_erangeflag_p ()); + mpfr_set_inf (x, +1); + mpfr_clear_erangeflag (); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + MPFR_ASSERTN(mpfr_erangeflag_p ()); + mpfr_set_nan (x); + mpfr_clear_erangeflag (); + MPFR_ASSERTN(mpfr_cmp_q (x, y) == 0); + MPFR_ASSERTN(mpfr_erangeflag_p ()); + mpq_clear (y); mpfr_clear (x); mpfr_clear (z); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-04-27 12:54:17.862594948 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-04-27 12:54:17.894594642 +0000 @@ -0,0 +1 @@ +io-null-stream diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-04-27 12:52:13.911782747 +0000 +++ mpfr-4.0.1-b/VERSION 2018-04-27 12:54:17.894594642 +0000 @@ -1 +1 @@ -4.0.1-p5 +4.0.1-p6 diff -Naurd mpfr-4.0.1-a/src/inp_str.c mpfr-4.0.1-b/src/inp_str.c --- mpfr-4.0.1-a/src/inp_str.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/inp_str.c 2018-04-27 12:54:17.882594757 +0000 @@ -37,9 +37,6 @@ int retval; size_t nread; - if (stream == NULL) - stream = stdin; - alloc_size = 100; str = (unsigned char *) mpfr_allocate_func (alloc_size); str_size = 0; diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-04-27 12:52:13.907782785 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-04-27 12:54:17.894594642 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p5" +#define MPFR_VERSION_STRING "4.0.1-p6" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/out_str.c mpfr-4.0.1-b/src/out_str.c --- mpfr-4.0.1-a/src/out_str.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/out_str.c 2018-04-27 12:54:17.882594757 +0000 @@ -43,10 +43,6 @@ MPFR_ASSERTN (base >= 2 && base <= 62); - /* when stream=NULL, output to stdout */ - if (stream == NULL) - stream = stdout; - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (op))) { if (MPFR_IS_NAN (op)) diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-04-27 12:52:13.911782747 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-04-27 12:54:17.894594642 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p5"; + return "4.0.1-p6"; } diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-10 15:29:07.937776370 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-10 15:29:08.017776303 +0000 @@ -0,0 +1 @@ +tstckintc-casts diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-04-27 12:54:17.894594642 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-10 15:29:08.017776303 +0000 @@ -1 +1 @@ -4.0.1-p6 +4.0.1-p7 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-04-27 12:54:17.894594642 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-10 15:29:08.013776307 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p6" +#define MPFR_VERSION_STRING "4.0.1-p7" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-04-27 12:54:17.894594642 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-10 15:29:08.017776303 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p6"; + return "4.0.1-p7"; } diff -Naurd mpfr-4.0.1-a/tests/tstckintc.c mpfr-4.0.1-b/tests/tstckintc.c --- mpfr-4.0.1-a/tests/tstckintc.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tstckintc.c 2018-07-10 15:29:07.989776327 +0000 @@ -32,6 +32,22 @@ #define ALIGNED(s) (((s) + sizeof (long) - 1) / sizeof (long) * sizeof (long)) +/* This code ensures alignment to "long". However, this might not be + sufficient on some platforms. GCC's -Wcast-align=strict option can + be useful, but this needs successive casts to help GCC, e.g. + + newx = (mpfr_ptr) (long *) (void *) old_stack; + + successively casts old_stack (of type char *) to + - void *: avoid a false positive for the following cast to long * + (as the code takes care of alignment to "long"); + - long *: this corresponds to the alignment checked by MPFR; coming + from void *, it will not trigger a warning (even if incorrect); + - mpfr_ptr: -Wcast-align=strict will emit a warning if mpfr_ptr has + an alignment requirement stronger than long *. In such a case, + the code will have to be fixed. +*/ + static void * new_st (size_t s) { @@ -94,12 +110,14 @@ void *mantissa = mpfr_custom_get_significand (x); size_t size_mantissa = mpfr_custom_get_size (mpfr_get_prec (x)); mpfr_ptr newx; + long *newx2; memmove (old_stack, x, sizeof (mpfr_t)); memmove (old_stack + ALIGNED (sizeof (mpfr_t)), mantissa, size_mantissa); - newx = (mpfr_ptr) old_stack; - mpfr_custom_move (newx, old_stack + ALIGNED (sizeof (mpfr_t))); - stack = old_stack + ALIGNED (sizeof (mpfr_t)) + ALIGNED (size_mantissa); + newx = (mpfr_ptr) (long *) (void *) old_stack; + newx2 = (long *) (void *) (old_stack + ALIGNED (sizeof (mpfr_t))); + mpfr_custom_move (newx, newx2); + stack = (char *) newx2 + ALIGNED (size_mantissa); return newx; } @@ -113,7 +131,7 @@ memmove (old_stack, x, sizeof (mpfr_t)); memmove (old_stack + ALIGNED (sizeof (mpfr_t)), mantissa, size_mantissa); - newx = (mpfr_ptr) old_stack; + newx = (mpfr_ptr) (long *) (void *) old_stack; (mpfr_custom_move) (newx, old_stack + ALIGNED (sizeof (mpfr_t))); stack = old_stack + ALIGNED (sizeof (mpfr_t)) + ALIGNED (size_mantissa); return newx; @@ -127,7 +145,7 @@ mpfr_ptr x, y; reset_stack (); - org = (long *) stack; + org = (long *) (void *) stack; x = new_mpfr (p); y = new_mpfr (p); @@ -277,7 +295,7 @@ long *a, *b, *c; reset_stack (); - org = (long *) stack; + org = (long *) (void *) stack; a = dummy_set_si (42); b = dummy_set_si (17); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-10 15:33:45.189532503 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-10 15:33:45.265532432 +0000 @@ -0,0 +1 @@ +set_d64-ternary diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-10 15:29:08.017776303 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-10 15:33:45.261532435 +0000 @@ -1 +1 @@ -4.0.1-p7 +4.0.1-p8 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-10 15:29:08.013776307 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-10 15:33:45.261532435 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p7" +#define MPFR_VERSION_STRING "4.0.1-p8" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/set_d64.c mpfr-4.0.1-b/src/set_d64.c --- mpfr-4.0.1-a/src/set_d64.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/set_d64.c 2018-07-10 15:33:45.237532458 +0000 @@ -425,7 +425,7 @@ 1 character for terminating \0. */ decimal64_to_string (s, d); - return mpfr_set_str (r, s, 10, rnd_mode); + return mpfr_strtofr (r, s, NULL, 10, rnd_mode); } #endif /* MPFR_WANT_DECIMAL_FLOATS */ diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-10 15:29:08.017776303 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-10 15:33:45.261532435 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p7"; + return "4.0.1-p8"; } diff -Naurd mpfr-4.0.1-a/tests/tget_set_d64.c mpfr-4.0.1-b/tests/tget_set_d64.c --- mpfr-4.0.1-a/tests/tget_set_d64.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tget_set_d64.c 2018-07-10 15:33:45.237532458 +0000 @@ -381,6 +381,66 @@ mpfr_clear (x); } +static void +powers_of_10 (void) +{ + mpfr_t x1, x2; + _Decimal64 d[2]; + int i, rnd; + unsigned int neg; + + mpfr_inits2 (200, x1, x2, (mpfr_ptr) 0); + for (i = 0, d[0] = 1, d[1] = 1; i < 150; i++, d[0] *= 10, d[1] /= 10) + for (neg = 0; neg <= 3; neg++) + RND_LOOP_NO_RNDF (rnd) + { + int inex1, inex2; + mpfr_flags_t flags1, flags2; + mpfr_rnd_t rx1; + _Decimal64 dd; + + inex1 = mpfr_set_si (x1, (neg >> 1) ? -i : i, MPFR_RNDN); + MPFR_ASSERTN (inex1 == 0); + + rx1 = (neg & 1) ? + MPFR_INVERT_RND ((mpfr_rnd_t) rnd) : (mpfr_rnd_t) rnd; + mpfr_clear_flags (); + inex1 = mpfr_exp10 (x1, x1, rx1); + flags1 = __gmpfr_flags; + + dd = d[neg >> 1]; + + if (neg & 1) + { + MPFR_SET_NEG (x1); + inex1 = -inex1; + dd = -dd; + } + + mpfr_clear_flags (); + inex2 = mpfr_set_decimal64 (x2, dd, (mpfr_rnd_t) rnd); + flags2 = __gmpfr_flags; + + if (!(mpfr_equal_p (x1, x2) && + SAME_SIGN (inex1, inex2) && + flags1 == flags2)) + { + printf ("Error in powers_of_10 for i=%d, neg=%d, %s\n", + i, neg, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd)); + 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, (mpfr_ptr) 0); +} + int main (void) { @@ -401,6 +461,7 @@ check_overflow (); #endif check_tiny (); + powers_of_10 (); tests_end_mpfr (); return 0; diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-10 15:46:13.596797032 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-10 15:46:13.676796949 +0000 @@ -0,0 +1 @@ +buffer_sandwich diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-10 15:33:45.261532435 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-10 15:46:13.676796949 +0000 @@ -1 +1 @@ -4.0.1-p8 +4.0.1-p9 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-10 15:33:45.261532435 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-10 15:46:13.672796954 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p8" +#define MPFR_VERSION_STRING "4.0.1-p9" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/vasprintf.c mpfr-4.0.1-b/src/vasprintf.c --- mpfr-4.0.1-a/src/vasprintf.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/vasprintf.c 2018-07-10 15:46:13.648796978 +0000 @@ -683,20 +683,31 @@ else { const size_t step = 3; - const size_t size = len + tz; - const size_t r = size % step == 0 ? step : size % step; - const size_t q = size % step == 0 ? size / step - 1 : size / step; - const size_t fullsize = size + q; - size_t i; + size_t size, q, r, fullsize; + /* check that len + tz does not overflow */ + if (len > (size_t) -1 - tz) + return 1; + + size = len + tz; /* number of digits */ MPFR_ASSERTD (size > 0); + q = (size - 1) / step; /* number of separators C */ + r = ((size - 1) % step) + 1; /* number of digits in the leftmost block */ + + /* check that size + q does not overflow */ + if (size > (size_t) -1 - q) + return 1; + + fullsize = size + q; /* number of digits and separators */ + if (buffer_incr_len (b, fullsize)) return 1; if (b->size != 0) { char *oldcurr; + size_t i; MPFR_ASSERTD (*b->curr == '\0'); MPFR_ASSERTN (b->size < ((size_t) -1) - fullsize); @@ -705,11 +716,21 @@ MPFR_DBGRES (oldcurr = b->curr); - /* first R significant digits */ - memcpy (b->curr, str, r); + /* first r significant digits (leftmost block) */ + if (r <= len) + { + memcpy (b->curr, str, r); + str += r; + len -= r; + } + else + { + MPFR_ASSERTD (r > len); + memcpy (b->curr, str, len); + memset (b->curr + len, '0', r - len); + len = 0; + } b->curr += r; - str += r; - len -= r; /* blocks of thousands. Warning: STR might end in the middle of a block */ for (i = 0; i < q; ++i) @@ -722,6 +743,7 @@ { memcpy (b->curr, str, step); len -= step; + str += step; } else /* last digits in STR, fill up thousand block with zeros */ @@ -736,7 +758,6 @@ memset (b->curr, '0', step); b->curr += step; - str += step; } MPFR_ASSERTD (b->curr - oldcurr == fullsize); @@ -1920,8 +1941,14 @@ /* integral part (may also be "nan" or "inf") */ MPFR_ASSERTN (np.ip_ptr != NULL); /* never empty */ if (MPFR_UNLIKELY (np.thousands_sep)) - buffer_sandwich (buf, np.ip_ptr, np.ip_size, np.ip_trailing_zeros, - np.thousands_sep); + { + if (buffer_sandwich (buf, np.ip_ptr, np.ip_size, np.ip_trailing_zeros, + np.thousands_sep)) + { + buf->len = -1; + goto clear_and_exit; + } + } else { buffer_cat (buf, np.ip_ptr, np.ip_size); diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-10 15:33:45.261532435 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-10 15:46:13.672796954 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p8"; + return "4.0.1-p9"; } diff -Naurd mpfr-4.0.1-a/tests/tprintf.c mpfr-4.0.1-b/tests/tprintf.c --- mpfr-4.0.1-a/tests/tprintf.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tprintf.c 2018-07-10 15:46:13.648796978 +0000 @@ -31,6 +31,10 @@ #include #include +#ifdef HAVE_LOCALE_H +#include +#endif + #include "mpfr-intmax.h" #include "mpfr-test.h" #define STDOUT_FILENO 1 @@ -474,30 +478,41 @@ mpfr_clear (x); } -#ifdef HAVE_LOCALE_H - -#include - -const char * const tab_locale[] = { - "en_US", - "en_US.iso88591", - "en_US.iso885915", - "en_US.utf8" -}; +#if defined(HAVE_LOCALE_H) && defined(HAVE_SETLOCALE) static void test_locale (void) { + const char * const tab_locale[] = { + "en_US", + "en_US.iso88591", + "en_US.iso885915", + "en_US.utf8" + }; int i; - char *s = NULL; mpfr_t x; int count; + char v[] = "99999999999999999999999.5"; - for(i = 0; i < numberof(tab_locale) && s == NULL; i++) - s = setlocale (LC_ALL, tab_locale[i]); + for (i = 0; i < numberof(tab_locale); i++) + { + char *s; - if (s == NULL || MPFR_THOUSANDS_SEPARATOR != ',') - return; + s = setlocale (LC_ALL, tab_locale[i]); + + if (s != NULL && MPFR_THOUSANDS_SEPARATOR == ',') + break; + } + + if (i == numberof(tab_locale)) + { + if (getenv ("MPFR_CHECK_LOCALES") == NULL) + return; + + fprintf (stderr, "Cannot find a locale with ',' thousands separator.\n" + "Please install one of the en_US based locales.\n"); + exit (1); + } mpfr_init2 (x, 113); mpfr_set_ui (x, 10000, MPFR_RNDN); @@ -507,6 +522,50 @@ count = mpfr_printf ("(2) 10000=%'Rf \n", x); check_length (10001, count, 25, d); + mpfr_set_ui (x, 1000, MPFR_RNDN); + count = mpfr_printf ("(3) 1000=%'Rf \n", x); + check_length (10002, count, 23, d); + + for (i = 1; i <= sizeof (v) - 3; i++) + { + mpfr_set_str (x, v + sizeof (v) - 3 - i, 10, MPFR_RNDN); + count = mpfr_printf ("(4) 10^i=%'.0Rf \n", x); + check_length (10002 + i, count, 12 + i + i/3, d); + } + +#define N0 20 + + for (i = 1; i <= N0; i++) + { + char s[N0+4]; + int j, rnd; + + s[0] = '1'; + for (j = 1; j <= i; j++) + s[j] = '0'; + s[i+1] = '\0'; + + mpfr_set_str (x, s, 10, MPFR_RNDN); + + RND_LOOP (rnd) + { + count = mpfr_printf ("(5) 10^i=%'.0R*f \n", (mpfr_rnd_t) rnd, x); + check_length (11000 + 10 * i + rnd, count, 12 + i + i/3, d); + } + + strcat (s + (i + 1), ".5"); + count = mpfr_printf ("(5) 10^i=%'.0Rf \n", x); + check_length (11000 + 10 * i + 9, count, 12 + i + i/3, d); + } + + mpfr_set_str (x, "1000", 10, MPFR_RNDN); + count = mpfr_printf ("%'012.3Rg\n", x); + check_length (12000, count, 13, d); + count = mpfr_printf ("%'012.4Rg\n", x); + check_length (12001, count, 13, d); + count = mpfr_printf ("%'013.4Rg\n", x); + check_length (12002, count, 14, d); + mpfr_clear (x); } @@ -515,7 +574,11 @@ static void test_locale (void) { - /* Nothing */ + if (getenv ("MPFR_CHECK_LOCALES") != NULL) + { + fprintf (stderr, "Cannot test locales.\n"); + exit (1); + } } #endif diff -Naurd mpfr-4.0.1-a/tests/tsprintf.c mpfr-4.0.1-b/tests/tsprintf.c --- mpfr-4.0.1-a/tests/tsprintf.c 2018-01-10 10:15:30.000000000 +0000 +++ mpfr-4.0.1-b/tests/tsprintf.c 2018-07-10 15:46:13.648796978 +0000 @@ -380,7 +380,7 @@ check_sprintf ("1.00 ", "%-#20.3RG", x); check_sprintf ("0.9999 ", "%-#20.4RG", x); - /* multiple of 10 */ + /* powers of 10 */ mpfr_set_str (x, "1e17", 10, MPFR_RNDN); check_sprintf ("1e+17", "%Re", x); check_sprintf ("1.000e+17", "%.3Re", x); @@ -402,7 +402,7 @@ check_sprintf ("1", "%.0RUf", x); check_sprintf ("1", "%.0RYf", x); - /* multiple of 10 with 'g' style */ + /* powers of 10 with 'g' style */ mpfr_set_str (x, "10", 10, MPFR_RNDN); check_sprintf ("10", "%Rg", x); check_sprintf ("1e+01", "%.0Rg", x); @@ -419,6 +419,12 @@ check_sprintf ("1e+03", "%.0Rg", x); check_sprintf ("1e+03", "%.3Rg", x); check_sprintf ("1000", "%.4Rg", x); + check_sprintf ("1e+03", "%.3Rg", x); + check_sprintf ("1000", "%.4Rg", x); + check_sprintf (" 1e+03", "%9.3Rg", x); + check_sprintf (" 1000", "%9.4Rg", x); + check_sprintf ("00001e+03", "%09.3Rg", x); + check_sprintf ("000001000", "%09.4Rg", x); mpfr_ui_div (x, 1, x, MPFR_RNDN); check_sprintf ("0.001", "%Rg", x); @@ -430,6 +436,10 @@ check_sprintf ("1e+05", "%.0Rg", x); check_sprintf ("1e+05", "%.5Rg", x); check_sprintf ("100000", "%.6Rg", x); + check_sprintf (" 1e+05", "%17.5Rg", x); + check_sprintf (" 100000", "%17.6Rg", x); + check_sprintf ("0000000000001e+05", "%017.5Rg", x); + check_sprintf ("00000000000100000", "%017.6Rg", x); mpfr_ui_div (x, 1, x, MPFR_RNDN); check_sprintf ("1e-05", "%Rg", x); @@ -857,6 +867,12 @@ "%.*Zi, %R*e, %Lf", 20, mpz, rnd, x, d); #endif + /* check invalid spec.spec */ + check_vsprintf ("%,", "%,"); + + /* check empty format */ + check_vsprintf ("%", "%"); + mpf_clear (mpf); mpq_clear (mpq); mpz_clear (mpz); @@ -864,12 +880,12 @@ return 0; } -#if MPFR_LCONV_DPTS +#if defined(HAVE_LOCALE_H) && defined(HAVE_SETLOCALE) && MPFR_LCONV_DPTS /* Check with locale "da_DK". On most platforms, decimal point is ',' and thousands separator is '.'; the test is not performed if this is not the case or if the locale doesn't exist. */ -static int +static void locale_da_DK (void) { mpfr_prec_t p = 128; @@ -878,7 +894,16 @@ if (setlocale (LC_ALL, "da_DK") == 0 || localeconv()->decimal_point[0] != ',' || localeconv()->thousands_sep[0] != '.') - return 0; + { + setlocale (LC_ALL, "C"); + + if (getenv ("MPFR_CHECK_LOCALES") == NULL) + return; + + fprintf (stderr, + "Cannot test the da_DK locale (not found or inconsistent).\n"); + exit (1); + } mpfr_init2 (x, p); @@ -917,10 +942,11 @@ check_sprintf ("100" S2 "0000", "%'.4Rf", x); mpfr_clear (x); - return 0; + + setlocale (LC_ALL, "C"); } -#endif /* MPFR_LCONV_DPTS */ +#endif /* ... && MPFR_LCONV_DPTS */ /* check concordance between mpfr_asprintf result with a regular mpfr float and with a regular double float */ @@ -1425,6 +1451,117 @@ exit (1); } +#if defined(HAVE_LOCALE_H) && defined(HAVE_SETLOCALE) + +/* The following tests should be equivalent to those from test_locale() + in tprintf.c (remove the \n at the end of the test strings). */ + +static void +test_locale (void) +{ + const char * const tab_locale[] = { + "en_US", + "en_US.iso88591", + "en_US.iso885915", + "en_US.utf8" + }; + int i; + mpfr_t x; + char v[] = "99999999999999999999999.5"; + + for (i = 0; i < numberof(tab_locale); i++) + { + char *s; + + s = setlocale (LC_ALL, tab_locale[i]); + + if (s != NULL && MPFR_THOUSANDS_SEPARATOR == ',') + break; + } + + if (i == numberof(tab_locale)) + { + if (getenv ("MPFR_CHECK_LOCALES") == NULL) + return; + + fprintf (stderr, "Cannot find a locale with ',' thousands separator.\n" + "Please install one of the en_US based locales.\n"); + exit (1); + } + + mpfr_init2 (x, 113); + mpfr_set_ui (x, 10000, MPFR_RNDN); + + check_sprintf ("(1) 10000=10,000 ", "(1) 10000=%'Rg ", x); + check_sprintf ("(2) 10000=10,000.000000 ", "(2) 10000=%'Rf ", x); + + mpfr_set_ui (x, 1000, MPFR_RNDN); + check_sprintf ("(3) 1000=1,000.000000 ", "(3) 1000=%'Rf ", x); + + for (i = 1; i <= sizeof (v) - 3; i++) + { + char buf[64]; + int j; + + strcpy (buf, "(4) 10^i=1"); + for (j = i; j > 0; j--) + strcat (buf, ",0" + (j % 3 != 0)); + strcat (buf, " "); + mpfr_set_str (x, v + sizeof (v) - 3 - i, 10, MPFR_RNDN); + check_sprintf (buf, "(4) 10^i=%'.0Rf ", x); + } + +#define N0 20 + + for (i = 1; i <= N0; i++) + { + char s[N0+4], buf[64]; + int j; + + s[0] = '1'; + for (j = 1; j <= i; j++) + s[j] = '0'; + s[i+1] = '\0'; + + strcpy (buf, "(5) 10^i=1"); + for (j = i; j > 0; j--) + strcat (buf, ",0" + (j % 3 != 0)); + strcat (buf, " "); + + mpfr_set_str (x, s, 10, MPFR_RNDN); + + check_sprintf (buf, "(5) 10^i=%'.0RNf ", x); + check_sprintf (buf, "(5) 10^i=%'.0RZf ", x); + check_sprintf (buf, "(5) 10^i=%'.0RUf ", x); + check_sprintf (buf, "(5) 10^i=%'.0RDf ", x); + check_sprintf (buf, "(5) 10^i=%'.0RYf ", x); + + strcat (s + (i + 1), ".5"); + check_sprintf (buf, "(5) 10^i=%'.0Rf ", x); + } + + mpfr_set_str (x, "1000", 10, MPFR_RNDN); + check_sprintf ("00000001e+03", "%'012.3Rg", x); + check_sprintf ("00000001,000", "%'012.4Rg", x); + check_sprintf ("000000001,000", "%'013.4Rg", x); + + mpfr_clear (x); +} + +#else + +static void +test_locale (void) +{ + if (getenv ("MPFR_CHECK_LOCALES") != NULL) + { + fprintf (stderr, "Cannot test locales.\n"); + exit (1); + } +} + +#endif + int main (int argc, char **argv) { @@ -1446,12 +1583,14 @@ binary (); decimal (); -#if defined(HAVE_LOCALE_H) && defined(HAVE_SETLOCALE) -#if MPFR_LCONV_DPTS +#if defined(HAVE_LOCALE_H) && defined(HAVE_SETLOCALE) && MPFR_LCONV_DPTS locale_da_DK (); - /* Avoid a warning by doing the setlocale outside of this #if */ -#endif - setlocale (LC_ALL, "C"); +#else + if (getenv ("MPFR_CHECK_LOCALES") != NULL) + { + fprintf (stderr, "Cannot test locales.\n"); + exit (1); + } #endif } @@ -1462,6 +1601,7 @@ snprintf_size (); percent_n (); mixed (); + test_locale (); if (getenv ("MPFR_CHECK_LIBC_PRINTF")) { diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-23 08:30:20.241198970 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-23 08:30:20.269198598 +0000 @@ -0,0 +1 @@ +tconst_pi diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-10 15:46:13.676796949 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-23 08:30:20.269198598 +0000 @@ -1 +1 @@ -4.0.1-p9 +4.0.1-p10 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-10 15:46:13.672796954 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-23 08:30:20.269198598 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p9" +#define MPFR_VERSION_STRING "4.0.1-p10" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-10 15:46:13.672796954 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-23 08:30:20.269198598 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p9"; + return "4.0.1-p10"; } diff -Naurd mpfr-4.0.1-a/tests/tconst_pi.c mpfr-4.0.1-b/tests/tconst_pi.c --- mpfr-4.0.1-a/tests/tconst_pi.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tconst_pi.c 2018-07-23 08:30:20.261198705 +0000 @@ -81,9 +81,11 @@ # define RUN_PTHREAD_TEST() \ (MPFR_ASSERTN(mpfr_buildopt_sharedcache_p() == 1), run_pthread_test()) + #else -# define RUN_PTHREAD_TEST() \ - (MPFR_ASSERTN(mpfr_buildopt_sharedcache_p() == 0)) + +# define RUN_PTHREAD_TEST() ((void) 0) + #endif /* tconst_pi [prec] [rnd] [0 = no print] */ diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-23 08:39:40.951246327 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-23 08:39:40.983246047 +0000 @@ -0,0 +1 @@ +erf diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-23 08:30:20.269198598 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-23 08:39:40.983246047 +0000 @@ -1 +1 @@ -4.0.1-p10 +4.0.1-p11 diff -Naurd mpfr-4.0.1-a/src/erf.c mpfr-4.0.1-b/src/erf.c --- mpfr-4.0.1-a/src/erf.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/erf.c 2018-07-23 08:39:40.971246153 +0000 @@ -186,10 +186,8 @@ { mpfr_prec_t n, m; mpfr_exp_t nuk, sigmak; - double tauk; mpfr_t y, s, t, u; unsigned int k; - int log2tauk; int inex; MPFR_GROUP_DECL (group); MPFR_ZIV_DECL (loop); @@ -204,10 +202,14 @@ MPFR_ZIV_INIT (loop, m); for (;;) { + mpfr_t tauk; + mpfr_exp_t log2tauk; + mpfr_mul (y, x, x, MPFR_RNDU); /* err <= 1 ulp */ mpfr_set_ui (s, 1, MPFR_RNDN); mpfr_set_ui (t, 1, MPFR_RNDN); - tauk = 0.0; + mpfr_init2 (tauk, 53); + mpfr_set_ui (tauk, 0, MPFR_RNDU); for (k = 1; ; k++) { @@ -222,12 +224,13 @@ sigmak -= MPFR_GET_EXP(s); nuk = MPFR_GET_EXP(u) - MPFR_GET_EXP(s); - if ((nuk < - (mpfr_exp_t) m) && ((double) k >= xf2)) + if ((nuk < - (mpfr_exp_t) m) && (k >= xf2)) break; /* tauk <- 1/2 + tauk * 2^sigmak + (1+8k)*2^nuk */ - tauk = 0.5 + mul_2exp (tauk, sigmak) - + mul_2exp (1.0 + 8.0 * (double) k, nuk); + mpfr_mul_2si (tauk, tauk, sigmak, MPFR_RNDU); + mpfr_add_d (tauk, tauk, 0.5 + mul_2exp (1.0 + 8.0 * (double) k, nuk), + MPFR_RNDU); } mpfr_mul (s, x, s, MPFR_RNDU); @@ -236,8 +239,14 @@ mpfr_const_pi (t, MPFR_RNDZ); mpfr_sqrt (t, t, MPFR_RNDZ); mpfr_div (s, s, t, MPFR_RNDN); - tauk = 4.0 * tauk + 11.0; /* final ulp-error on s */ - log2tauk = __gmpfr_ceil_log2 (tauk); + /* tauk = 4 * tauk + 11: final ulp-error on s */ + mpfr_mul_2ui (tauk, tauk, 2, MPFR_RNDU); + mpfr_add_ui (tauk, tauk, 11, MPFR_RNDU); + /* In practice, one should not get an infinity, as it would require + a huge precision and a lot of time, but just in case... */ + MPFR_ASSERTN (!MPFR_IS_INF (tauk)); + log2tauk = MPFR_GET_EXP (tauk); + mpfr_clear (tauk); if (MPFR_LIKELY (MPFR_CAN_ROUND (s, m - log2tauk, n, rnd_mode))) break; diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-23 08:30:20.269198598 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-23 08:39:40.983246047 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p10" +#define MPFR_VERSION_STRING "4.0.1-p11" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-23 08:30:20.269198598 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-23 08:39:40.983246047 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p10"; + return "4.0.1-p11"; } diff -Naurd mpfr-4.0.1-a/tests/terf.c mpfr-4.0.1-b/tests/terf.c --- mpfr-4.0.1-a/tests/terf.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/terf.c 2018-07-23 08:39:40.971246153 +0000 @@ -638,6 +638,20 @@ mpfr_set_emax (emax); } +/* Similar to a bug reported by Naoki Shibata: + https://sympa.inria.fr/sympa/arc/mpfr/2018-07/msg00028.html +*/ +static void +bug20180723 (void) +{ + mpfr_t x; + + mpfr_init2 (x, 256); + mpfr_set_ui (x, 28, MPFR_RNDN); + mpfr_erfc (x, x, MPFR_RNDN); + mpfr_clear (x); +} + int main (int argc, char *argv[]) { @@ -648,9 +662,10 @@ large_arg (); test_erfc (); reduced_expo_range (); + bug20180723 (); - test_generic_erf (MPFR_PREC_MIN, 100, 15); - test_generic_erfc (MPFR_PREC_MIN, 100, 15); + test_generic_erf (MPFR_PREC_MIN, 300, 150); + test_generic_erfc (MPFR_PREC_MIN, 300, 150); data_check ("data/erf", mpfr_erf, "mpfr_erf"); data_check ("data/erfc", mpfr_erfc, "mpfr_erfc"); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-31 11:21:08.747432173 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-31 11:21:08.887430085 +0000 @@ -0,0 +1 @@ +vasprintf-error-ub diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-23 08:39:40.983246047 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-31 11:21:08.887430085 +0000 @@ -1 +1 @@ -4.0.1-p11 +4.0.1-p12 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-23 08:39:40.983246047 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-31 11:21:08.887430085 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p11" +#define MPFR_VERSION_STRING "4.0.1-p12" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/vasprintf.c mpfr-4.0.1-b/src/vasprintf.c --- mpfr-4.0.1-a/src/vasprintf.c 2018-07-10 15:46:13.648796978 +0000 +++ mpfr-4.0.1-b/src/vasprintf.c 2018-07-31 11:21:08.771431816 +0000 @@ -33,11 +33,22 @@ So, for the time being, we return a negative value and set the erange flag, and set errno to EOVERFLOW in POSIX system. */ -/* Note: Due to limitations from the C standard and GMP, if - size_t < unsigned int (which is allowed by the C standard but unlikely - to occur on any platform), the behavior is undefined for output that - would reach SIZE_MAX = (size_t) -1 (if the result cannot be delivered, - there should be an assertion failure, but this could not be tested). */ +/* Notes about limitations on some platforms: + + Due to limitations from the C standard and GMP, if size_t < unsigned int + (which is allowed by the C standard but unlikely to occur on any + platform), the behavior is undefined for output that would reach + SIZE_MAX = (size_t) -1 (if the result cannot be delivered, there should + be an assertion failure, but this could not be tested). + + The stdarg(3) Linux man page says: + On some systems, va_end contains a closing '}' matching a '{' in + va_start, so that both macros must occur in the same function, + and in a way that allows this. + However, the only requirement from ISO C is that both macros must be + invoked in the same function (MPFR uses va_copy instead of va_start, + but the requirement is the same). Here, MPFR just follows ISO C. +*/ #ifdef HAVE_CONFIG_H # include "config.h" @@ -2273,44 +2284,42 @@ va_end (ap2); - if (buf.len > INT_MAX) /* overflow */ - buf.len = -1; + if (buf.len == -1 || buf.len > INT_MAX) /* overflow */ + goto overflow; - if (buf.len != -1) - { - nbchar = buf.len; - MPFR_ASSERTD (nbchar >= 0); + nbchar = buf.len; + MPFR_ASSERTD (nbchar >= 0); - if (ptr != NULL) /* implement mpfr_vasprintf */ + if (ptr != NULL) /* implement mpfr_vasprintf */ + { + MPFR_ASSERTD (nbchar == strlen (buf.start)); + *ptr = (char *) mpfr_reallocate_func (buf.start, buf.size, nbchar + 1); + } + else if (size != 0) /* implement mpfr_vsnprintf */ + { + if (nbchar < size) { - MPFR_ASSERTD (nbchar == strlen (buf.start)); - *ptr = (char *) - mpfr_reallocate_func (buf.start, buf.size, nbchar + 1); + strncpy (Buf, buf.start, nbchar); + Buf[nbchar] = '\0'; } - else if (size > 0) /* implement mpfr_vsnprintf */ + else { - if (nbchar < size) - { - strncpy (Buf, buf.start, nbchar); - Buf[nbchar] = '\0'; - } - else - { - strncpy (Buf, buf.start, size - 1); - Buf[size-1] = '\0'; - } - mpfr_free_func (buf.start, buf.size); + strncpy (Buf, buf.start, size - 1); + Buf[size-1] = '\0'; } - - MPFR_SAVE_EXPO_FREE (expo); - return nbchar; /* return the number of characters that would have - been written had 'size' been sufficiently large, - not counting the terminating null character */ + mpfr_free_func (buf.start, buf.size); } + MPFR_SAVE_EXPO_FREE (expo); + return nbchar; /* return the number of characters that would have + been written had 'size' been sufficiently large, + not counting the terminating null character */ + error: + va_end (ap2); if (buf.len == -1) /* overflow */ { + overflow: MPFR_LOG_MSG (("Overflow\n", 0)); MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_ERANGE); #ifdef EOVERFLOW @@ -2320,8 +2329,10 @@ } MPFR_SAVE_EXPO_FREE (expo); - *ptr = NULL; - mpfr_free_func (buf.start, buf.size); + if (ptr != NULL) /* implement mpfr_vasprintf */ + *ptr = NULL; + if (ptr != NULL || size != 0) + mpfr_free_func (buf.start, buf.size); return -1; } diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-23 08:39:40.983246047 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-31 11:21:08.887430085 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p11"; + return "4.0.1-p12"; } diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-07-31 11:26:25.226690969 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-07-31 11:26:25.262690428 +0000 @@ -0,0 +1 @@ +vasprintf-p-length-modifier diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-31 11:21:08.887430085 +0000 +++ mpfr-4.0.1-b/VERSION 2018-07-31 11:26:25.262690428 +0000 @@ -1 +1 @@ -4.0.1-p12 +4.0.1-p13 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-31 11:21:08.887430085 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-07-31 11:26:25.262690428 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p12" +#define MPFR_VERSION_STRING "4.0.1-p13" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/vasprintf.c mpfr-4.0.1-b/src/vasprintf.c --- mpfr-4.0.1-a/src/vasprintf.c 2018-07-31 11:21:08.771431816 +0000 +++ mpfr-4.0.1-b/src/vasprintf.c 2018-07-31 11:26:25.246690670 +0000 @@ -259,6 +259,8 @@ } } +/* Note: additional flags should be added to the MPFR_PREC_ARG code + for gmp_asprintf (when supported). */ static const char * parse_flags (const char *format, struct printf_spec *specinfo) { @@ -2220,7 +2222,7 @@ /* output mpfr_prec_t variable */ { char *s; - char format[MPFR_PREC_FORMAT_SIZE + 6]; /* see examples below */ + char format[MPFR_PREC_FORMAT_SIZE + 12]; /* e.g. "%0#+ -'*.*ld\0" */ size_t length; mpfr_prec_t prec; @@ -2232,14 +2234,14 @@ start = fmt; /* construct format string, like "%*.*hd" "%*.*d" or "%*.*ld" */ - format[0] = '%'; - format[1] = '*'; - format[2] = '.'; - format[3] = '*'; - format[4] = '\0'; - strcat (format, MPFR_PREC_FORMAT_TYPE); - format[4 + MPFR_PREC_FORMAT_SIZE] = spec.spec; - format[5 + MPFR_PREC_FORMAT_SIZE] = '\0'; + sprintf (format, "%%%s%s%s%s%s%s*.*" MPFR_PREC_FORMAT_TYPE "%c", + spec.pad == '0' ? "0" : "", + spec.alt ? "#" : "", + spec.showsign ? "+" : "", + spec.space ? " " : "", + spec.left ? "-" : "", + spec.group ? "'" : "", + spec.spec); length = gmp_asprintf (&s, format, spec.width, spec.prec, prec); MPFR_ASSERTN (length >= 0); /* guaranteed by GMP 6 */ buffer_cat (&buf, s, length); diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-31 11:21:08.887430085 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-07-31 11:26:25.262690428 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p12"; + return "4.0.1-p13"; } diff -Naurd mpfr-4.0.1-a/tests/tsprintf.c mpfr-4.0.1-b/tests/tsprintf.c --- mpfr-4.0.1-a/tests/tsprintf.c 2018-07-10 15:46:13.648796978 +0000 +++ mpfr-4.0.1-b/tests/tsprintf.c 2018-07-31 11:26:25.250690608 +0000 @@ -232,6 +232,22 @@ /* specifier 'P' for precision */ check_vsprintf ("128", "%Pu", p); check_vsprintf ("00128", "%.5Pu", p); + check_vsprintf (" 128", "%5Pu", p); + check_vsprintf ("000128", "%06Pu", p); + check_vsprintf ("128 :", "%-7Pu:", p); + check_vsprintf ("000128:", "%-2.6Pd:", p); + check_vsprintf (" 000128:", "%8.6Pd:", p); + check_vsprintf ("000128 :", "%-8.6Pd:", p); + check_vsprintf ("+128:", "%+d:", p); + check_vsprintf (" 128:", "% d:", p); + check_vsprintf ("80:", "% x:", p); + check_vsprintf ("0x80:", "% #x:", p); + check_vsprintf ("0x80:", "%0#+ -x:", p); + check_vsprintf ("0200:", "%0#+ -o:", p); + check_vsprintf ("+0000128 :", "%0+ *.*Pd:", -9, 7, p); + check_vsprintf ("+12345 :", "%0+ -*.*Pd:", -9, -3, (mpfr_prec_t) 12345); + /* Do not add a test like "%05.1Pd" as MS Windows is buggy: when + a precision is given, the '0' flag must be ignored. */ /* special numbers */ mpfr_set_inf (x, 1); diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2018-12-31 11:18:05.589387869 +0000 +++ mpfr-4.0.1-b/PATCHES 2018-12-31 11:18:05.621387341 +0000 @@ -0,0 +1 @@ +set_1_2 diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-07-31 11:26:25.262690428 +0000 +++ mpfr-4.0.1-b/VERSION 2018-12-31 11:18:05.621387341 +0000 @@ -1 +1 @@ -4.0.1-p13 +4.0.1-p14 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-07-31 11:26:25.262690428 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2018-12-31 11:18:05.617387407 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p13" +#define MPFR_VERSION_STRING "4.0.1-p14" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/set.c mpfr-4.0.1-b/src/set.c --- mpfr-4.0.1-a/src/set.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/set.c 2018-12-31 11:18:05.609387539 +0000 @@ -80,15 +80,20 @@ return mpfr_set4 (a, b, rnd_mode, MPFR_SIGN_POS); } -/* Round (u, inex) into s with rounding mode rnd, where inex is the ternary - value associated to u with the *same* rounding mode. +/* Round (u, inex) into s with rounding mode rnd_mode, where inex is + the ternary value associated with u, which has been obtained using + the *same* rounding mode rnd_mode. Assumes PREC(u) = 2*PREC(s). - The main algorithm is the following: - rnd=RNDZ: inex2 = mpfr_set (s, u, rnd_mode); return inex2 | inex; - (a negative value, if any, is preserved in inex2 | inex) - rnd=RNDA: idem - rnd=RNDN: inex2 = mpfr_set (s, u, rnd_mode); - if (inex2) return inex2; else return inex; */ + The generic algorithm is the following: + 1. inex2 = mpfr_set (s, u, rnd_mode); + 2. if (inex2 != 0) return inex2; else return inex; + except in the double-rounding case: in MPFR_RNDN, when u is in the + middle of two consecutive PREC(s)-bit numbers, if inex and inex2 + are both > 0 (resp. both < 0), we correct s to nextbelow(s) (resp. + nextabove(s)), and return the opposite of inex. + Note: this function can be called with rnd_mode == MPFR_RNDF, in + which case, the rounding direction and the returned ternary value + are unspecified. */ int mpfr_set_1_2 (mpfr_ptr s, mpfr_srcptr u, mpfr_rnd_t rnd_mode, int inex) { @@ -106,9 +111,9 @@ return inex; } - MPFR_ASSERTD(MPFR_PREC(u) == 2 * MPFR_PREC(s)); + MPFR_ASSERTD(MPFR_PREC(u) == 2 * p); - if (MPFR_PREC(s) < GMP_NUMB_BITS) + if (p < GMP_NUMB_BITS) { mask = MPFR_LIMB_MASK(sh); @@ -194,6 +199,19 @@ /* general case PREC(s) >= GMP_NUMB_BITS */ inex2 = mpfr_set (s, u, rnd_mode); - return (rnd_mode != MPFR_RNDN) ? inex | inex2 - : (inex2) ? inex2 : inex; + /* Check the double-rounding case, i.e. with u = middle of two + consecutive PREC(s)-bit numbers, which is equivalent to u being + exactly representable on PREC(s) + 1 bits but not on PREC(s) bits. + Moreover, since PREC(u) = 2*PREC(s), u and s cannot be identical + (as pointers), thus u was not changed. */ + if (rnd_mode == MPFR_RNDN && inex * inex2 > 0 && + mpfr_min_prec (u) == p + 1) + { + if (inex > 0) + mpfr_nextbelow (s); + else + mpfr_nextabove (s); + return -inex; + } + return inex2 != 0 ? inex2 : inex; } diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-07-31 11:26:25.262690428 +0000 +++ mpfr-4.0.1-b/src/version.c 2018-12-31 11:18:05.621387341 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p13"; + return "4.0.1-p14"; } diff -Naurd mpfr-4.0.1-a/tests/tfmma.c mpfr-4.0.1-b/tests/tfmma.c --- mpfr-4.0.1-a/tests/tfmma.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tfmma.c 2018-12-31 11:18:05.609387539 +0000 @@ -469,6 +469,69 @@ mpfr_clears (x, y, z, (mpfr_ptr) 0); } +/* 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. */ +static void +double_rounding (void) +{ + int p; + + for (p = 4; p < 4 * GMP_NUMB_BITS; p++) + { + mpfr_t a, b, c, d, z1, z2; + int neg; + + mpfr_inits2 (p, a, b, c, d, z1, z2, (mpfr_ptr) 0); + /* Choose a, b, c, d such that a*b = 2^p and c*d = 1 + epsilon */ + mpfr_set_ui (a, 2, MPFR_RNDN); + mpfr_set_ui_2exp (b, 1, p - 1, MPFR_RNDN); + mpfr_set_ui (c, 1, MPFR_RNDN); + mpfr_nextabove (c); + /* c = 1 + ulp(1) = 1 + 2 * ulp(1/2) */ + + for (neg = 0; neg <= 1; neg++) + { + int inex1, inex2; + + /* Set d = 1 - (1 + neg) * ulp(1/2), thus + * c*d = 1 + (1 - neg) * ulp(1/2) - 2 * (1 + neg) * ulp(1/2)^2, + * so that a*b + c*d rounds to 2^p + 1 and epsilon has the + * same sign as (-1)^neg. + */ + mpfr_set_ui (d, 1, MPFR_RNDN); + mpfr_nextbelow (d); + mpfr_set_ui_2exp (z1, 1, p, MPFR_RNDN); + if (neg) + { + mpfr_nextbelow (d); + inex1 = -1; + } + else + { + mpfr_nextabove (z1); + inex1 = 1; + } + + inex2 = mpfr_fmma (z2, a, b, c, d, MPFR_RNDN); + if (!(mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2))) + { + printf ("Error in double_rounding for precision %d, neg=%d\n", + p, neg); + printf ("Expected "); + mpfr_dump (z1); + printf ("with inex = %d\n", inex1); + printf ("Got "); + mpfr_dump (z2); + printf ("with inex = %d\n", inex2); + exit (1); + } + } + + mpfr_clears (a, b, c, d, z1, z2, (mpfr_ptr) 0); + } +} + int main (int argc, char *argv[]) { @@ -480,6 +543,7 @@ overflow_tests (); half_plus_half (); bug20170405 (); + double_rounding (); tests_end_mpfr (); return 0; diff -Naurd mpfr-4.0.1-a/tests/tset.c mpfr-4.0.1-b/tests/tset.c --- mpfr-4.0.1-a/tests/tset.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tset.c 2018-12-31 11:18:05.609387539 +0000 @@ -256,6 +256,93 @@ mpfr_clear (y); } +static void +test_set_1_2 (void) +{ + mpfr_t u, v, zz, z; + int inex; + + /* (8,16)-bit test */ + mpfr_inits2 (16, u, v, zz, (mpfr_ptr) 0); + mpfr_init2 (z, 8); + mpfr_set_str_binary (u, "0.1100001100011010E-1"); + mpfr_set_str_binary (v, "0.1100010101110010E0"); + /* u + v = 1.0010011011111111 */ + inex = mpfr_add (zz, u, v, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + mpfr_set_str_binary (u, "1.001001110000000"); + MPFR_ASSERTN(mpfr_equal_p (zz, u)); + inex = mpfr_set_1_2 (z, zz, MPFR_RNDN, inex); + /* we should have z = 1.0010011 and inex < 0 */ + MPFR_ASSERTN(inex < 0); + mpfr_set_str_binary (u, "1.0010011"); + MPFR_ASSERTN(mpfr_equal_p (z, u)); + mpfr_clears (u, v, zz, z, (mpfr_ptr) 0); + + /* (16,32)-bit test: + * take for v a random 32-bit number in [1/2,1), here 2859611790/2^32 + * take for z a random 16-bit number in [1,2), less than 2*v, + with last bit 0, here we take z = 40900/2^15 + * take u = z-v-1/2^16-1/2^32 */ + mpfr_inits2 (32, u, v, zz, (mpfr_ptr) 0); + mpfr_init2 (z, 16); + mpfr_set_str_binary (u, "0.10010101000101001100100101110001"); + mpfr_set_str_binary (v, "0.10101010011100100011011010001110"); + /* u + v = 1.00111111100001101111111111111111 */ + inex = mpfr_add (zz, u, v, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + mpfr_set_str_binary (u, "1.0011111110000111"); + MPFR_ASSERTN(mpfr_equal_p (zz, u)); + inex = mpfr_set_1_2 (z, zz, MPFR_RNDN, inex); + /* we should have z = 1.001111111000011 and inex < 0 */ + MPFR_ASSERTN(inex < 0); + mpfr_set_str_binary (u, "1.001111111000011"); + MPFR_ASSERTN(mpfr_equal_p (z, u)); + mpfr_clears (u, v, zz, z, (mpfr_ptr) 0); + + /* (32,64)-bit test: + * take for v a random 64-bit number in [1/2,1), + here v = 13687985014345662879/2^64 + * take for z a random 32-bit number in [1,2), less than 2*v, + with last bit 0, here we take z = 2871078774/2^31 + * take u = z-v-1/2^32-1/2^64 */ + mpfr_inits2 (64, u, v, zz, (mpfr_ptr) 0); + mpfr_init2 (z, 32); + mpfr_set_str_binary (u, "0.10011000010011001110000100010001110010010000111001111110011"); + mpfr_set_str_binary (v, "0.1011110111110101011111011101100100110110111100011000000110011111"); + /* u + v = 1.0101011001000010010111101110101011111111111111111111111111111111 */ + inex = mpfr_add (zz, u, v, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + mpfr_set_str_binary (u, "1.01010110010000100101111011101011"); + MPFR_ASSERTN(mpfr_equal_p (zz, u)); + inex = mpfr_set_1_2 (z, zz, MPFR_RNDN, inex); + /* we should have z = 1.0101011001000010010111101110101 and inex < 0 */ + MPFR_ASSERTN(inex < 0); + mpfr_set_str_binary (u, "1.0101011001000010010111101110101"); + MPFR_ASSERTN(mpfr_equal_p (z, u)); + mpfr_clears (u, v, zz, z, (mpfr_ptr) 0); + + /* (64,128)-bit test: + * take for v a random 128-bit number in [1/2,1), + here v = 322263811942091240216761391118876232409/2^128 + * take for z a random 64-bit number in [1,2), less than 2*v, + with last bit 0, here we take z = 16440347967874738276/2^63 + * take u = z-v-1/2^64-1/2^128 */ + mpfr_inits2 (128, u, v, zz, (mpfr_ptr) 0); + mpfr_init2 (z, 64); + mpfr_set_str_binary (u, "0.1101010111011101111100100001011111111000010011011001000101111010110101101101011011100110101001010001101011011110101101010010011"); + mpfr_set_str_binary (v, "0.11110010011100011100000010100110100010011010110010111111010011000010100100101001000110010101101011100101001000010100101011011001"); + inex = mpfr_add (zz, u, v, MPFR_RNDN); + MPFR_ASSERTN(inex > 0); + mpfr_set_str_binary (u, "1.1100100001001111101100101011111010000001111110100101000011000111"); + MPFR_ASSERTN(mpfr_equal_p (zz, u)); + inex = mpfr_set_1_2 (z, zz, MPFR_RNDN, inex); + MPFR_ASSERTN(inex < 0); + mpfr_set_str_binary (u, "1.1100100001001111101100101011111010000001111110100101000011000110"); + MPFR_ASSERTN(mpfr_equal_p (z, u)); + mpfr_clears (u, v, zz, z, (mpfr_ptr) 0); +} + #define TEST_FUNCTION mpfr_set #include "tgeneric.c" @@ -268,6 +355,8 @@ tests_start_mpfr (); + test_set_1_2 (); + /* Default : no error */ error = 0; diff -Naurd mpfr-4.0.1-a/PATCHES mpfr-4.0.1-b/PATCHES --- mpfr-4.0.1-a/PATCHES 2019-01-31 12:05:55.995705752 +0000 +++ mpfr-4.0.1-b/PATCHES 2019-01-31 12:05:56.091705514 +0000 @@ -0,0 +1 @@ +strtofr diff -Naurd mpfr-4.0.1-a/VERSION mpfr-4.0.1-b/VERSION --- mpfr-4.0.1-a/VERSION 2018-12-31 11:18:05.621387341 +0000 +++ mpfr-4.0.1-b/VERSION 2019-01-31 12:05:56.091705514 +0000 @@ -1 +1 @@ -4.0.1-p14 +4.0.1-p15 diff -Naurd mpfr-4.0.1-a/src/mpfr.h mpfr-4.0.1-b/src/mpfr.h --- mpfr-4.0.1-a/src/mpfr.h 2018-12-31 11:18:05.617387407 +0000 +++ mpfr-4.0.1-b/src/mpfr.h 2019-01-31 12:05:56.079705544 +0000 @@ -27,7 +27,7 @@ #define MPFR_VERSION_MAJOR 4 #define MPFR_VERSION_MINOR 0 #define MPFR_VERSION_PATCHLEVEL 1 -#define MPFR_VERSION_STRING "4.0.1-p14" +#define MPFR_VERSION_STRING "4.0.1-p15" /* User macros: MPFR_USE_FILE: Define it to make MPFR define functions dealing diff -Naurd mpfr-4.0.1-a/src/strtofr.c mpfr-4.0.1-b/src/strtofr.c --- mpfr-4.0.1-a/src/strtofr.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/src/strtofr.c 2019-01-31 12:05:56.051705613 +0000 @@ -36,79 +36,80 @@ allocated for the mantissa field. */ size_t prec; /* length of mant (zero for +/-0) */ size_t alloc; /* allocation size of mantissa */ - mpfr_exp_t exp_base; /* number of digits before the point */ - mpfr_exp_t exp_bin; /* exponent in case base=2 or 16, and the pxxx - format is used (i.e., exponent is given in - base 10) */ + mpfr_exp_t exp_base; /* number of digits before the point, + exponent + except in case of binary exponent (exp_bin) */ + mpfr_exp_t exp_bin; /* binary exponent of the pxxx format for + base = 2 or 16 */ }; /* This table has been generated by the following program. For 2 <= b <= MPFR_MAX_BASE, RedInvLog2Table[b-2][0] / RedInvLog2Table[b-2][1] - is an upper approximation of log(2)/log(b). + is an upper approximation to log(2)/log(b), no larger than 1. + Note: these numbers must fit on 16 bits, thus unsigned int is OK. */ -static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = { - {1UL, 1UL}, - {53UL, 84UL}, - {1UL, 2UL}, - {4004UL, 9297UL}, - {53UL, 137UL}, - {2393UL, 6718UL}, - {1UL, 3UL}, - {665UL, 2108UL}, - {4004UL, 13301UL}, - {949UL, 3283UL}, - {53UL, 190UL}, - {5231UL, 19357UL}, - {2393UL, 9111UL}, - {247UL, 965UL}, - {1UL, 4UL}, - {4036UL, 16497UL}, - {665UL, 2773UL}, - {5187UL, 22034UL}, - {4004UL, 17305UL}, - {51UL, 224UL}, - {949UL, 4232UL}, - {3077UL, 13919UL}, - {53UL, 243UL}, - {73UL, 339UL}, - {5231UL, 24588UL}, - {665UL, 3162UL}, - {2393UL, 11504UL}, - {4943UL, 24013UL}, - {247UL, 1212UL}, - {3515UL, 17414UL}, - {1UL, 5UL}, - {4415UL, 22271UL}, - {4036UL, 20533UL}, - {263UL, 1349UL}, - {665UL, 3438UL}, - {1079UL, 5621UL}, - {5187UL, 27221UL}, - {2288UL, 12093UL}, - {4004UL, 21309UL}, - {179UL, 959UL}, - {51UL, 275UL}, - {495UL, 2686UL}, - {949UL, 5181UL}, - {3621UL, 19886UL}, - {3077UL, 16996UL}, - {229UL, 1272UL}, - {53UL, 296UL}, - {109UL, 612UL}, - {73UL, 412UL}, - {1505UL, 8537UL}, - {5231UL, 29819UL}, - {283UL, 1621UL}, - {665UL, 3827UL}, - {32UL, 185UL}, - {2393UL, 13897UL}, - {1879UL, 10960UL}, - {4943UL, 28956UL}, - {409UL, 2406UL}, - {247UL, 1459UL}, - {231UL, 1370UL}, - {3515UL, 20929UL} }; +static const unsigned int RedInvLog2Table[MPFR_MAX_BASE-1][2] = { + {1, 1}, + {53, 84}, + {1, 2}, + {4004, 9297}, + {53, 137}, + {2393, 6718}, + {1, 3}, + {665, 2108}, + {4004, 13301}, + {949, 3283}, + {53, 190}, + {5231, 19357}, + {2393, 9111}, + {247, 965}, + {1, 4}, + {4036, 16497}, + {665, 2773}, + {5187, 22034}, + {4004, 17305}, + {51, 224}, + {949, 4232}, + {3077, 13919}, + {53, 243}, + {73, 339}, + {5231, 24588}, + {665, 3162}, + {2393, 11504}, + {4943, 24013}, + {247, 1212}, + {3515, 17414}, + {1, 5}, + {4415, 22271}, + {4036, 20533}, + {263, 1349}, + {665, 3438}, + {1079, 5621}, + {5187, 27221}, + {2288, 12093}, + {4004, 21309}, + {179, 959}, + {51, 275}, + {495, 2686}, + {949, 5181}, + {3621, 19886}, + {3077, 16996}, + {229, 1272}, + {53, 296}, + {109, 612}, + {73, 412}, + {1505, 8537}, + {5231, 29819}, + {283, 1621}, + {665, 3827}, + {32, 185}, + {2393, 13897}, + {1879, 10960}, + {4943, 28956}, + {409, 2406}, + {247, 1459}, + {231, 1370}, + {3515, 20929} }; #if 0 #define N 8 int main () @@ -377,6 +378,13 @@ res = 1; MPFR_ASSERTD (pstr->exp_base >= 0); + /* FIXME: In the code below (both cases), if the exponent from the + string is large, it will be replaced by MPFR_EXP_MIN or MPFR_EXP_MAX, + i.e. it will have a different value. This may not change the result + in most cases, but there is no guarantee on very long strings when + mpfr_exp_t is a 32-bit type, as the exponent could be brought back + to the current exponent range. */ + /* an optional exponent (e or E, p or P, @) */ if ( (*str == '@' || (base <= 10 && (*str == 'e' || *str == 'E'))) && (!isspace((unsigned char) str[1])) ) @@ -445,112 +453,174 @@ static int parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd) { - mpfr_prec_t prec; + mpfr_prec_t precx, prec, ysize_bits, pstr_size; mpfr_exp_t exp; - mpfr_exp_t ysize_bits; - mp_limb_t *y, *result; + mp_limb_t *result; int count, exact; - size_t pstr_size; - mp_size_t ysize, real_ysize; + mp_size_t ysize, real_ysize, diff_ysize; int res, err; + const int extra_limbs = GMP_NUMB_BITS >= 12 ? 1 : 2; /* see below */ MPFR_ZIV_DECL (loop); MPFR_TMP_DECL (marker); /* initialize the working precision */ - prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); + precx = MPFR_GET_PREC (x); + prec = precx + MPFR_INT_CEIL_LOG2 (precx); - /* compute the value y of the leading characters as long as rounding is not - possible */ + /* Compute the value y of the leading characters as long as rounding is not + possible. + Note: We have some integer overflow checking using MPFR_EXP_MIN and + MPFR_EXP_MAX in this loop. Thanks to the large margin between these + extremal values of the mpfr_exp_t type and the valid minimum/maximum + exponents, such integer overflows would correspond to real underflow + or overflow on the result (possibly except in huge precisions, which + are disregarded here; anyway, in practice, such issues could occur + only with 32-bit precision and exponent types). Such checks could be + extended to real early underflow/overflow checking, in order to avoid + useless computations in such cases; in such a case, be careful that + the approximation errors need to be taken into account. */ MPFR_TMP_MARK(marker); MPFR_ZIV_INIT (loop, prec); for (;;) { - /* Set y to the value of the ~prec most significant bits of pstr->mant - (as long as we guarantee correct rounding, we don't need to get - exactly prec bits). */ + mp_limb_t *y0, *y; + + /* y will be regarded as a number with precision prec. */ ysize = MPFR_PREC2LIMBS (prec); /* prec bits corresponds to ysize limbs */ - ysize_bits = ysize * GMP_NUMB_BITS; - /* and to ysize_bits >= prec > MPFR_PREC (x) bits */ - /* we need to allocate one more limb to work around bug + ysize_bits = (mpfr_prec_t) ysize * GMP_NUMB_BITS; + MPFR_ASSERTD (ysize_bits >= prec); + /* and to ysize_bits >= prec > precx bits. */ + /* We need to allocate one more limb as specified by mpn_set_str + (a limb may be written in rp[rn]). Note that the manual of GMP + up to 5.1.3 was incorrect on this point. + See the following discussion: https://gmplib.org/list-archives/gmp-bugs/2013-December/003267.html */ - y = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 2); - y += ysize; /* y has (ysize+2) allocated limbs */ + y0 = MPFR_TMP_LIMBS_ALLOC (2 * ysize + extra_limbs + 1); + y = y0 + ysize; /* y has (ysize + extra_limbs + 1) allocated limbs */ - /* pstr_size is the number of characters we read in pstr->mant - to have at least ysize full limbs. + /* pstr_size is the number of bytes we want to read from pstr->mant + to fill at least ysize full limbs with mpn_set_str. We must have base^(pstr_size-1) >= (2^(GMP_NUMB_BITS))^ysize (in the worst case, the first digit is one and all others are zero). i.e., pstr_size >= 1 + ysize*GMP_NUMB_BITS/log2(base) Since ysize ~ prec/GMP_NUMB_BITS and prec < Umax/2 => ysize*GMP_NUMB_BITS can not overflow. We compute pstr_size = 1 + ceil(ysize_bits * Num / Den) - where Num/Den >= 1/log2(base) - It is not exactly ceil(1/log2(base)) but could be one more (base 2) + where 1/log2(base) <= Num/Den <= 1 + It is not exactly ceil(1/log2(base)) but could be one more (base 2). Quite ugly since it tries to avoid overflow: let Num = RedInvLog2Table[pstr->base-2][0] and Den = RedInvLog2Table[pstr->base-2][1], and ysize_bits = a*Den+b, then ysize_bits * Num/Den = a*Num + (b * Num)/Den, thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den + + Note: denoting m = pstr_size and n = ysize_bits, assuming we have + m = 1 + ceil(n/log2(b)), i.e., b^(m-1) >= 2^n > b^(m-2), then + b^(m-1)/2^n < b, and since we consider m characters of the input, + the corresponding part is less than b^m < b^2*2^n. + This implies that if b^2 < 2^GMP_NUMB_BITS, which for b <= 62 holds + for GMP_NUMB_BITS >= 12, we have real_ysize <= ysize+1 below + (this also implies that for GMP_NUMB_BITS >= 13, the number of bits + of y[real_ysize-1] below is less than GMP_NUMB_BITS, thus + count < GMP_NUMB_BITS). + Warning: for GMP_NUMB_BITS=8, we can have real_ysize = ysize + 2! + Hence the allocation above for ysize + extra_limbs limbs. */ { - unsigned long Num = RedInvLog2Table[pstr->base-2][0]; - unsigned long Den = RedInvLog2Table[pstr->base-2][1]; - pstr_size = ((ysize_bits / Den) * Num) - + (((ysize_bits % Den) * Num + Den - 1) / Den) + unsigned int Num = RedInvLog2Table[pstr->base-2][0]; + unsigned int Den = RedInvLog2Table[pstr->base-2][1]; + MPFR_ASSERTD (Num <= Den && Den <= 65535); /* thus no overflow */ + pstr_size = (ysize_bits / Den) * Num + + ((unsigned long) (ysize_bits % Den) * Num + Den - 1) / Den + 1; } - /* since pstr_size corresponds to at least ysize_bits full bits, - and ysize_bits > prec, the weight of the neglected part of - pstr->mant (if any) is < ulp(y) < ulp(x) */ + /* Since pstr_size corresponds to at least ysize_bits bits, + and ysize_bits >= prec, the weight of the neglected part of + pstr->mant (if any) is < ulp(y) < ulp(x). */ - /* if the number of wanted characters is more than what we have in - pstr->mant, round it down */ - if (pstr_size >= pstr->prec) + /* If the number of wanted bytes is more than what is available + in pstr->mant, i.e. pstr->prec, reduce it to pstr->prec. */ + if (pstr_size > pstr->prec) pstr_size = pstr->prec; - MPFR_ASSERTD (pstr_size == (mpfr_exp_t) pstr_size); - /* convert str into binary: note that pstr->mant is big endian, - thus no offset is needed */ + /* Convert str (potentially truncated to pstr_size) into binary. + Note that pstr->mant is big endian, thus no offset is needed. */ real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base); - MPFR_ASSERTD (real_ysize <= ysize+1); - /* normalize y: warning we can even get ysize+1 limbs! */ + /* See above for the explanation of the following assertion. */ + MPFR_ASSERTD (real_ysize <= ysize + extra_limbs); + + /* The Boolean "exact" will attempt to track exactness of the result: + If it is true, then this means that the result is exact, allowing + termination, even though the rounding test may not succeed. + Conversely, if the result is exact, then "exact" will not + necessarily be true at the end of the Ziv loop, but we will need + to make sure that at some point, "exact" will be true in order to + guarantee termination. FIXME: check that. */ + /* First, consider the part of the input string that has been ignored. + Note that the trailing zeros have been removed in parse_string, so + that if something has been ignored, it must be non-zero. */ + exact = pstr_size == pstr->prec; + + /* Normalize y and set the initial value of its exponent exp, which + is 0 when y is not shifted. + Since pstr->mant was normalized, mpn_set_str guarantees that + the most significant limb is non-zero. */ MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ count_leading_zeros (count, y[real_ysize - 1]); - /* exact means that the number of limbs of the output of mpn_set_str - is less or equal to ysize */ - exact = real_ysize <= ysize; - if (exact) /* shift y to the left in that case y should be exact */ + diff_ysize = ysize - real_ysize; + MPFR_LOG_MSG (("diff_ysize = %ld\n", (long) diff_ysize)); + if (diff_ysize >= 0) { - /* we have enough limbs to store {y, real_ysize} */ - /* shift {y, num_limb} for count bits to the left */ + /* We have enough limbs to store {y, real_ysize} exactly + in {y, ysize}, so that we can do a left shift, without + losing any information ("exact" will not change). */ if (count != 0) - mpn_lshift (y + ysize - real_ysize, y, real_ysize, count); - if (real_ysize != ysize) + mpn_lshift (y + diff_ysize, y, real_ysize, count); + if (diff_ysize > 0) { if (count == 0) - mpn_copyd (y + ysize - real_ysize, y, real_ysize); - MPN_ZERO (y, ysize - real_ysize); + mpn_copyd (y + diff_ysize, y, real_ysize); + MPN_ZERO (y, diff_ysize); } - /* for each bit shift decrease exponent of y */ - /* (This should not overflow) */ - exp = - ((ysize - real_ysize) * GMP_NUMB_BITS + count); + /* exp = negation of the total shift count, avoiding overflows. */ + exp = - ((mpfr_exp_t) diff_ysize * GMP_NUMB_BITS + count); } - else /* shift y to the right, by doing this we might lose some - bits from the result of mpn_set_str (in addition to the - characters neglected from pstr->mant) */ + else { - /* shift {y, num_limb} for (GMP_NUMB_BITS - count) bits - to the right. FIXME: can we prove that count cannot be zero here, - since mpn_rshift does not accept a shift of GMP_NUMB_BITS? */ - MPFR_ASSERTD (count != 0); - exact = mpn_rshift (y, y, real_ysize, GMP_NUMB_BITS - count) == - MPFR_LIMB_ZERO; - /* for each bit shift increase exponent of y */ - exp = GMP_NUMB_BITS - count; + /* Shift {y, real_ysize} for (GMP_NUMB_BITS - count) bits to the + right, and put the ysize most significant limbs into {y, ysize}. + We have either real_ysize = ysize + 1 or real_ysize = ysize + 2 + (only possible with extra_limbs == 2). */ + MPFR_ASSERTD (diff_ysize == -1 || + (extra_limbs == 2 && diff_ysize == -2)); + if (count != 0) + { + /* Before doing the shift, consider the limb that will entirely + be lost if real_ysize = ysize + 2. */ + exact = exact && (diff_ysize == -1 || y[0] == MPFR_LIMB_ZERO); + /* mpn_rshift allows overlap, provided destination <= source */ + /* FIXME: The bits lost due to mpn_rshift are not taken + into account in the error analysis below! */ + if (mpn_rshift (y, y - (diff_ysize + 1), real_ysize, + GMP_NUMB_BITS - count) != MPFR_LIMB_ZERO) + exact = 0; /* some non-zero bits have been shifted out */ + } + else + { + /* the case real_ysize = ysize + 2 with count = 0 cannot happen + even with GMP_NUMB_BITS = 8 since 62^2 < 256^2/2 */ + MPFR_ASSERTD (diff_ysize == -1); + exact = exact && y[0] == MPFR_LIMB_ZERO; + /* copy {y+real_ysize-ysize, ysize} to {y, ysize} */ + mpn_copyi (y, y + 1, real_ysize - 1); + } + /* exp = shift count */ + /* TODO: add some explanations about what exp means exactly. */ + exp = GMP_NUMB_BITS * (- diff_ysize) - count; } /* compute base^(exp_base - pstr_size) on n limbs */ @@ -560,6 +630,8 @@ int pow2; mpfr_exp_t tmp; + MPFR_LOG_MSG (("case 1 (base = power of 2)\n", 0)); + count_leading_zeros (pow2, (mp_limb_t) pstr->base); pow2 = GMP_NUMB_BITS - pow2 - 1; /* base = 2^pow2 */ MPFR_ASSERTD (0 < pow2 && pow2 <= 5); @@ -596,10 +668,12 @@ mp_limb_t *z; mpfr_exp_t exp_z; + MPFR_LOG_MSG (("case 2 (exp_base > pstr_size)\n", 0)); + result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1); /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */ - z = y - ysize; + z = y0; /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */ err = mpfr_mpn_exp (z, &exp_z, pstr->base, pstr->exp_base - pstr_size, ysize); @@ -646,8 +720,7 @@ /* if the low ysize limbs of {result, 2*ysize} are all zero, then the result is still "exact" (if it was before) */ - exact = exact && (mpn_scan1 (result, 0) - >= (unsigned long) ysize_bits); + exact = exact && (mpn_scan1 (result, 0) >= ysize_bits); result += ysize; } /* case exp_base < pstr_size */ @@ -656,11 +729,12 @@ mp_limb_t *z; mpfr_exp_t exp_z; + MPFR_LOG_MSG (("case 3 (exp_base < pstr_size)\n", 0)); + result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1); - /* set y to y * K^ysize */ - y = y - ysize; /* we have allocated ysize limbs at y - ysize */ - MPN_ZERO (y, ysize); + /* y0 = y * K^ysize */ + MPN_ZERO (y0, ysize); /* pstr_size - pstr->exp_base can overflow */ MPFR_SADD_OVERFLOW (exp_z, (mpfr_exp_t) pstr_size, -pstr->exp_base, @@ -668,33 +742,35 @@ MPFR_EXP_MIN, MPFR_EXP_MAX, goto underflow, goto overflow); - /* (z, exp_z) = base^(exp_base-pstr_size) */ + /* (z, exp_z) = base^(pstr_size - exp_base) */ z = result + 2*ysize + 1; err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize); - /* Since we want y/z rounded toward zero, we must get an upper - bound of z. If err >= 0, the error on z is bounded by 2^err. */ - if (err >= 0) - { - mp_limb_t cy; - unsigned long h = err / GMP_NUMB_BITS; - unsigned long l = err - h * GMP_NUMB_BITS; - if (h >= ysize) /* not enough precision in z */ - goto next_loop; - cy = mpn_add_1 (z, z, ysize - h, MPFR_LIMB_ONE << l); - if (cy != 0) /* the code below requires z on ysize limbs */ - goto next_loop; - } - exact = exact && (err == -1); + /* Now {z, ysize} * 2^(exp_z_out - ysize_bits) is an approximation + to base^exp_z_in (denoted b^e below), rounded toward zero, with: + * if err = -1, the result is exact; + * if err = -2, an overflow occurred in the computation of exp_z; + * otherwise the error is bounded by 2^err ulps. + Thus the exact value of b^e is between z and z + 2^err, where + z is {z, ysize} properly scaled by a power of 2. Then the error + will be: + y/b^e - trunc(y/z) = eps1 + eps2 + with + eps1 = y/b^e - y/z <= 0 + eps2 = y/z - trunc(y/z) >= 0 + thus the errors will (partly) compensate, giving a bound + max(|eps1|,|eps2|). + In addition, there is a 3rd error eps3 since y might be the + conversion of only a part of the character string, and/or y + might be truncated by the mpn_rshift call above: + eps3 = exact_y/b^e - y/b^e >= 0. + */ if (err == -2) goto underflow; /* FIXME: Sure? */ - if (err == -1) - err = 0; - - /* compute y / z */ - /* result will be put into result + n, and remainder into result */ - mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y, - 2 * ysize, z, ysize); + else if (err == -1) + err = 0; /* see the note below */ + else + exact = 0; /* exp -= exp_z + ysize_bits with overflow checking and check that we can add/subtract 2 to exp without overflow */ @@ -706,7 +782,59 @@ mpfr_exp_t, mpfr_uexp_t, MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, goto overflow, goto underflow); - err += 2; + + /* Compute the integer division y/z rounded toward zero. + The quotient will be put at result + ysize (size: ysize + 1), + and the remainder at result (size: ysize). + Both the dividend {y, 2*ysize} and the divisor {z, ysize} are + normalized, i.e., the most significant bit of their most + significant limb is 1. */ + MPFR_ASSERTD (MPFR_LIMB_MSB (y0[2 * ysize - 1]) != 0); + MPFR_ASSERTD (MPFR_LIMB_MSB (z[ysize - 1]) != 0); + mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y0, + 2 * ysize, z, ysize); + + /* The truncation error of the mpn_tdiv_qr call (eps2 above) is at + most 1 ulp. Idem for the error eps3, which has the same sign, + thus eps2 + eps3 <= 2 ulps. + FIXME: For eps3, this is not obvious and should be explained. + For the error eps1 coming from the approximation to b^e, + we have (still up to a power-of-2 normalization): + y/z - y/b^e = y * (b^e-z) / (z * b^e) <= y * 2^err / (z * b^e). + We have to convert that error in terms of ulp(trunc(y/z)). + We first have ulp(trunc(y/z)) = ulp(y/z). + + FIXME: There must be some discussion about the exponents, + because up to a power of 2, 1/2 <= |y/z| < 1 and + 1 <= |y/z| < 2 are equivalent and give no information. + Moreover 1/2 <= b^e < 1 has not been explained and may + hide mistakes since one may have 1/2 <= z < 1 < b^e. + + Since both y and z are normalized, the quotient + {result+ysize, ysize+1} has exactly ysize limbs, plus maybe one + bit (this corresponds to the MPFR_ASSERTD below): + * if the quotient has exactly ysize limbs, then 1/2 <= |y/z| < 1 + (up to a power of 2) and since 1/2 <= b^e < 1, the error is at + most 2^(err+1) ulps; + * if the quotient has one extra bit, then 1 <= |y/z| < 2 + (up to a power of 2) and since 1/2 <= b^e < 1, the error is at + most 2^(err+2) ulps; but since we will shift the result right + below by one bit, the final error will be at most 2^(err+1) ulps + too. + + Thus the error is: + * at most 2^(err+1) ulps for eps1 + * at most 2 ulps for eps2 + eps3, which is of opposite sign + and we can bound the error by 2^(err+1) ulps in all cases. + + Note: If eps1 was 0, the error would be bounded by 2 ulps, + thus replacing err = -1 by err = 0 above was the right thing + to do, since 2^(0+1) = 2. + */ + MPFR_ASSERTD (result[2 * ysize] <= 1); + + err += 1; /* see above for the explanation of the +1 term */ + /* if the remainder of the division is zero, then the result is still "exact" if it was before */ exact = exact && (mpn_popcount (result, ysize) == 0); @@ -726,17 +854,15 @@ /* case exp_base = pstr_size: no multiplication or division needed */ else { + MPFR_LOG_MSG (("case 4 (exp_base = pstr_size)\n", 0)); + /* base^(exp-pr) = 1 nothing to compute */ result = y; err = 0; } - /* If result is exact, we still have to consider the neglected part - of the input string. For a directed rounding, in that case we could - still correctly round, since the neglected part is less than - one ulp, but that would make the code more complex, and give a - speedup for rare cases only. */ - exact = exact && (pstr_size == pstr->prec); + MPFR_LOG_MSG (("exact = %d, err = %d, precx = %Pu\n", + exact, err, precx)); /* at this point, result is an approximation rounded toward zero of the pstr_size most significant digits of pstr->mant, with @@ -744,24 +870,22 @@ /* test if rounding is possible, and if so exit the loop. Note: we also need to be able to determine the correct ternary value, - thus we use the MPFR_PREC(x) + (rnd == MPFR_RNDN) trick. + thus we use the precx + (rnd == MPFR_RNDN) trick. For example if result = xxx...xxx111...111 and rnd = RNDN, then we know the correct rounding is xxx...xx(x+1), but we cannot know the correct ternary value. */ if (exact || mpfr_round_p (result, ysize, ysize_bits - err - 1, - MPFR_PREC(x) + (rnd == MPFR_RNDN))) + precx + (rnd == MPFR_RNDN))) break; - next_loop: /* update the prec for next loop */ MPFR_ZIV_NEXT (loop, prec); } /* loop */ MPFR_ZIV_FREE (loop); /* round y */ - if (mpfr_round_raw (MPFR_MANT (x), result, - ysize_bits, - pstr->negative, MPFR_PREC(x), rnd, &res )) + if (mpfr_round_raw (MPFR_MANT (x), result, ysize_bits, + pstr->negative, precx, rnd, &res)) { /* overflow when rounding y */ MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT; @@ -769,12 +893,11 @@ exp ++; } - if (res == 0) /* fix ternary value */ - { - exact = exact && (pstr_size == pstr->prec); - if (!exact) - res = (pstr->negative) ? 1 : -1; - } + /* Note: if exact <> 0, then the approximation {result, ysize} is exact, + thus no double-rounding can occur: + (a) either the ternary value res is non-zero, and it is the correct + ternary value that we should return + (b) or the ternary value res is zero, and we should return 0. */ /* Set sign of x before exp since check_range needs a valid sign */ (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x); diff -Naurd mpfr-4.0.1-a/src/version.c mpfr-4.0.1-b/src/version.c --- mpfr-4.0.1-a/src/version.c 2018-12-31 11:18:05.621387341 +0000 +++ mpfr-4.0.1-b/src/version.c 2019-01-31 12:05:56.091705514 +0000 @@ -25,5 +25,5 @@ const char * mpfr_get_version (void) { - return "4.0.1-p14"; + return "4.0.1-p15"; } diff -Naurd mpfr-4.0.1-a/tests/tstrtofr.c mpfr-4.0.1-b/tests/tstrtofr.c --- mpfr-4.0.1-a/tests/tstrtofr.c 2018-01-09 12:30:58.000000000 +0000 +++ mpfr-4.0.1-b/tests/tstrtofr.c 2019-01-31 12:05:56.051705613 +0000 @@ -22,6 +22,12 @@ #include "mpfr-test.h" +/* The implicit \0 is useless, but we do not write num_to_text[62] otherwise + g++ complains. */ +static const char num_to_text36[] = "0123456789abcdefghijklmnopqrstuvwxyz"; +static const char num_to_text62[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" + "abcdefghijklmnopqrstuvwxyz"; + static void check_special (void) { @@ -898,6 +904,18 @@ { printf ("Check overflow failed (2) with:\n s='%s'\n x=", s); mpfr_dump (x); +#if defined(__GNUC__) + printf ("This failure is triggered by GCC bug 86554:\n" + " https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86554\n" + " https://gcc.gnu.org/bugzilla/show_bug.cgi?id=87276 " + "(about this test)\nWorkaround: disable code hoisting " + "with -fno-code-hoisting in CFLAGS.\n"); + /* Note: In Debian, this error is obtained with gcc-snapshot from + 20180908-1 to 20181127-1. With gcc-snapshot from 20181209-1 to + 20190102-1 (at least), the MPFR build no longer seems affected + in general, but using --with-gmp-build=... together with + --enable-assert still triggers this failure. */ +#endif exit (1); } mpfr_strtofr (x, "123456789E170141183460469231731687303715884105728", @@ -938,6 +956,20 @@ mpfr_dump (x); exit (1); } + mpfr_strtofr (x, "1@2305843009213693951", &s, 16, MPFR_RNDN); + if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x)) + { + printf ("Check overflow failed (8) with:\n s=%s\n x=", s); + mpfr_dump (x); + exit (1); + } + mpfr_strtofr (x, "1@2305843009213693951", &s, 17, MPFR_RNDN); + if (s[0] != 0 || !MPFR_IS_INF (x) || !MPFR_IS_POS (x)) + { + printf ("Check overflow failed (9) with:\n s=%s\n x=", s); + mpfr_dump (x); + exit (1); + } /* Check underflow */ mpfr_strtofr (x, "123456789E-2147483646", &s, 0, MPFR_RNDN); @@ -969,6 +1001,13 @@ mpfr_dump (x); exit (1); } + mpfr_strtofr (x, "1@-2305843009213693952", &s, 16, MPFR_RNDN); + if (s[0] != 0 || !MPFR_IS_ZERO (x) || !MPFR_IS_POS (x) ) + { + printf ("Check underflow failed (8) with:\n s='%s'\n x=", s); + mpfr_dump (x); + exit (1); + } mpfr_clear (x); } @@ -1238,12 +1277,225 @@ int inex; emin = mpfr_get_emin (); - mpfr_set_emin (-1073); - mpfr_set_emin (emin); mpfr_init2 (z, 53); + mpfr_set_emin (-1073); + /* with emin = -1073, the smallest positive number is 0.5*2^emin = 2^(-1074), + thus str should be rounded to 2^(-1074) with inex > 0 */ + inex = mpfr_strtofr (z, str, NULL, 10, MPFR_RNDN); + MPFR_ASSERTN(inex > 0 && mpfr_cmp_ui_2exp (z, 1, -1074) == 0); + mpfr_set_emin (-1074); + /* with emin = -1074, str should be rounded to 2^(-1075) with inex < 0 */ inex = mpfr_strtofr (z, str, NULL, 10, MPFR_RNDN); MPFR_ASSERTN(inex < 0 && mpfr_cmp_ui_2exp (z, 1, -1075) == 0); mpfr_clear (z); + mpfr_set_emin (emin); +} + +/* r13299 fails with 8-bit limbs (micro-gmp/8). */ +static void +bug20181127 (void) +{ + char s[] = "9.Z6nrLVSMG1bDcCF2ONJdX@-183295525"; /* base 58 */ + mpfr_t x, y; + + mpfr_inits2 (6, x, y, (mpfr_ptr) 0); + mpfr_set_ui_2exp (x, 5, -1073741701, MPFR_RNDN); + mpfr_strtofr (y, s, NULL, 58, MPFR_RNDZ); + if (! mpfr_equal_p (x, y)) + { + printf ("Error in bug20181127 on %s (base 58)\n", s); + printf ("Expected x = "); + mpfr_dump (x); + printf ("Got y = "); + mpfr_dump (y); + printf ("*Note* In base 58, x ~= "); + mpfr_out_str (stdout, 58, 8, x, MPFR_RNDN); + printf ("\n"); + exit (1); + } + mpfr_clears (x, y, (mpfr_ptr) 0); +} + +static void +coverage (void) +{ +#if _MPFR_EXP_FORMAT >= 3 && _MPFR_PREC_FORMAT == 3 && MPFR_PREC_BITS == 64 + char str3[] = "1@-2112009130072406892"; + char str31[] = "1@-511170973314085831"; + mpfr_t x; + int inex; + mpfr_exp_t emin; + + /* exercise assertion cy == 0 around line 698 of strtofr.c */ + emin = mpfr_get_emin (); + mpfr_set_emin (mpfr_get_emin_min ()); + /* emin = -4611686018427387903 on a 64-bit machine */ + mpfr_init2 (x, 1); + inex = mpfr_strtofr (x, str3, NULL, 3, MPFR_RNDN); + /* 3^-2112009130072406892 is slightly larger than (2^64)^-52303988630398057 + thus we should get inex < 0 */ + MPFR_ASSERTN(inex < 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, -52303988630398057 * 64) == 0); + inex = mpfr_strtofr (x, str31, NULL, 31, MPFR_RNDN); + /* 31^-511170973314085831 is slightly smaller than (2^64)^-39569396093273623 + thus we should get inex > 0 */ + MPFR_ASSERTN(inex > 0); + MPFR_ASSERTN(mpfr_cmp_ui_2exp (x, 1, -39569396093273623 * 64) == 0); + mpfr_clear (x); + mpfr_set_emin (emin); +#endif +} + +#define BSIZE 512 + +static void +random_tests (void) +{ + char s0[BSIZE], s1[BSIZE], s2[BSIZE+64]; + mpfr_t x0, x1, x2; + int prec, i; + + for (prec = MPFR_PREC_MIN; prec < 300; prec++) + { + mpfr_inits2 (prec, x0, x1, x2, (mpfr_ptr) 0); + + for (i = 0; i < 20; i++) + { + const char *num_to_text; + mpfr_exp_t e0, e1; + int base, j, k, neg; + int noteq = 0; + char d; + + /* We want the same exponent for x0 and its successor x1. + This is not possible for precision 1 in base 2. */ + do + base = 2 + (randlimb () % 61); + while (prec == 1 && base == 2); + + num_to_text = base <= 36 ? num_to_text36 : num_to_text62; + + do + { + /* Let's consider only positive numbers. We should test + negative numbers, but they should be built later, just + for the test itself. */ + tests_default_random (x0, 0, + mpfr_get_emin (), mpfr_get_emax (), 1); + mpfr_set (x1, x0, MPFR_RNDN); + mpfr_nextabove (x1); + mpfr_get_str (s0, &e0, base, BSIZE - 1, x0, MPFR_RNDU); + mpfr_get_str (s1, &e1, base, BSIZE - 1, x1, MPFR_RNDD); + } + while (! (mpfr_regular_p (x0) && mpfr_regular_p (x1) && e0 == e1)); + + /* 0 < x0 <= (s0,e) <= (s1,e) <= x1 with e = e0 = e1. + Let's build a string s2 randomly formed by: + - the common prefix of s0 and s1; + - some of the following digits of s0 (possibly none); + - the next digit of s0 + 1; + - some of the following digits of s1 (possibly none). + Then 0 < x0 <= (s0,e) < (s2,e) <= (s1,e) <= x1, and with + a very high probability that (s2,e) < (s1,e); noteq is + set to true in this case. + For instance, if: + s0 = 123456789 + s1 = 124012345 + one can have, e.g.: + s2 = 12345734 + s2 = 123556789 + s2 = 124 + s2 = 124012 + s2 is not taken completely randomly between s0 and s1, but it + will be built rather easily, and with the randomness of x0, + we should cover all cases, with s2 very close to s0, s2 very + close to s1, or not too close to either. */ + + neg = randlimb () & 1; + s2[0] = neg ? '-' : '+'; + s2[1] = '.'; + + for (j = 0; + MPFR_ASSERTN (s0[j] != 0 && s1[j] != 0), s0[j] == s1[j]; + j++) + s2[j+2] = s0[j]; + + /* k is the position of the first differing digit. */ + k = j; + + while (j < BSIZE - 2 && randlimb () % 8 != 0) + { + MPFR_ASSERTN (s0[j] != 0); + s2[j+2] = s0[j]; + j++; + } + + MPFR_ASSERTN (s0[j] != 0); + /* We will increment the next digit. Thus while s0[j] is the + maximum digit, go back until this is no longer the case + (the first digit after the common prefix cannot be the + maximum digit, so that we will stop early enough). */ + while ((d = s0[j]) == num_to_text[base - 1]) + j--; + noteq = j != k; + s2[j+2] = d = *(strchr (num_to_text, d) + 1); + if (d != s1[j]) + noteq = 1; + j++; + + while (j < BSIZE - 1 && randlimb () % 8 != 0) + { + MPFR_ASSERTN (s1[j] != 0); + s2[j+2] = s1[j]; + j++; + } + + sprintf (s2 + (j+2), "@%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) e0); + + while (noteq == 0 && j < BSIZE - 1) + { + if (s1[j] != '0') + noteq = 1; + j++; + } + + if (neg) + { + mpfr_neg (x0, x0, MPFR_RNDN); + mpfr_neg (x1, x1, MPFR_RNDN); + } + + if (noteq) + { + mpfr_strtofr (x2, s2, NULL, base, MPFR_RNDZ); + if (! mpfr_equal_p (x2, x0)) + { + printf ("Error in random_tests for prec=%d i=%d base=%d\n", + prec, i, base); + printf ("s0 = %s\ns1 = %s\ns2 = %s\n", s0, s1, s2); + printf ("x0 = "); + mpfr_dump (x0); + printf ("x2 = "); + mpfr_dump (x2); + exit (1); + } + } + + mpfr_strtofr (x2, s2, NULL, base, MPFR_RNDA); + if (! mpfr_equal_p (x2, x1)) + { + printf ("Error in random_tests for prec=%d i=%d base=%d\n", + prec, i, base); + printf ("s0 = %s\ns1 = %s\ns2 = %s\n", s0, s1, s2); + printf ("x1 = "); + mpfr_dump (x1); + printf ("x2 = "); + mpfr_dump (x2); + exit (1); + } + } + mpfr_clears (x0, x1, x2, (mpfr_ptr) 0); + } } int @@ -1251,6 +1503,7 @@ { tests_start_mpfr (); + coverage (); check_special(); check_reftable (); check_parse (); @@ -1262,6 +1515,8 @@ bug20120829 (); bug20161217 (); bug20170308 (); + bug20181127 (); + random_tests (); tests_end_mpfr (); return 0;