[read_data.c] added random tests to check that MPC does not *clear* MPFR flags

	      (some flags might be *set* internally however)
[pow.c,exp.c,norm.c,div.c] fix issues with MPFR flags which were cleared
			   but not reset to their initial value
[div.dat] updated current result of some commented out failing tests


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

6 files changed, 164 insertions, 13 deletions

diff --git a/src/div.c b/src/div.c
index 0dbf29c..3186333 100644
--- a/src/div.c
+++ b/src/div.c
@@ -1,6 +1,6 @@
 /* mpc_div -- Divide two complex numbers.
 
-Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011 INRIA
+Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA
 
 This file is part of GNU MPC.
 
@@ -20,18 +20,23 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 
 #include "mpc-impl.h"
 
+/* this routine deals with the case where w is zero */
 static int
 mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
 /* Assumes w==0, implementation according to C99 G.5.1.8 */
 {
    int sign = MPFR_SIGNBIT (mpc_realref (w));
    mpfr_t infty;
+
+   mpfr_init2 (infty, MPFR_PREC_MIN);
    mpfr_set_inf (infty, sign);
    mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd));
    mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd));
+   mpfr_clear (infty);
    return MPC_INEX (0, 0); /* exact */
 }
 
+/* this routine deals with the case where z is infinite and w finite */
 static int
 mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
 /* Assumes w finite and non-zero and z infinite; implementation
@@ -42,9 +47,13 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
    a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0);
    b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0);
 
+   /* a is -1 if Re(z) = -Inf, 1 if Re(z) = +Inf, 0 if Re(z) is finite
+      b is -1 if Im(z) = -Inf, 1 if Im(z) = +Inf, 0 if Im(z) is finite */
+
    /* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */
    /* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */
    if (a == 0 || b == 0) {
+     /* only one of a or b can be zero, since z is infinite */
       x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w));
       y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w));
    }
@@ -53,8 +62,9 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
          considerations and comparisons. Since operations with non-finite
          numbers are not considered time-critical, we let mpfr do the work. */
       mpfr_t sign;
+
       mpfr_init2 (sign, 2);
-         /* This is enough to determine the sign of sums and differences. */
+      /* This is enough to determine the sign of sums and differences. */
 
       if (a == 1)
          if (b == 1) {
@@ -98,6 +108,7 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
 }
 
 
+/* this routine deals with the case where z if finite and w infinite */
 static int
 mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
 /* Assumes z finite and w infinite; implementation according to
@@ -230,6 +241,7 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
    int underflow_norm, overflow_norm, underflow_prod, overflow_prod;
    int underflow_re, overflow_re, underflow_im = 0, overflow_im = 0;
    mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd);
+   int saved_underflow, saved_overflow;
 
    /* According to the C standard G.3, there are three types of numbers:   */
    /* finite (both parts are usual real numbers; contains 0), infinite     */
@@ -265,6 +277,10 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
    mpc_imagref (c_conj)[0] = mpc_imagref (c)[0];
    MPFR_CHANGE_SIGN (mpc_imagref (c_conj));
 
+   /* save the underflow or overflow flags from MPFR */
+   saved_underflow = mpfr_underflow_p ();
+   saved_overflow = mpfr_overflow_p ();
+
    do {
       loops ++;
       prec += loops <= 2 ? mpc_ceil_log2 (prec) + 5 : prec / 2;
@@ -388,5 +404,11 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
    mpc_clear (res);
    mpfr_clear (q);
 
+   /* restore underflow and overflow flags from MPFR */
+   if (saved_underflow)
+     mpfr_set_underflow ();
+   if (saved_overflow)
+     mpfr_set_overflow ();
+
    return (MPC_INEX (inexact_re, inexact_im));
 }
diff --git a/src/exp.c b/src/exp.c
index f427474..3646225 100644
--- a/src/exp.c
+++ b/src/exp.c
@@ -27,6 +27,7 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   mpfr_prec_t prec;
   int ok = 0;
   int inex_re, inex_im;
+  int saved_underflow, saved_overflow;
 
   /* special values */
   if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
@@ -142,6 +143,10 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   mpfr_init2 (y, 2);
   mpfr_init2 (z, 2);
 
