mul.c: corrected uncleared variable from previous commit

mul.dat: test cases for overflow inside Karatsuba


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

3 files changed, 28 insertions, 5 deletions

diff --git a/TODO b/TODO
index 5618ceb..bd2401d 100644
--- a/TODO
+++ b/TODO
@@ -1,3 +1,8 @@
+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
+requires analysis of error propagation
+
 From Andreas Enge 30 August 2011:
 move static tests from tdiv, tpow_fr, tadd_fr, tpow_z, ttan, tui_div
 to data files
diff --git a/src/mul.c b/src/mul.c
index 7926488..e5a2274 100644
--- a/src/mul.c
+++ b/src/mul.c
@@ -413,22 +413,26 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
   /* now |a| >= |b| and |c| >= |d| */
   prec = MPC_MAX_PREC(rop);
 
-  mpfr_init2 (u, 2);
   mpfr_init2 (v, prec_v = mpfr_get_prec (a) + mpfr_get_prec (d));
   mpfr_init2 (w, prec_w = mpfr_get_prec (b) + mpfr_get_prec (c));
+  mpfr_init2 (u, 2);
   mpfr_init2 (x, 2);
 
   inexact = mpfr_mul (v, a, d, GMP_RNDN);
-  if (inexact)
+  if (inexact) {
      /* over- or underflow */
-     return mpc_mul_naive (rop, op1, op2, rnd);
+    ok = 0;
+    goto clear;
+  }
   if (mul_a == -1)
     mpfr_neg (v, v, GMP_RNDN);
 
   inexact = mpfr_mul (w, b, c, GMP_RNDN);
-  if (inexact)
+  if (inexact) {
      /* over- or underflow */
-     return mpc_mul_naive (rop, op1, op2, rnd);
+     ok = 0;
+     goto clear;
+  }
   if (mul_c == -1)
     mpfr_neg (w, w, GMP_RNDN);
 
@@ -565,6 +569,7 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
       mpc_set (rop, result, MPC_RNDNN);
    }
 
+clear:
    mpfr_clear (u);
    mpfr_clear (v);
    mpfr_clear (w);
diff --git a/tests/mul.dat b/tests/mul.dat
index 36b5b91..9ed14bd 100755
--- a/tests/mul.dat
+++ b/tests/mul.dat
@@ -156,6 +156,19 @@
 + - 10 -0b1.111111111e1073741822 10 0b1.111111111e1073741822  10 0x1@125000750 10 0x3@125000750  10 0x1@225000750 10 0x2@225000750 U D
 + - 10 -0   10 +0    10 0x1@-125000750 10 0x3@-125000750  10 0x1@-225000750 10 0x2@-225000750 U D
 - + 10 -0b1e-1073741824   10 0b1e-1073741824    10 0x1@-125000750 10 0x3@-125000750  10 0x1@-225000750 10 0x2@-225000750 D U
+# starting with Karatsuba
+- + 10000 -inf 100 +inf  100 0x1@125000750 100 0x3@125000750  100 0x1@225000750 100 0x2@225000750 N N
+- + 10000 -inf 100 +inf  100 0x1@225000750 100 0x2@225000750  100 0x1@125000750 100 0x3@125000750 N N
+- - 10000 -inf 100 -inf  100 0x1@225000750 100 -0x2@225000750  100 0x1@125000750 100 -0x3@125000750 N N
++ - 10000 -0   100 +0    100 0x1@-125000750 100 0x3@-125000750  100 0x1@-225000750 100 0x2@-225000750 N N
++ - 10000 -0   100 +0    100 0x1@-225000750 100 0x2@-225000750  100 0x1@-125000750 100 0x3@-125000750 N N
+- - 10000 +0   100 +0    100 0x3@-125000750 100 0x1@-125000750  100 0x1@-225000750 100 0x2@-225000750 N N
+- - 10000 +0   100 +0    100 0x4@-125000750 100 0x1@-125000750  100 0x1@-225000750 100 0x2@-225000750 N N
++ - 10000 -0   100 +0    100 0x2@-225000750 100 0x1@-225000750  100 0x1@-125000750 100 0x3@-125000750 N N
+0 - 10000 +0   100 +0    100 0x1@-225000750 100 0x1@-225000750  100 0x1@-125000750 100 0x1@-125000750 N N
+0 + 10000 +0   100 +inf  100 0x1@125000750 100 0x1@125000750  100 0x1@225000750 100 0x1@225000750 N N
++ 0 10000 +inf 100 +0    100 0x1@125000750 100 0x1@125000750  100 0x1@225000750 100 -0x1@225000750 N N
++ + 10000 +inf 100 +inf  100 0x1@125000750 100 0x3@125000750  100 0x1@143434706 100 0x2@143434705 N N
 
 # improve code coverage: case where sign_x==0 in mpc_mul_karatsuba
 0 0 2000 6 2000 8  2000 4 2000 2 2000 2 2000 1 N N