[src/sqr.c] Fixed a bug in mpfr_sqr_1n in a rare case near underflow.

[tests/tsqr.c] Added tests, including non-regression for above bug
(manually patched src/sqr.c since r12398 had other, unrelated changes;
merged changesets r12398-12399 on tests/tsqr.c from the trunk)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/4.0@12629 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 163 insertions, 9 deletions

diff --git a/src/sqr.c b/src/sqr.c
index bfae01853..dd71f2355 100644
--- a/src/sqr.c
+++ b/src/sqr.c
@@ -161,15 +161,18 @@ mpfr_sqr_1n (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
   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 --git a/tests/tsqr.c b/tests/tsqr.c
index 6fb9806b8..3701c146d 100644
--- a/tests/tsqr.c
+++ b/tests/tsqr.c
@@ -188,6 +188,156 @@ check_mpn_sqr (void)
 #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 @@ main (void)
 
   tests_start_mpfr ();
 
+  coverage (1024);
   check_mpn_sqr ();
   check_special ();
   test_underflow ();