improve mpc_asin for tiny imaginary part

2 files changed, 37 insertions, 4 deletions

diff --git a/src/asin.c b/src/asin.c
index 96749db..4cf2353 100644
--- a/src/asin.c
+++ b/src/asin.c
@@ -87,7 +87,7 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   mpfr_rnd_t rnd_re, rnd_im;
   mpc_t z1;
   int inex, inex_re, inex_im, loop = 0;
-  mpfr_exp_t saved_emin, saved_emax, err;
+  mpfr_exp_t saved_emin, saved_emax, err, olderr;
 
   /* special values */
   if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
@@ -206,13 +206,14 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   rnd_im = MPC_RND_IM(rnd);
   p = p_re >= p_im ? p_re : p_im;
   mpc_init2 (z1, p);
-  err = 0; /* number of lost bits */
+  olderr = err = 0; /* number of lost bits */
   while (1)
   {
     mpfr_exp_t ex, ey;
 
     loop ++;
-    p += err; /* add number of lost bits in previous loop */
+    p += err - olderr; /* add extra number of lost bits in previous loop */
+    olderr = err;
     p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2;
     mpfr_set_prec (mpc_realref(z1), p);
     mpfr_set_prec (mpc_imagref(z1), p);
diff --git a/tests/tacos.c b/tests/tacos.c
index 2f764b1..8463bbb 100644
--- a/tests/tacos.c
+++ b/tests/tacos.c
@@ -44,7 +44,6 @@ test_reduced (void)
   mpfr_set_flt (mpc_realref (x), -0x0.01f28c10p0f, MPFR_RNDN);
   mpfr_set_flt (mpc_imagref (x), -0x6.79cdd0p-68f, MPFR_RNDN);
   mpc_acos (z, x, MPC_RNDNN);
-  /* square root is 0.00038652126908433 - 1.37634353022868e12*I */
   mpfr_set_flt (mpc_realref (zr), 0x1.941242p0f, MPFR_RNDN);
   mpfr_set_flt (mpc_imagref (zr), 0x6.79da18p-68f, MPFR_RNDN);
   if (mpc_cmp (z, zr))
@@ -61,12 +60,45 @@ test_reduced (void)
   mpfr_set_emax (emax);
 }
 
+/* another test with reduced exponent range */
+static void
+test_reduced2 (void)
+{
+  mpfr_exp_t emin = mpfr_get_emin ();
+  mpfr_exp_t emax = mpfr_get_emax ();
+  mpc_t x, z, zr;
+
+  mpfr_set_emin (-1073);
+  mpfr_set_emax (1024);
+  mpc_init2 (x, 53);
+  mpc_init2 (z, 53);
+  mpc_init2 (zr, 53);
+  mpfr_set_d (mpc_realref (x), -2.2694475687286223e-15, MPFR_RNDN);
+  mpfr_set_d (mpc_imagref (x), 2.7236935900137536e-309, MPFR_RNDN);
+  mpc_acos (z, x, MPC_RNDNN);
+  mpfr_set_d (mpc_realref (zr), 1.5707963267948988, MPFR_RNDN);
+  mpfr_set_d (mpc_imagref (zr), -2.7236935900137536e-309, MPFR_RNDN);
+  if (mpc_cmp (z, zr))
+    {
+      printf ("Incorrect acos 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 ();
 
   test_reduced ();
+  test_reduced2 ();
 
   data_check_template ("acos.dsc", "acos.dat");