fixed sqrt for reduced exponent range

3 files changed, 49 insertions, 1 deletions

diff --git a/src/sqrt.c b/src/sqrt.c
index 048a54c..8f58cbc 100644
--- a/src/sqrt.c
+++ b/src/sqrt.c
@@ -48,6 +48,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
      /* we need to know the sign of Im(b) when it is +/-0 */
   const mpfr_rnd_t r = im_sgn ? MPFR_RNDD : MPFR_RNDU;
      /* rounding mode used when computing t */
+  mpfr_exp_t saved_emin, saved_emax;
 
   /* special values */
   if (!mpc_fin_p (b)) {
@@ -202,6 +203,11 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
       }
    }
 
+  saved_emin = mpfr_get_emin ();
+  saved_emax = mpfr_get_emax ();
+  mpfr_set_emin (mpfr_get_emin_min ());
+  mpfr_set_emax (mpfr_get_emax_max ());
+
   do
     {
       loops ++;
@@ -363,5 +369,11 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
   mpfr_clear (w);
   mpfr_clear (t);
 
+  /* restore the exponent range, and check the range of results */
+  mpfr_set_emin (saved_emin);
+  mpfr_set_emax (saved_emax);
+  inex_re = mpfr_check_range (mpc_realref (a), inex_re, MPC_RND_RE (rnd));
+  inex_im = mpfr_check_range (mpc_imagref (a), inex_im, MPC_RND_IM (rnd));
+
   return MPC_INEX (inex_re, inex_im);
 }
diff --git a/tests/sqrt.dat b/tests/sqrt.dat
index 0e5be87..432b5d9 100644
--- a/tests/sqrt.dat
+++ b/tests/sqrt.dat
@@ -113,8 +113,10 @@
 - - 53 0x5a827999fcef30p-54 53 -0x16a09e667f3bcdp-52    25 +0 67 -4 Z D
 + - 53 0x16a09e667f3bcdp-52 53 -0x16a09e667f3bcdp-52    23 -0 68 -4 U D
 
-# x=(-1.95050937e+25,-2.65087379e-39)
+# z=(-1.95050937e+25,-2.65087379e-39)
 + - 24 0x7.2fb038p-164 24 -0x1.011258p32 24 -0x1.0225d4p64 24 -0xe.6ec760p-132 N N
+# z=(-1.89432151e+24,-1.06397210e+09)
+? ? 24 0x0.001954c04p0 24 -0x1.407478p40 24 -0x1.912360p80 24 -0x3.f6aed0p28 N N
 
 # bugs fixed in r160 2008-07-15
 - + 19 0b11101001001001001100p+39 19 -0b1010110101100111011p-236    19 0b1.101010001010100000p+117 19 -0b1.001110111101100001p-158 N Z
diff --git a/tests/tsqrt.c b/tests/tsqrt.c
index 17d05c3..15f9888 100644
--- a/tests/tsqrt.c
+++ b/tests/tsqrt.c
@@ -28,11 +28,45 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 #include "data_check.tpl"
 #include "tgeneric.tpl"
 
+/* check with reduced exponent range */
+static void
+bug20200207 (void)
+{
+  mpfr_exp_t emin = mpfr_get_emin ();
+  mpfr_exp_t emax = mpfr_get_emax ();
+  mpc_t x, z, zr;
+
+  mpfr_set_emin (-148);
+  mpfr_set_emax (128);
+  mpc_init2 (x, 24);
+  mpc_init2 (z, 24);
+  mpc_init2 (zr, 24);
+  mpfr_set_d (mpc_realref (x), -1.89432151320234e24, MPFR_RNDN);
+  mpfr_set_d (mpc_imagref (x), -1.06397209600000e9, MPFR_RNDN);
+  mpc_sqrt (z, x, MPC_RNDNN);
+  /* square root is 0.00038652126908433 - 1.37634353022868e12*I */
+  mpfr_set_d (mpc_realref (zr), 0.00038652126908433, MPFR_RNDN);
+  mpfr_set_d (mpc_imagref (zr), -1.37634353022868e12, MPFR_RNDN);
+  if (mpc_cmp (z, zr))
+    {
+      printf ("Incorrect square root with reduced exponent range:\n");
+      mpfr_printf ("Expected (%Re,%Re)\n", mpc_realref (zr), mpc_imagref (zr));
+      mpfr_printf ("Got      (%Re,%Re)\n", mpc_realref (z), mpc_imagref (z));
+      exit (1);
+    }
+  mpc_clear (x);
+  mpc_clear (z);
+  mpc_clear (zr);
+  mpfr_set_emin (emin);
+  mpfr_set_emax (emax);
+}
+
 int
 main (void)
 {
   test_start ();
 
+  bug20200207 ();
   data_check_template ("sqrt.dsc", "sqrt.dat");
 
   tgeneric_template ("sqrt.dsc", 2, 1024, 7, 256);