div.c: implemented return value

div.dat: updated accordingly


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

5 files changed, 18 insertions, 25 deletions

diff --git a/NEWS b/NEWS
index 9053487..f894389 100644
--- a/NEWS
+++ b/NEWS
@@ -4,8 +4,10 @@ Changes in version 1.0:
   - Bug fixes:
     - mul and norm now return correct results even if there are over- or
       underflows during the computation
-    - asin, proj, sqr: wrong result when input variable has infinite part and
-      equals output variable
+    - div now returns a value indicating the effective rounding direction,
+      as the other functions.
+    - asin, proj, sqr: Wrong result when input variable has infinite part and
+      equals output variable is corrected.
 
 Changes in version 0.9:
   - New functions:
diff --git a/TODO b/TODO
index 9896e4e..1809a6d 100644
--- a/TODO
+++ b/TODO
@@ -1,9 +1,3 @@
-From Paul Zimmermann 1st September 2011:
-return the exact ternary value in mpc_div:
-- for directed rounding modes it is easy (assuming we are able to detect
-  exact cases)
-- rounding to nearest to n bits is equivalent to directed rounding to n+1 bits
-
 From Andreas Enge 31 August 2011:
 implement mul_karatsuba with three multiplications at precision around p,
 instead of two at precision 2*p and one at precision p
diff --git a/doc/mpc.texi b/doc/mpc.texi
index f8c6f42..8e400ef 100644
--- a/doc/mpc.texi
+++ b/doc/mpc.texi
@@ -881,11 +881,6 @@ Set @var{rop} to the square of @var{op} rounded according to @var{rnd}.
 @deftypefunx int mpc_ui_div (mpc_t @var{rop}, unsigned long int @var{op1}, mpc_t @var{op2}, mpc_rnd_t @var{rnd})
 @deftypefunx int mpc_fr_div (mpc_t @var{rop}, mpfr_t @var{op1}, mpc_t @var{op2}, mpc_rnd_t @var{rnd})
 Set @var{rop} to @var{op1}/@var{op2} rounded according to @var{rnd}.
-For @code{mpc_div}, @code{mpc_ui_div} and @code{mpc_fr_div}, the usual
-semantics of the return value is not yet implemented. We only guarantee
-that a return value of 0 indicates an exact result.
-But exact results may lead to a non-zero return value, and the sign of a
-non-zero return value is not significant.
 @end deftypefun
 
 @deftypefun int mpc_neg (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd})
@@ -902,7 +897,6 @@ if @var{rop} and @var{op} are the same variable.
 @deftypefun int mpc_abs (mpfr_t @var{rop}, mpc_t @var{op}, mpfr_rnd_t @var{rnd})
 Set the floating-point number @var{rop} to the absolute value of @var{op},
 rounded in the direction @var{rnd}.
-The returned value is zero iff the result is exact.
 @end deftypefun
 
 @deftypefun int mpc_norm (mpfr_t @var{rop}, mpc_t @var{op}, mpfr_rnd_t @var{rnd})
diff --git a/src/div.c b/src/div.c
index 2ed8585..ffb260b 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 INRIA
+Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011 INRIA
 
 This file is part of GNU MPC.
 
@@ -227,9 +227,10 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
    mpc_t res, c_conj;
    mpfr_t q;
    mpfr_prec_t prec;
-   int inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
+   int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
    int underflow_norm, overflow_norm;
    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);
 
    /* According to the C standard G.3, there are three types of numbers:   */
    /* finite (both parts are usual real numbers; contains 0), infinite     */
@@ -299,7 +300,7 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
          overflow_re = mpfr_overflow_p ();
          ok_re = !inexact_re || underflow_re || overflow_re
                  || mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDZ,
-                                    MPC_RND_RE(rnd), MPC_PREC_RE(a));
+                    GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN));
 
          if (ok_re) /* compute imaginary part */ {
             mpfr_clear_underflow ();
@@ -308,8 +309,8 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
             underflow_im = mpfr_underflow_p ();
             overflow_im = mpfr_overflow_p ();
             ok_im = !inexact_im || underflow_im || overflow_im
-                    || mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDZ,
-                                       MPC_RND_IM(rnd), MPC_PREC_IM(a));
+                    || mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDN,
+                       GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN));
          }
       }
       else {
@@ -331,7 +332,7 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
          overflow_re = mpfr_overflow_p ();
          ok_re = !inexact_re || underflow_re || overflow_re
                  || mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDZ,
-                                    MPC_RND_RE(rnd), MPC_PREC_RE(a));
+                    GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN));
 
          if (ok_re) /* compute imaginary part */ {
             mpfr_clear_underflow ();
@@ -342,14 +343,16 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
             overflow_im = mpfr_overflow_p ();
             ok_im = !inexact_im || underflow_im || overflow_im
                     || mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDZ,
-                                       MPC_RND_IM(rnd), MPC_PREC_IM(a));
+                       GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN));
          }
       }
    }
    while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm);
 
-   mpc_set (a, res, rnd);
-
+   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) {
       mpfr_set_inf (MPC_RE (a), mpfr_sgn (MPC_RE (res)));
diff --git a/tests/div.dat b/tests/div.dat
index 4d7c1b5..ff2e1f3 100644
--- a/tests/div.dat
+++ b/tests/div.dat
@@ -2428,13 +2428,13 @@
 - + 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
+- - 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
 0 0 2 2 2 -1   2 4 2 3  2 1 2 2 N N
 + - 4 1.375 4 1.25   4 15 4 14  4 11 4 0 N N
 
 # Bug 20080923
-#+ - 4 0b11@527 4 -0b111@-489   4 -0b11@-206 4 0  4 -0b1@-733 4 -0b101@-1750 N Z
++ + 4 0b11@527 4 -0b111@-489   4 -0b11@-206 4 0  4 -0b1@-733 4 -0b101@-1750 N Z
 
 # potential intermediate over- or underflow
 0 0 10 1 10 0  10 0 10 0b1@536870912 10 0 10 0b1@536870912  N N