sqr: copied code for imaginary part from naive computation also to the

   Karatsuba algorithm
sqr.dat: added examples with imaginary part previously miscomputed by
   Karatsuba


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

2 files changed, 15 insertions, 21 deletions

diff --git a/src/sqr.c b/src/sqr.c
index 109e314..47a3f7e 100644
--- a/src/sqr.c
+++ b/src/sqr.c
@@ -180,6 +180,7 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
    mpfr_prec_t prec;
    int inex_re, inex_im, inexact;
    mpfr_exp_t emin, emax;
+   int saved_underflow;
 
    /* special values: NaN and infinities */
    if (!mpc_fin_p (op)) {
@@ -252,20 +253,10 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
          exactly with the standard formulae instead, even if this means an
          additional multiplication. Using the approach copied from mul, over-
          and underflows are also handled correctly. */
-      int saved_underflow;
 
       inex_re = mpfr_fmma (rop->re, x, x, op->im, op->im, -1,
                            MPC_RND_RE (rnd));
 
-      saved_underflow = mpfr_underflow_p ();
-      mpfr_clear_underflow ();
-      inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd));
-      if (!mpfr_underflow_p ())
-         inex_im |= mpfr_mul_2exp (rop->im, rop->im, 1, MPC_RND_IM (rnd));
-         /* We must not multiply by 2 if rop->im has been set to the smallest
-            representable number. */
-      if (saved_underflow)
-         mpfr_set_underflow ();
    }
    else {
       /* Karatsuba squaring */
@@ -367,19 +358,18 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 
       mpfr_clear (u);
       mpfr_clear (v);
-
-      if (mpfr_get_exp (x) + mpfr_get_exp(mpc_imagref (op)) <= emin + 1) {
-         mpfr_mul_2ui (x, x, 1, GMP_RNDN);
-         inex_im = mpfr_mul (mpc_imagref (rop), x, mpc_imagref (op), MPC_RND_IM (rnd));
-      }
-      else {
-         inex_im = mpfr_mul (mpc_imagref (rop), x, mpc_imagref (op), MPC_RND_IM (rnd));
-         /* checks that no underflow occurs (an overflow can occur here) */
-         MPC_ASSERT (mpfr_zero_p (mpc_imagref (rop)) == 0);
-         mpfr_mul_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
-      }
    }
 
+   saved_underflow = mpfr_underflow_p ();
+   mpfr_clear_underflow ();
+   inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd));
+   if (!mpfr_underflow_p ())
+      inex_im |= mpfr_mul_2exp (rop->im, rop->im, 1, MPC_RND_IM (rnd));
+      /* We must not multiply by 2 if rop->im has been set to the smallest
+         representable number. */
+   if (saved_underflow)
+      mpfr_set_underflow ();
+
    if (rop == op)
       mpfr_clear (x);
 
diff --git a/tests/sqr.dat b/tests/sqr.dat
index ee59b0e..9d56f0e 100644
--- a/tests/sqr.dat
+++ b/tests/sqr.dat
@@ -160,6 +160,10 @@
 # another interesting one with not exactly the same behaviour
 - - 100 -inf 100 -inf  100 -0xf@192058806 100 0x1@192058873 N N
 0 + 100 0 100 inf  100 0x1@192058806  100 0x1@192058806 N N
+# Re(op)*Im(op) can be computed, but multiplication by 2 triggers overflow
+0 + 100 0 100 inf  100 0b1@536870911  100 0b1@536870911 N N
+0 - 10 0 10 0b1.111111111e1073741822  100 0b1@536870911  100 0b1@536870911 N D
+0 - 10 0 10 0b1.111111111e1073741822  100 0b1@536870912  100 0b1@536870912 N D
 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