increase code coverage to 99.9%

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1199 211d60ee-9f03-0410-a15a-8952a2c7a4e4

8 files changed, 84 insertions, 60 deletions

diff --git a/src/div.c b/src/div.c
index c96832d..83584b8 100644
--- a/src/div.c
+++ b/src/div.c
@@ -312,45 +312,53 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
             hopefully, the side-effects of mpc_mul do indeed raise the
             mpfr exceptions */
       if (overflow_prod) {
+         int isinf = 0;
          tmpsgn = mpfr_sgn (mpc_realref(res));
          if (tmpsgn > 0)
-            mpfr_nextabove (mpc_realref(res));
+           {
+             mpfr_nextabove (mpc_realref(res));
+             isinf = mpfr_inf_p (mpc_realref(res));
+             mpfr_nextbelow (mpc_realref(res));
+           }
          else if (tmpsgn < 0)
-            mpfr_nextbelow (mpc_realref(res));
-         if (mpfr_inf_p (mpc_realref(res)))
+           {
+             mpfr_nextbelow (mpc_realref(res));
+             isinf = mpfr_inf_p (mpc_realref(res));
+             mpfr_nextabove (mpc_realref(res));
+           }
+         if (isinf)
            {
              mpfr_set_inf (mpc_realref(res), tmpsgn);
              overflow_re = 1;
            }
-         else
-            if (tmpsgn > 0)
-               mpfr_nextbelow (mpc_realref(res));
-            else if (tmpsgn < 0)
-               mpfr_nextabove (mpc_realref(res));
          tmpsgn = mpfr_sgn (mpc_imagref(res));
+         isinf = 0;
          if (tmpsgn > 0)
-            mpfr_nextabove (mpc_imagref(res));
+           {
+             mpfr_nextabove (mpc_imagref(res));
+             isinf = mpfr_inf_p (mpc_imagref(res));
+             mpfr_nextbelow (mpc_imagref(res));
+           }
          else if (tmpsgn < 0)
-            mpfr_nextbelow (mpc_imagref(res));
-         if (mpfr_inf_p (mpc_imagref(res)))
+           {
+             mpfr_nextbelow (mpc_imagref(res));
+             isinf = mpfr_inf_p (mpc_imagref(res));
+             mpfr_nextabove (mpc_imagref(res));
+           }
+         if (isinf)
            {
              mpfr_set_inf (mpc_imagref(res), tmpsgn);
              overflow_im = 1;
            }
-         else
-            if (tmpsgn > 0)
-               mpfr_nextbelow (mpc_imagref(res));
-            else if (tmpsgn < 0)
-               mpfr_nextabove (mpc_imagref(res));
          mpc_set (a, res, rnd);
          goto end;
       }
 
       /* divide the product by the norm */
       if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) {
-         /* The division has good chances to be exact in at least one part.   */
-         /* Since this can cause problems when not rounding to the nearest,   */
-         /* we use the division code of mpfr, which handles the situation.    */
+         /* The division has good chances to be exact in at least one part.  */
+         /* Since this can cause problems when not rounding to the nearest,  */
+         /* we use the division code of mpfr, which handles the situation.   */
          mpfr_clear_underflow ();
          mpfr_clear_overflow ();
          inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ);
@@ -416,16 +424,16 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
       inexact_re = mpfr_sgn (mpc_realref (res));
    }
    else if (underflow_re || (overflow_norm && !overflow_prod)) {
-      mpfr_set_zero (mpc_realref (a), mpfr_sgn (mpc_realref (res)));
-      inexact_re = -mpfr_sgn (mpc_realref (res));
+      inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1;
+      mpfr_set_zero (mpc_realref (a), -inexact_re);
    }
    if (overflow_im || (underflow_norm && !underflow_prod)) {
       mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res)));
       inexact_im = mpfr_sgn (mpc_imagref (res));
    }
    else if (underflow_im || (overflow_norm && !overflow_prod)) {
-      mpfr_set_zero (mpc_imagref (a), mpfr_sgn (mpc_imagref (res)));
-      inexact_im = -mpfr_sgn (mpc_imagref (res));
+      inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1;
+      mpfr_set_zero (mpc_imagref (a), -inexact_im);
    }
 
    mpc_clear (res);
