testing/mpfr
SHA256
3
0
Fork 0
forked from pool/mpfr
Code Pull Requests Activity
mpfr/mpfr-4.0.1-cummulative-patch.patch

1975 lines
69 KiB
Diff
Raw Normal View History

- Add mpfr-4.0.1-cummulative-patch.patch. Fixes * A subtraction of two numbers of the same sign or addition of two numbers of different signs can be rounded incorrectly (and the ternary value can be incorrect) when one of the two inputs is reused as the output (destination) and all these MPFR numbers have exactly GMP_NUMB_BITS bits of precision (typically, 32 bits on 32-bit machines, 64 bits on 64-bit machines). * The mpfr_fma and mpfr_fms functions can behave incorrectly in case of internal overflow or underflow. * The result of the mpfr_sqr function can be rounded incorrectly in a rare case near underflow when the destination has exactly GMP_NUMB_BITS bits of precision (typically, 32 bits on 32-bit machines, 64 bits on 64-bit machines) and the input has at most GMP_NUMB_BITS bits of precision. * The behavior and documentation of the mpfr_get_str function are inconsistent concerning the minimum precision (this is related to the change of the minimum precision from 2 to 1 in MPFR 4.0.0). The get_str patch fixes this issue in the following way: the value 1 can now be provided for n (4th argument of mpfr_get_str); if n = 0, then the number of significant digits in the output string can now be 1, as already implied by the documentation (but the code was increasing it to 2). * The mpfr_cmp_q function can behave incorrectly when the rational (mpq_t) number has a null denominator. * The mpfr_inp_str and mpfr_out_str functions might behave incorrectly when the stream is a null pointer: the stream is replaced by stdin and stdout, respectively. This behavior is useless, not documented (thus incorrect in case a null pointer would have a special meaning), and not consistent with other input/output functions. OBS-URL: https://build.opensuse.org/package/show/devel:libraries:c_c++/mpfr?expand=0&rev=51
2018-05-02 09:02:53 +00:00
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.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/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/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";
}
Copy Permalink