+  /* save the underflow or overflow flags from MPFR */
+  saved_underflow = mpfr_underflow_p ();
+  saved_overflow = mpfr_overflow_p ();
+
   do
     {
       prec += mpc_ceil_log2 (prec) + 5;
@@ -187,5 +192,11 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   mpfr_clear (y);
   mpfr_clear (z);
 
+  /* restore underflow and overflow flags from MPFR */
+  if (saved_underflow)
+    mpfr_set_underflow ();
+  if (saved_overflow)
+    mpfr_set_overflow ();
+
   return MPC_INEX(inex_re, inex_im);
 }
diff --git a/src/norm.c b/src/norm.c
index 64e0f4a..ab413b6 100644
--- a/src/norm.c
+++ b/src/norm.c
@@ -27,6 +27,7 @@ int
 mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd)
 {
    int inexact;
+   int saved_underflow, saved_overflow;
 
    /* handling of special values; consistent with abs in that
       norm = abs^2; so norm (+-inf, xxx) = norm (xxx, +-inf) = +inf */
@@ -54,6 +55,10 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd)
       mpfr_init (v);
       mpfr_init (res);
 
+      /* save the underflow or overflow flags from MPFR */
+      saved_underflow = mpfr_underflow_p ();
+      saved_overflow = mpfr_overflow_p ();
+
       loops = 0;
       mpfr_clear_underflow ();
       mpfr_clear_overflow ();
@@ -162,6 +167,12 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd)
       else /* no problems, ternary value due to mpfr_can_round trick */
          inexact = mpfr_set (a, res, rnd);
 
+      /* restore underflow and overflow flags from MPFR */
+      if (saved_underflow)
+        mpfr_set_underflow ();
+      if (saved_overflow)
+        mpfr_set_overflow ();
+
       mpfr_clear (u);
       mpfr_clear (v);
       mpfr_clear (res);
diff --git a/src/pow.c b/src/pow.c
index 3c4e8dc..92abca7 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -483,6 +483,11 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
   int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0;
   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 ();
 
   x_real = mpfr_zero_p (mpc_imagref(x));
   y_real = mpfr_zero_p (mpc_imagref(y));
@@ -674,7 +679,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
          /* under- and overflow flags are set by mpc_exp */
          mpc_set (z, u, MPC_RNDNN);
          ret = ret_exp;
-         goto clear_t_and_u;
+         goto exact;
       }
 
       /* Since the error bound is global, we have to take into account the
@@ -804,12 +809,12 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
   mpc_clear (t);
   mpc_clear (u);
 
- end:
-  return ret;
-
+  /* restore underflow and overflow flags from MPFR */
+  if (saved_underflow)
+    mpfr_set_underflow ();
+  if (saved_overflow)
+    mpfr_set_overflow ();
 
-clear_t_and_u:
-  mpc_clear (t);
-  mpc_clear (u);
+ end:
   return ret;
 }
diff --git a/tests/div.dat b/tests/div.dat
index 8bbebf7..4617dde 100644
--- a/tests/div.dat
+++ b/tests/div.dat
@@ -2457,13 +2457,13 @@
 #0 0 10 1 10 0  10 0b1@536870912 10 0b1@-536870913 10 0b1@536870912 10 0b1@-536870913  N N
 # cases that should yield i, but cannot be handled due to intermediate
 # over- or underflows
-# current result: (0 @NaN@)
+# current result: (0 @Inf@)
 #0 0 10 0 10 1  10 -0b1@536870912 10 1 10 1 10 0b1@536870912  N N
-# current result: (@NaN@ @Inf@)
+# current result: (@NaN@ 1.0)
 #0 0 10 0 10 1  10 -0b1@-536870913 10 1 10 1 10 0b1@-536870913  N N
-# current result: (0 @NaN@)
+# current result: (0 @Inf@)
 #0 0 10 0 10 1  10 -0b1@536870912 10 0b1@536870912 10 0b1@536870912 10 0b1@536870912  N N
-# current result: (@NaN@ @NaN@)
+# current result: (@NaN@ 0)
 #0 0 10 0 10 1  10 -0b1@-536870913 10 0b1@-536870913 10 0b1@-536870913 10 0b1@-536870913  N N
 # current result: (@NaN@ @Inf@)
 #0 0 10 0 10 1  10 -0b1@-536870913 10 0b1@536870912 10 0b1@536870912 10 0b1@-536870913  N N