diff --git a/src/log10.c b/src/log10.c
index e5972cc..10a2566 100644
--- a/src/log10.c
+++ b/src/log10.c
@@ -254,7 +254,7 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
           mpfr_integer_p (mpc_imagref (op)))
         {
           mpz_t x, y;
-          unsigned long s;
+          unsigned long s, v;
 
           mpz_init (x);
           mpz_init (y);
@@ -263,23 +263,23 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
           mpz_mul (x, x, x);
           mpz_mul (y, y, y);
           mpz_add (x, x, y); /* x^2+y^2 */
-          s = mpz_sizeinbase (x, 10); /* always exact or 1 too big */
-          /* since s is the number of digits of x in base 10 (or one more),
-             we subtract 1 since we want to check whether x = 10^s */
-          s --;
-          mpz_ui_pow_ui (y, 10, s);
-          if (mpz_cmp (y, x) > 0) /* might be 1 too big */
+          v = mpz_scan1 (x, 0);
+          /* if x = 10^s then necessarily s = v */
+          s = mpz_sizeinbase (x, 10);
+          /* since s is either the number of digits of x or one more,
+             then x = 10^(s-1) or 10^(s-2) */
+          if (s == v + 1 || s == v + 2)
             {
-              mpz_divexact_ui (y, y, 10);
-              s --;
-            }
-          if (mpz_cmp (y, x) == 0) /* Re(log10(x+i*y)) is exactly s/2 */
-            {
-              /* we reset the precision of Re(ww) so that s can be
-                 represented exactly */
-              mpfr_set_prec (mpc_realref (ww), sizeof(unsigned long)*CHAR_BIT);
-              mpfr_set_ui_2exp (mpc_realref (ww), s, -1, GMP_RNDN); /* exact */
-              ok = 1;
+              mpz_div_2exp (x, x, v);
+              mpz_ui_pow_ui (y, 5, v);
+              if (mpz_cmp (y, x) == 0) /* Re(log10(x+i*y)) is exactly v/2 */
+                {
+                  /* we reset the precision of Re(ww) so that v can be
+                     represented exactly */
+                  mpfr_set_prec (mpc_realref (ww), sizeof(unsigned long)*CHAR_BIT);
+                  mpfr_set_ui_2exp (mpc_realref (ww), v, -1, GMP_RNDN); /* exact */
+                  ok = 1;
+                }
             }
           mpz_clear (x);
           mpz_clear (y);
diff --git a/src/pow.c b/src/pow.c
index 92abca7..892f467 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -202,7 +202,7 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
   t = mpz_scan1 (my, 0);
   ey += (mpfr_exp_t) t;
   mpz_tdiv_q_2exp (my, my, t);
-  /* y = my*2^ey */
+  /* y = my*2^ey with my odd */
 
   if (x_imag)
     {
@@ -229,7 +229,6 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
       mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
       if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
         goto end;
-      ed = ec;
     }
   else if (ed < ec)
     {
@@ -484,7 +483,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
   mpc_t t, u;
   mpfr_prec_t p, pr, pi, maxprec;
   int saved_underflow, saved_overflow;
-
+  
   /* save the underflow or overflow flags from MPFR */
   saved_underflow = mpfr_underflow_p ();
   saved_overflow = mpfr_overflow_p ();
@@ -614,8 +613,8 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
          (a) x is negative and y is half-an-integer
          (b) x = -1 and Re(y) is half-an-integer
       */
-      if (mpfr_cmp_ui (mpc_realref(x), 0) < 0 && is_odd (mpc_realref(y), 1) &&
-          (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
+      if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1)
+         && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
         z_imag = 1;
     }
   else /* x non real */
diff --git a/src/sqrt.c b/src/sqrt.c
index 9250c27..dd2ff60 100644
--- a/src/sqrt.c
+++ b/src/sqrt.c
@@ -1,6 +1,6 @@
 /* mpc_sqrt -- Take the square root of a complex number.
 
-Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA
+Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 INRIA
 
 This file is part of GNU MPC.
 
@@ -329,11 +329,11 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
                inex_im = inex_t;
             else {
                inex_im = -inex_t;
-               if (   (im_cmp > 0 && r == GMP_RNDD)
-                   || (im_cmp < 0 && r == GMP_RNDU))
-                  MPFR_ADD_ONE_ULP (mpc_imagref (a));
-               else
-                  MPFR_SUB_ONE_ULP (mpc_imagref (a));
+               /* im_cmp > 0 implies that Im(b) > 0, thus im_sgn = 0
+                  and r = GMP_RNDU.
+                  im_cmp < 0 implies that Im(b) < 0, thus im_sgn = -1
+                  and r = GMP_RNDD. */
+               MPFR_SUB_ONE_ULP (mpc_imagref (a));
             }
          }
          else if (im_cmp > 0) {
@@ -341,21 +341,17 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
                inex_re = inex_t;
             else {
                inex_re = -inex_t;
-               if (r == GMP_RNDD)
-                  MPFR_ADD_ONE_ULP (mpc_realref (a));
-               else
-                  MPFR_SUB_ONE_ULP (mpc_realref (a));
+               /* im_cmp > 0 implies r = GMP_RNDU (see above) */
+               MPFR_SUB_ONE_ULP (mpc_realref (a));
             }
          }
