div.c: some heuristics for intermediate under-/overflows

div.dat: added (and commented out) test cases that should give 1 or i, but
   fail due to intermediate under-/overflows; returned results are Nan or
   infinite
tui_div currently fails


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

2 files changed, 54 insertions, 11 deletions

diff --git a/src/div.c b/src/div.c
index 7d8c811..31d7656 100644
--- a/src/div.c
+++ b/src/div.c
@@ -18,7 +18,7 @@ You should have received a copy of the GNU Lesser General Public License
 along with this program. If not, see http://www.gnu.org/licenses/ .
 */
 
-#include <stdio.h>
+#include <stdio.h> /* for MPC_ASSERT */
 #include "mpc-impl.h"
 
 static int
@@ -228,7 +228,7 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
    mpfr_t q;
    mpfr_prec_t prec;
    int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
-   int underflow_norm, overflow_norm;
+   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);
 
@@ -284,9 +284,27 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
          /* to obtain divisions by 0 later on */
 
       /* now compute b*conjugate(c) */
+      mpfr_clear_underflow ();
+      mpfr_clear_overflow ();
       inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ);
       inexact_re = MPC_INEX_RE (inexact_prod);
       inexact_im = MPC_INEX_IM (inexact_prod);
+      underflow_prod = mpfr_underflow_p ();
+      overflow_prod = mpfr_overflow_p ();
+         /* unfortunately, does not distinguish between under-/overflow
+            in real or imaginary parts
+            hopefully, the side-effects of mpc_mul do indeed raise the
+            mpfr exceptions
+            FIXME: handling of under- and overflows here and in the
+            following is not totally rigorous and just a best effort to
+            salvage almost hopeless situations                          */
+      if (overflow_prod) {
+         if (!mpfr_zero_p (MPC_RE (res)))
+            mpfr_set_inf (MPC_RE (res), mpfr_sgn (MPC_RE (res)));
+         if (!mpfr_zero_p (MPC_IM (res)))
+            mpfr_set_inf (MPC_IM (res), mpfr_sgn (MPC_IM (res)));
+      }
+            
 
       /* divide the product by the norm */
       if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) {
@@ -316,7 +334,6 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
       else {
          /* The division is inexact, so for efficiency reasons we invert q */
          /* only once and multiply by the inverse. */
-         /* We do not decide about the sign of the difference. */
          if (mpfr_ui_div (q, 1ul, q, GMP_RNDZ) || inexact_norm) {
              /* if 1/q is inexact, the approximations of the real and
                 imaginary part below will be inexact, unless RE(res)
@@ -347,26 +364,28 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
          }
       }
    }
-   while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm);
+   while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm
+                             && !underflow_prod && !overflow_prod);
 
    inex = mpc_set (a, res, rnd);
    inexact_re = MPC_INEX_RE (inex);
    inexact_im = MPC_INEX_IM (inex);
    
    /* fix values and inexact flags in case of overflow/underflow */
-   if (overflow_re || underflow_norm) {
+   /* FIXME: heuristic, certainly does not cover all cases */
+   if (overflow_re || (underflow_norm && !underflow_prod)) {
       mpfr_set_inf (MPC_RE (a), mpfr_sgn (MPC_RE (res)));
       inexact_re = mpfr_sgn (MPC_RE (res));
    }
-   else if (underflow_re || overflow_norm) {
+   else if (underflow_re || (overflow_norm && !overflow_prod)) {
       mpfr_set_zero (MPC_RE (a), mpfr_sgn (MPC_RE (res)));
       inexact_re = -mpfr_sgn (MPC_RE (res));
    }
-   if (overflow_im || underflow_norm) {
+   if (overflow_im || (underflow_norm && !underflow_prod)) {
       mpfr_set_inf (MPC_IM (a), mpfr_sgn (MPC_IM (res)));
       inexact_im = mpfr_sgn (MPC_IM (res));
    }
-   else if (underflow_im || overflow_norm) {
+   else if (underflow_im || (overflow_norm && !overflow_prod)) {
       mpfr_set_zero (MPC_IM (a), mpfr_sgn (MPC_IM (res)));
       inexact_im = -mpfr_sgn (MPC_IM (res));
    }
diff --git a/tests/div.dat b/tests/div.dat
index ff2e1f3..2459b59 100644
--- a/tests/div.dat
+++ b/tests/div.dat
@@ -2424,9 +2424,6 @@
 0 0 7    1 7    1   7    1 7    1  7    1 7   +0 N N
 0 0 7    1 7   +0   7    1 7    1  7    1 7    1 N N
 
-# overflow (reported by Emmanuel Thome)
-- + 250 -inf 250 +inf 250 1 250 0 250 -1e-164895850 250 -1e-164895850 N N
-
 # small exact/inexact examples
 - - 10 0b1.01010001 10 -0b1.1000101@-6   10 973 10 964  10 725 10 745 N N
 0 0 10 -14 10 9   10 -837 10 637  10 63 10 -5 N N
@@ -2438,5 +2435,32 @@
 
 # potential intermediate over- or underflow
 0 0 10 1 10 0  10 0 10 0b1@536870912 10 0 10 0b1@536870912  N N
+0 0 10 1 10 0  10 0b1@536870912 10 0 10 0b1@536870912 10 0 N N
+
+# overflow (reported by Emmanuel Thome)
+- + 250 -inf 250 +inf 250 1 250 0 250 -1e-164895850 250 -1e-164895850 N N
+
+# cases that should yield 1, but cannot be handled due to intermediate
+# over- or underflows
+# current result: (@NaN@ 0)
 #0 0 10 1 10 0  10 1 10 0b1@536870912 10 1 10 0b1@536870912  N N
+# current result: (@Inf@ 0)
 #0 0 10 1 10 0  10 1 10 0b1@-536870913 10 1 10 0b1@-536870913  N N
+# current result: (@NaN@ 0)
+#0 0 10 1 10 0  10 0b1@536870912 10 0b1@536870912 10 0b1@536870912 10 0b1@536870912  N N
+# current result: (@NaN@ 0)
+#0 0 10 1 10 0  10 0b1@-536870913 10 0b1@-536870913 10 0b1@-536870913 10 0b1@-536870913  N N
+# current result: (@Inf@ 0)
+#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@)
+#0 0 10 0 10 1  10 -0b1@536870912 10 1 10 1 10 0b1@536870912  N N
+# current result: (@NaN@ @Inf@)
+#0 0 10 0 10 1  10 -0b1@-536870913 10 1 10 1 10 0b1@-536870913  N N
+# current result: (0 @NaN@)
+#0 0 10 0 10 1  10 -0b1@536870912 10 0b1@536870912 10 0b1@536870912 10 0b1@536870912  N N
+# current result: (@NaN@ @NaN@)
+#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