sqr.c: simplification of Karatsuba; call naive function for real part in case

   of intermediate under-/overflow
   potentially slightly slower in these corner cases, but more likely to
   be correct


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

2 files changed, 13 insertions, 67 deletions

diff --git a/src/sqr.c b/src/sqr.c
index f71a6ad..f7fa27a 100644
--- a/src/sqr.c
+++ b/src/sqr.c
@@ -245,7 +245,6 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
          and underflows are also handled correctly. */
 
       inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd));
-
    }
    else {
       /* Karatsuba squaring: we compute the real part as (x+y)*(x-y) and the
@@ -268,42 +267,33 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
          /* The error is bounded above by 1 ulp.                             */
          /* We first let inexact be 1 if the real part is not computed       */
          /* exactly and determine the sign later.                            */
-         inexact = ROUND_AWAY (mpfr_add (u, x, mpc_imagref (op), MPFR_RNDA), u)
+         inexact =  ROUND_AWAY (mpfr_add (u, x, mpc_imagref (op), MPFR_RNDA), u)
                   | ROUND_AWAY (mpfr_sub (v, x, mpc_imagref (op), MPFR_RNDA), v);
 
          /* compute the real part as u*v, rounded away                    */
          /* determine also the sign of inex_re                            */
-         if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0)
-         {
+
+         if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0) {
             /* as we have rounded away, the result is exact */
             mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
             inex_re = 0;
             ok = 1;
          }
-         else if (mpfr_sgn (u) * mpfr_sgn (v) > 0)
-         {
-            inexact |= mpfr_mul (u, u, v, GMP_RNDU); /* error 5 */
-            /* checks that no overflow nor underflow occurs: since u*v > 0
-               and we round up, an overflow will give +Inf, but an underflow
-               will give 0.5*2^emin */
-            if (mpfr_get_exp (u) == emin) /* underflow */
-            {
-                inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd));
-                ok = 1;
-            }
-            else if (mpfr_inf_p (u))
-            {
-               /* let mpc_realref(rop) be a "correctly rounded overflow" */
-               inex_re = mpfr_set_ui_2exp (mpc_realref (rop), 1, emax,
-                                           MPC_RND_RE (rnd));
+         else {
+            mpfr_rnd_t rnd_away;
+               /* FIXME: can be replaced by MPFR_RNDA in mpfr >= 3 */
+            rnd_away = (mpfr_sgn (u) * mpfr_sgn (v) > 0 ? GMP_RNDU : GMP_RNDD);
+            inexact |= ROUND_AWAY (mpfr_mul (u, u, v, MPFR_RNDA), u); /* error 5 */
+            if (mpfr_get_exp (u) == emin || mpfr_inf_p (u)) {
+               /* under- or overflow */
+               inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd));
                ok = 1;
             }
             else {
                ok = (!inexact) | mpfr_can_round (u, prec - 3,
-                     GMP_RNDU, GMP_RNDZ,
+                     rnd_away, GMP_RNDZ,
                      MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == GMP_RNDN));
-               if (ok)
-               {
+               if (ok) {
                   inex_re = mpfr_set (mpc_realref (rop), u, MPC_RND_RE (rnd));
                   if (inex_re == 0)
                      /* remember that u was already rounded */
@@ -311,47 +301,6 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
                }
             }
          }
-         else
-         {
-            mpfr_rnd_t rnd_re = MPC_RND_RE (rnd);
-            inexact |= mpfr_mul (u, u, v, GMP_RNDD); /* error 5 */
-            /* if an overflow occurs: since u*v < 0 and we round down,
-               the result is -Inf or -MAXDBL */
-            if (mpfr_inf_p (u))
-               {
-               /* replace by "correctly rounded overflow" */
-               mpfr_set_si (u, -1, GMP_RNDN);
-               mpfr_mul_2ui (u, u, mpfr_get_emax (), rnd_re);
-               ok = 1;
-               }
-            else
-               {
-               /* if an underflow happens (i.e., u = -0.5*2^emin since we round
-                  away from zero), the result will be an underflow */
-               if (mpfr_get_exp (u) == emin)
-                  {
-                     if (rnd_re == GMP_RNDZ || rnd_re == GMP_RNDN ||
-                        rnd_re == GMP_RNDU)
-                     {
-                        mpfr_set_ui (mpc_realref (rop), 0, rnd_re);
-                        inex_re = 1;
-                     }
-                     else /* round down or away from zero */ {
-                     mpfr_set (mpc_realref (rop), u, rnd_re);
-                     inex_re = -1;
-                     }
-                     break;
-                  }
-               ok = (!inexact) | mpfr_can_round (u, prec - 3, GMP_RNDD,
-                  GMP_RNDZ, MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == GMP_RNDN));
-               }
-            if (ok)
-               {
-               inex_re = mpfr_set (mpc_realref (rop), u, MPC_RND_RE (rnd));
-               if (inex_re == 0)
-                  inex_re = inexact;
-               }
-         }
       }
       while (!ok);
 
diff --git a/tests/sqr.dat b/tests/sqr.dat
index 095a9eb..72bfe07 100644
--- a/tests/sqr.dat
+++ b/tests/sqr.dat
@@ -167,7 +167,4 @@
 0 0 10 0 10 0b1e-1073741823  100 0b1@-536870912  100 0b1@-536870912 N N
 0 - 10 0 10 0  100 0b1@-536870913  100 0b1@-536870913 N N
 0 + 10 0 10 0b1@-1073741824  100 0b1@-536870913  100 0b1@-536870913 N U
-
-# wrong ternary value for real part with naive algorithm,
-# sqr.c:297: MPC assertion failed: mpfr_get_exp (u) != emin for Karatsuba
 + - 10 0b1e-1073741824 10 0  100 0b1@-536870912  100 0b1@-536870913 N N