mul.c: reverted previous, too hasty commit; needs more thought

mpc-impl.h, mul.c, sqr.c: for the time being, exported mpfr_fmma to be shared
   between mul.c and sqr.c
sqr.dat: added commented out test that currently fails with naive and
   Karatsuba algorithms alike


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

3 files changed, 8 insertions, 151 deletions

diff --git a/src/mpc-impl.h b/src/mpc-impl.h
index 8080479..874902b 100644
--- a/src/mpc-impl.h
+++ b/src/mpc-impl.h
@@ -1,6 +1,6 @@
 /* mpc-impl.h -- Internal include file for mpc.
 
-Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011 INRIA
+Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA
 
 This file is part of GNU MPC.
 
@@ -149,6 +149,8 @@ __MPC_DECLSPEC int mpfr_regular_p __MPC_PROTO ((mpfr_srcptr));
 extern "C" {
 #endif
 
+
+__MPC_DECLSPEC int mpfr_fmma __MPC_PROTO ((mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_srcptr d, int sign, mpfr_rnd_t rnd));
 __MPC_DECLSPEC int  mpc_mul_naive __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t));
 __MPC_DECLSPEC int  mpc_mul_karatsuba __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t));
 __MPC_DECLSPEC int  mpc_pow_usi __MPC_PROTO ((mpc_ptr, mpc_srcptr, unsigned long, int, mpc_rnd_t));
diff --git a/src/mul.c b/src/mul.c
index 6d64329..842d0ef 100644
--- a/src/mul.c
+++ b/src/mul.c
@@ -171,7 +171,7 @@ mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
 }
 
 
-static int
+int
 mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
            mpfr_srcptr d, int sign, mpfr_rnd_t rnd)
 {
@@ -201,9 +201,10 @@ mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
       mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN);
       inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd);
    }
-   else if (mpfr_inf_p (u) || mpfr_inf_p (v)
-            || mpfr_zero_p (u) || mpfr_zero_p (v)) {
-      /* There is at least one over- or underflow.
+   else if (mpfr_nan_p (z) || (mpfr_zero_p (u) && mpfr_zero_p (v))) {
+      /* In the first case, u and v are infinities with opposite signs.
+         In the second case, u and v are zeroes; their sum may be 0 or the
+         least representable number, with a sign to be determined.
          Redo the computations with mpz_t exponents */
       mpfr_exp_t ea, eb, ec, ed;
       mpz_t eu, ev;
diff --git a/src/sqr.c b/src/sqr.c
index b5be792..a1699e3 100644
--- a/src/sqr.c
+++ b/src/sqr.c
@@ -21,152 +21,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 #include <stdio.h>    /* for MPC_ASSERT */
 #include "mpc-impl.h"
 
