File mpfr-4.0.1-cummulative-patch.patch of Package mpfr

Overview Repositories Revisions Requests Users Attributes Meta

File mpfr-4.0.1-cummulative-patch.patch of Package mpfr (Revision f87873b9c4e8a5f1399ec4e76712ed49)

Currently displaying revision f87873b9c4e8a5f1399ec4e76712ed49 , Show latest

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";
 }

Places

File mpfr-4.0.1-cummulative-patch.patch of Package mpfr (Revision f87873b9c4e8a5f1399ec4e76712ed49)

Places