-         else {
+         else { /* im_cmp < 0 */
             if (rnd_t == r)
                inex_re = -inex_t;
             else {
                inex_re = inex_t;
-               if (r == GMP_RNDD)
-                  MPFR_SUB_ONE_ULP (mpc_realref (a));
-               else
-                  MPFR_ADD_ONE_ULP (mpc_realref (a));
+               /* im_cmp < 0 implies r = GMP_RNDD (see above) */
+               MPFR_SUB_ONE_ULP (mpc_realref (a));
             }
          }
       }
diff --git a/tests/div.dat b/tests/div.dat
index dbdc2b1..6f5e600 100644
--- a/tests/div.dat
+++ b/tests/div.dat
@@ -2478,3 +2478,9 @@
 0 + 2 0 2 inf 10 0x3ffp1073741813 10 0x3ffp1073741813 10 0x2abp-10 10 -0x2abp-10 N N
 #  negative overflow
 0 - 2 0 2 -inf 10 -0x3ffp1073741813 10 -0x3ffp1073741813 10 0x2abp-10 10 -0x2abp-10 N N
+
+# examples to exercise underflow
+# (1.5+i)*2^emin/(1-i) gives (0.25 + 1.25*i)*2^emin
+- - 2 0 2 0x1p-1073741823 2 0x3p-1073741824 2 0x1p-1073741823 2 1 2 -1 Z Z
+# (1.5+i)*2^emin/(1+i) gives (1.25 - 0.25*i)*2^emin
+- + 2 0x1p-1073741823 2 -0 2 0x3p-1073741824 2 0x1p-1073741823 2 1 2 1 Z Z
diff --git a/tests/pow.dat b/tests/pow.dat
index 15eab18..1105f59 100644
--- a/tests/pow.dat
+++ b/tests/pow.dat
@@ -454,3 +454,16 @@
 - + 2 +0 2 -0 2 4 2 3 28 -744261116 2 +0 N N
 - - 2 +0 2 +0 2 4 2 -3 28 -744261116 2 +0 N N
 
+# exact powers with non-integer exponent
+0 0 2 1 2 1 2 0 2 2 2 0.5 2 0 N N
+0 0 2 -2 2 2 2 0 2 2 2 1.5 2 0 N N
+0 0 2 1 2 64 12 -4095 2 128 2 0.5 2 0 N N
+0 0 3 5 2 3 2 16 4 30 2 0.5 2 0 N N
+0 0 7 97 7 99 6 -392 14 19206 2 0.5 2 0 N N
+0 0 6 63 6 61 5 248 18 7686 2 0.5 2 0 N N
+0 0 6 63 6 61 24 -59013092 17 3812256 2 0.25 2 0 N N
+
+0 + 2 0 2 0x3p-6 2 -1 2 0 2 0.5 2 1 N N
++ + 2 6 2 1 41 -0x2ce019e6f1e 36 0x1878418ba20 2 0.0625 2 0 N N
++ + 4 11 2 1 111 -0x73558286726957f922819cbeffff 109 0x1c484a8b32dbf409e966a8c00000 2 0x1p-5 2 0 N N
++ + 5 21 2 1 282 -0x24ea91ddba938e750d999f1075444e15d6ca0fff6a19c8cbefe6260261fd57effffffff 278 0x390aa828a3d933391ab999b0b0aa71aafbfc7b127fe30c84d107634940ba8000000000 2 0x1p-6 2 0 N N
diff --git a/tests/sqrt.dat b/tests/sqrt.dat
index 8a9e7ca..076f2af 100644
--- a/tests/sqrt.dat
+++ b/tests/sqrt.dat
@@ -123,6 +123,7 @@
 - + 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281177 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 N N
 + + 375 0x1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281177 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 U U
 - - 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef5fffffffffffffffffffffffffffffffffffffffffffep-202281177 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 D D
+- + 375 1 375 -0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef5fffffffffffffffffffffffffffffffffffffffffffep-202281177 375 1 375 -0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 D U
 - - 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef5fffffffffffffffffffffffffffffffffffffffffffep-202281177 375 1 375 0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 Z Z
 - - 375 1 375 -0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281177 375 1 375 -0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 N N
 + + 375 0x1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004 375  -0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef5fffffffffffffffffffffffffffffffffffffffffffep-202281177 375 1 375 -0xf.8a8aae3080b3dd665e316d262fd54c1ca22a83dc9acb92ef6p-202281176 U U
diff --git a/tests/tan.dat b/tests/tan.dat
index d9ce811..ce4c097 100644
--- a/tests/tan.dat
+++ b/tests/tan.dat
@@ -125,6 +125,7 @@
 
 # huge values
 + - 53 -0 53 -1  53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N
+- - 53 +0 53 -1  53 -0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N
 + + 53 -0 53 +1  53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N
 
 # some values taken from ttan.c