-#define mpz_add_si(z,x,y) do { \
-   if (y >= 0) \
-      mpz_add_ui (z, x, (long int) y); \
-   else \
-      mpz_sub_ui (z, x, (long int) (-y)); \
-   } while (0);
-
-static int
-mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
-           mpfr_srcptr d, int sign, mpfr_rnd_t rnd)
-{
-   /* Computes z = ab+cd if sign >= 0, or z = ab-cd if sign < 0.
-      Assumes that a, b, c, d are finite and non-zero; so any multiplication
-      of two of them yielding an infinity is an overflow, and a
-      multiplication yielding 0 is an underflow.
-      Assumes further that z is distinct from a, b, c, d. */
-
-   int inex;
-   mpfr_t u, v;
-
-   /* u=a*b, v=sign*c*d exactly */
-   mpfr_init2 (u, mpfr_get_prec (a) + mpfr_get_prec (b));
-   mpfr_init2 (v, mpfr_get_prec (c) + mpfr_get_prec (d));
-   mpfr_mul (u, a, b, GMP_RNDN);
-   mpfr_mul (v, c, d, GMP_RNDN);
-   if (sign < 0)
-      mpfr_neg (v, v, GMP_RNDN);
-
-   /* tentatively compute z as u+v; here we need z to be distinct
-      from a, b, c, d to not lose the latter */
-   inex = mpfr_add (z, u, v, rnd);
-
-   if (mpfr_inf_p (z)) {
-      /* replace by "correctly rounded overflow" */
-      mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN);
-      inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd);
-   }
-   else if (mpfr_inf_p (u) || mpfr_inf_p (v)
-            || mpfr_zero_p (u) || mpfr_zero_p (v)) {
-      /* There is at least one over- or underflow.
-         Redo the computations with mpz_t exponents */
-      mpfr_exp_t ea, eb, ec, ed;
-      mpz_t eu, ev;
-         /* cheat to work around the const qualifiers */
-
-      /* Normalise the input by shifting and keep track of the shifts in
-         the exponents of u and v */
-      ea = mpfr_get_exp (a);
-      eb = mpfr_get_exp (b);
-      ec = mpfr_get_exp (c);
-      ed = mpfr_get_exp (d);
-
-      mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0);
-      mpfr_set_exp ((mpfr_ptr) b, (mpfr_prec_t) 0);
-      mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0);
-      mpfr_set_exp ((mpfr_ptr) d, (mpfr_prec_t) 0);
-
-      mpz_init (eu);
-      mpz_init (ev);
-      mpz_set_si (eu, (long int) ea);
-      mpz_add_si (eu, eu, (long int) eb);
-      mpz_set_si (ev, (long int) ec);
-      mpz_add_si (ev, ev, (long int) ed);
-
-      /* recompute u and v and move exponents to eu and ev */
-      mpfr_mul (u, a, b, GMP_RNDN);
-      /* exponent of u is non-positive */
-      mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u)));
-      mpfr_set_exp (u, (mpfr_prec_t) 0);
-      mpfr_mul (v, c, d, GMP_RNDN);
-      if (sign < 0)
-         mpfr_neg (v, v, GMP_RNDN);
-      mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v)));
-      mpfr_set_exp (v, (mpfr_prec_t) 0);
-
-      if (mpfr_nan_p (z)) {
-         mpfr_exp_t emax = mpfr_get_emax ();
-         int overflow;
-         /* We have a = ma * 2^ea with 1/2 <= |ma| < 1 and ea <= emax, and
-            analogously for b. So eu <= 2*emax, and eu > emax since we have
-            an overflow. The same holds for ev. Shift u and v by as much as
-            possible so that one of them has exponent emax and the
-            remaining exponents in eu and ev are the same. Then carry out
-            the addition. Shifting u and v prevents an underflow. */
-         if (mpz_cmp (eu, ev) >= 0) {
-            mpfr_set_exp (u, emax);
-            mpz_sub_ui (eu, eu, (long int) emax);
-            mpz_sub (ev, ev, eu);
-            mpfr_set_exp (v, (mpfr_exp_t) mpz_get_ui (ev));
-               /* remaining common exponent is now in eu */
-         }
-         else {
-            mpfr_set_exp (v, emax);
-            mpz_sub_ui (ev, ev, (long int) emax);
-            mpz_sub (eu, eu, ev);
-            mpfr_set_exp (u, (mpfr_exp_t) mpz_get_ui (eu));
-            mpz_set (eu, ev);
-               /* remaining common exponent is now also in eu */
-         }
-         inex = mpfr_add (z, u, v, rnd);
-            /* Result is finite since u and v have different signs. */
-         overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd);
-         if (overflow)
-            inex = overflow;
-      }
-      else {
-         int underflow;
-         /* Addition of two zeroes with same sign. We have a = ma * 2^ea
-            with 1/2 <= |ma| < 1 and ea >= emin and similarly for b.
-            So 2*emin < 2*emin+1 <= eu < emin < 0, and analogously for v. */
-         mpfr_exp_t emin = mpfr_get_emin ();
-         if (mpz_cmp (eu, ev) <= 0) {
-            mpfr_set_exp (u, emin);
-            mpz_add_ui (eu, eu, (unsigned long int) (-emin));
-            mpz_sub (ev, ev, eu);
-            mpfr_set_exp (v, (mpfr_exp_t) mpz_get_si (ev));
-         }
-         else {
-            mpfr_set_exp (v, emin);
-            mpz_add_ui (ev, ev, (unsigned long int) (-emin));
-            mpz_sub (eu, eu, ev);
-            mpfr_set_exp (u, (mpfr_exp_t) mpz_get_si (eu));
-            mpz_set (eu, ev);
-         }
-         inex = mpfr_add (z, u, v, rnd);
-         mpz_neg (eu, eu);
-         underflow = mpfr_div_2ui (z, z, mpz_get_ui (eu), rnd);
-         if (underflow)
-            inex = underflow;
-      }
-
-      mpz_clear (eu);
-      mpz_clear (ev);
-
-      mpfr_set_exp ((mpfr_ptr) a, ea);
-      mpfr_set_exp ((mpfr_ptr) b, eb);
-      mpfr_set_exp ((mpfr_ptr) c, ec);
-      mpfr_set_exp ((mpfr_ptr) d, ed);
-         /* works also when some of a, b, c, d are not all distinct */
-   }
-
-   mpfr_clear (u);
-   mpfr_clear (v);
-
-   return inex;
-}
 
 int
 mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)