diff --git a/tests/read_data.c b/tests/read_data.c
index a6e1b64..ccfd0ed 100644
--- a/tests/read_data.c
+++ b/tests/read_data.c
@@ -554,6 +554,100 @@ read_ccs (FILE *fp, int *inex_re, int *inex_im, mpc_ptr expected,
   check_compatible (*inex_im, mpc_imagref(expected), MPC_RND_IM(*rnd), "imag");
 }
 
+/* set MPFR flags to random values */
+static void
+set_mpfr_flags (int counter)
+{
+  if (counter & 1)
+    mpfr_set_underflow ();
+  else
+    mpfr_clear_underflow ();
+  if (counter & 2)
+    mpfr_set_overflow ();
+  else
+    mpfr_clear_overflow ();
+  /* the divide-by-0 flag was added in MPFR 3.1.0 */
+#ifdef mpfr_set_divby0
+  if (counter & 4)
+    mpfr_set_divby0 ();
+  else
+    mpfr_clear_divby0 ();
+#endif
+  if (counter & 8)
+    mpfr_set_nanflag ();
+  else
+    mpfr_clear_nanflag ();
+  if (counter & 16)
+    mpfr_set_inexflag ();
+  else
+    mpfr_clear_inexflag ();
+  if (counter & 32)
+    mpfr_set_erangeflag ();
+  else
+    mpfr_clear_erangeflag ();
+}
+
+/* Check MPFR flags: we allow that some flags are set internally by MPC,
+   for example if MPC does internal computations (using MPFR) which yield
+   an overflow, even if the final MPC result fits in the exponent range.
+   However we don't allow MPC to *clear* the MPFR flags */
+static void
+check_mpfr_flags (int counter)
+{
+  int old, new;
+
+  old = (counter & 1) != 0;
+  new = mpfr_underflow_p () != 0;
+  if (old && (new == 0))
+    {
+      printf ("Error, underflow flag has been modified from %d to %d\n",
+              old, new);
+      exit (1);
+    }
+  old = (counter & 2) != 0;
+  new = mpfr_overflow_p () != 0;
+  if (old && (new == 0))
+    {
+      printf ("Error, overflow flag has been modified from %d to %d\n",
+              old, new);
+      exit (1);
+    }
+#ifdef mpfr_divby0_p
+  old = (counter & 4) != 0;
+  new = mpfr_divby0_p () != 0;
+  if (old && (new == 0))
+    {
+      printf ("Error, divby0 flag has been modified from %d to %d\n",
+              old, new);
+      exit (1);
+    }
+#endif
+  old = (counter & 8) != 0;
+  new = mpfr_nanflag_p () != 0;
+  if (old && (new == 0))
+    {
+      printf ("Error, nanflag flag has been modified from %d to %d\n",
+              old, new);
+      exit (1);
+    }
+  old = (counter & 16) != 0;
+  new = mpfr_inexflag_p () != 0;
+  if (old && (new == 0))
+    {
+      printf ("Error, inexflag flag has been modified from %d to %d\n",
+              old, new);
+      exit (1);
+    }
+  old = (counter & 32) != 0;
+  new = mpfr_erangeflag_p () != 0;
+  if (old && (new == 0))
+    {
+      printf ("Error, erangeflag flag has been modified from %d to %d\n",
+              old, new);
+      exit (1);
+    }
+}
+
 /* data_check (function, data_file_name) checks function results against
  precomputed data in a file.*/
 void
@@ -576,6 +670,8 @@ data_check (mpc_function function, const char *file_name)
   known_signs_t signs;
   int inex = 0;
 
+  static int rand_counter = 0;
+
   fp = open_data_file (file_name);
 
   /* 1. init needed variables */
@@ -615,6 +711,8 @@ data_check (mpc_function function, const char *file_name)
   nextchar = getc (fp);
   skip_whitespace_comments (fp);
   while (nextchar != EOF) {
+      set_mpfr_flags (rand_counter);
+
       /* for each kind of function prototype: */
       /* 3.1 read a line of data: expected result, parameters, rounding mode */
       /* 3.2 compute function at the same precision as the expected result */
@@ -919,6 +1017,10 @@ data_check (mpc_function function, const char *file_name)
           printf ("Unhandled function prototype %i in 'data_check'\n", function.type);
           exit (1);
         }
+
+      /* check MPFR flags were not modified */
+      check_mpfr_flags (rand_counter);
+      rand_counter ++;
     }
 
   /* 3. Clear used variables */