/* mpc_mul -- Multiply two complex numbers

Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012, 2016 INRIA

This file is part of GNU MPC.

GNU MPC is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.

You should have received a copy of the GNU Lesser General Public License
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);

/* compute z=x*y when x has an infinite part */
static int
mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y)
{
   /* Let x=xr+i*xi and y=yr+i*yi; extract the signs of the operands */
   int xrs = mpfr_signbit (mpc_realref (x)) ? -1 : 1;
   int xis = mpfr_signbit (mpc_imagref (x)) ? -1 : 1;
   int yrs = mpfr_signbit (mpc_realref (y)) ? -1 : 1;
   int yis = mpfr_signbit (mpc_imagref (y)) ? -1 : 1;

   int u, v;

   /* compute the sign of
      u = xrs * yrs * xr * yr - xis * yis * xi * yi
      v = xrs * yis * xr * yi + xis * yrs * xi * yr
      +1 if positive, -1 if negative, 0 if NaN */
   if (  mpfr_nan_p (mpc_realref (x)) || mpfr_nan_p (mpc_imagref (x))
      || mpfr_nan_p (mpc_realref (y)) || mpfr_nan_p (mpc_imagref (y))) {
      u = 0;
      v = 0;
   }
   else if (mpfr_inf_p (mpc_realref (x))) {
      /* x = (+/-inf) xr + i*xi */
      u = (   mpfr_zero_p (mpc_realref (y))
           || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_imagref (y)))
           || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_imagref (y)))
           || (   (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (y)))
              && xrs*yrs == xis*yis)
           ? 0 : xrs * yrs);
      v = (   mpfr_zero_p (mpc_imagref (y))
           || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_realref (y)))
           || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_realref (y)))
           || (   (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (x)))
               && xrs*yis != xis*yrs)
           ? 0 : xrs * yis);
   }
   else {
      /* x = xr + i*(+/-inf) with |xr| != inf */
      u = (   mpfr_zero_p (mpc_imagref (y))
           || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_realref (y)))
           || (mpfr_inf_p (mpc_realref (y)) && xrs*yrs == xis*yis)
           ? 0 : -xis * yis);
      v = (   mpfr_zero_p (mpc_realref (y))
           || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_imagref (y)))
           || (mpfr_inf_p (mpc_imagref (y)) && xrs*yis != xis*yrs)
           ? 0 : xis * yrs);
   }

   if (u == 0 && v == 0) {
      /* Naive result is NaN+i*NaN. Obtain an infinity using the algorithm
         given in Annex G.5.1 of the ISO C99 standard */
      int xr = (mpfr_zero_p (mpc_realref (x)) || mpfr_nan_p (mpc_realref (x)) ? 0
                : (mpfr_inf_p (mpc_realref (x)) ? 1 : 0));
      int xi = (mpfr_zero_p (mpc_imagref (x)) || mpfr_nan_p (mpc_imagref (x)) ? 0
                : (mpfr_inf_p (mpc_imagref (x)) ? 1 : 0));
      int yr = (mpfr_zero_p (mpc_realref (y)) || mpfr_nan_p (mpc_realref (y)) ? 0 : 1);
      int yi = (mpfr_zero_p (mpc_imagref (y)) || mpfr_nan_p (mpc_imagref (y)) ? 0 : 1);
      if (mpc_inf_p (y)) {
         yr = mpfr_inf_p (mpc_realref (y)) ? 1 : 0;
         yi = mpfr_inf_p (mpc_imagref (y)) ? 1 : 0;
      }

      u = xrs * xr * yrs * yr - xis * xi * yis * yi;
      v = xrs * xr * yis * yi + xis * xi * yrs * yr;
   }

   if (u == 0)
      mpfr_set_nan (mpc_realref (z));
   else
      mpfr_set_inf (mpc_realref (z), u);

   if (v == 0)
      mpfr_set_nan (mpc_imagref (z));
   else
      mpfr_set_inf (mpc_imagref (z), v);

   return MPC_INEX (0, 0); /* exact */
}


/* compute z = x*y for Im(y) == 0 */
static int
mul_real (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
{
   int xrs, xis, yrs, yis;
   int inex;

   /* save signs of operands */
   xrs = MPFR_SIGNBIT (mpc_realref (x));
   xis = MPFR_SIGNBIT (mpc_imagref (x));
   yrs = MPFR_SIGNBIT (mpc_realref (y));
   yis = MPFR_SIGNBIT (mpc_imagref (y));

   inex = mpc_mul_fr (z, x, mpc_realref (y), rnd);
   /* Signs of zeroes may be wrong. Their correction does not change the
      inexact flag. */
   if (mpfr_zero_p (mpc_realref (z)))
      mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == MPFR_RNDD
                     || (xrs != yrs && xis == yis), MPFR_RNDN);
   if (mpfr_zero_p (mpc_imagref (z)))
      mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == MPFR_RNDD
                     || (xrs != yis && xis != yrs), MPFR_RNDN);

   return inex;
}


/* compute z = x*y for Re(y) == 0, and Im(x) != 0 and Im(y) != 0 */
static int
mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
{
   int sign;
   int inex_re, inex_im;
   int overlap = z == x || z == y;
   mpc_t rop;

   if (overlap)
      mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z));
   else
      rop [0] = z[0];

   sign =    (MPFR_SIGNBIT (mpc_realref (y)) != MPFR_SIGNBIT (mpc_imagref (x)))
          && (MPFR_SIGNBIT (mpc_imagref (y)) != MPFR_SIGNBIT (mpc_realref (x)));

   inex_re = -mpfr_mul (mpc_realref (rop), mpc_imagref (x), mpc_imagref (y),
                        INV_RND (MPC_RND_RE (rnd)));
   mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); /* exact */
   inex_im = mpfr_mul (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y),
                       MPC_RND_IM (rnd));
   mpc_set (z, rop, MPC_RNDNN);

   /* Sign of zeroes may be wrong (note that Re(z) cannot be zero) */
   if (mpfr_zero_p (mpc_imagref (z)))
      mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == MPFR_RNDD
                     || sign, MPFR_RNDN);

   if (overlap)
      mpc_clear (rop);

   return MPC_INEX (inex_re, inex_im);
}

#define MPFR_MANT(x)      ((x)->_mpfr_d)
#define MPFR_PREC(x)      ((x)->_mpfr_prec)
#define MPFR_EXP(x)       ((x)->_mpfr_exp)
#define MPFR_LIMB_SIZE(x) ((MPFR_PREC (x) - 1) / GMP_NUMB_BITS + 1)

#if HAVE_MPFR_FMMA == 0
static int
mpc_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;
   mp_size_t an, bn, cn, dn;

   /* u=a*b, v=sign*c*d exactly */
   an = MPFR_LIMB_SIZE(a);
   bn = MPFR_LIMB_SIZE(b);
   cn = MPFR_LIMB_SIZE(c);
   dn = MPFR_LIMB_SIZE(d);
   MPFR_MANT(u) = malloc ((an + bn) * sizeof (mp_limb_t));
   MPFR_MANT(v) = malloc ((cn + dn) * sizeof (mp_limb_t));
   if (an >= bn)
     mpn_mul (MPFR_MANT(u), MPFR_MANT(a), an, MPFR_MANT(b), bn);
   else
     mpn_mul (MPFR_MANT(u), MPFR_MANT(b), bn, MPFR_MANT(a), an);
   if ((MPFR_MANT(u)[an + bn - 1] >> (GMP_NUMB_BITS - 1)) == 0)
     {
       mpn_lshift (MPFR_MANT(u), MPFR_MANT(u), an + bn, 1);
       MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b) - 1;
     }
   else
     MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b);
   if (cn >= dn)
     mpn_mul (MPFR_MANT(v), MPFR_MANT(c), cn, MPFR_MANT(d), dn);
   else
     mpn_mul (MPFR_MANT(v), MPFR_MANT(d), dn, MPFR_MANT(c), cn);
   if ((MPFR_MANT(v)[cn + dn - 1] >> (GMP_NUMB_BITS - 1)) == 0)
     {
       mpn_lshift (MPFR_MANT(v), MPFR_MANT(v), cn + dn, 1);
       MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d) - 1;
     }
   else
     MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d);
   MPFR_PREC(u) = (an + bn) * GMP_NUMB_BITS;
   MPFR_PREC(v) = (cn + dn) * GMP_NUMB_BITS;
   MPFR_SIGN(u) = MPFR_SIGN(a) * MPFR_SIGN(b);
   if (sign > 0)
     MPFR_SIGN(v) = MPFR_SIGN(c) * MPFR_SIGN(d);
   else
     MPFR_SIGN(v) = -MPFR_SIGN(c) * MPFR_SIGN(d);

   mpfr_check_range (u, 0, MPFR_RNDN);
   mpfr_check_range (v, 0, MPFR_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_regular_p(z) || !mpfr_regular_p(u) || !mpfr_regular_p(v))
     {
       if (mpfr_inf_p (z)) {
         /* replace by "correctly rounded overflow" */
         mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), MPFR_RNDN);
         inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd);
       }
       else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) {
         /* exactly u underflowed, determine inexact flag */
         inex = (mpfr_signbit (u) ? 1 : -1);
       }
       else if (mpfr_zero_p (v) && !mpfr_zero_p (u)) {
         /* exactly v underflowed, determine inexact flag */
         inex = (mpfr_signbit (v) ? 1 : -1);
       }
       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;
         /* 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, MPFR_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, MPFR_RNDN);
      if (sign < 0)
         mpfr_neg (v, v, MPFR_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 */
       }
     }

   free (MPFR_MANT(u));
   free (MPFR_MANT(v));
   return inex;
}
#endif

int
mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
{
   /* computes z=x*y by the schoolbook method, where x and y are assumed
      to be finite and without zero parts                                */
   int overlap, inex;
   mpc_t rop;

   MPC_ASSERT (   mpfr_regular_p (mpc_realref (x)) && mpfr_regular_p (mpc_imagref (x))
               && mpfr_regular_p (mpc_realref (y)) && mpfr_regular_p (mpc_imagref (y)));
   overlap = (z == x) || (z == y);
   if (overlap)
      mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z));
   else
      rop [0] = z [0];

#if HAVE_MPFR_FMMA
   inex = MPC_INEX (mpfr_fmms (mpc_realref (rop), mpc_realref (x), mpc_realref (y), mpc_imagref (x),
                               mpc_imagref (y), MPC_RND_RE (rnd)),
                    mpfr_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), mpc_imagref (x),
                               mpc_realref (y), MPC_RND_IM (rnd)));
#else
   inex = MPC_INEX (mpc_fmma (mpc_realref (rop), mpc_realref (x), mpc_realref (y), mpc_imagref (x),
                               mpc_imagref (y), -1, MPC_RND_RE (rnd)),
                    mpc_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), mpc_imagref (x),
                               mpc_realref (y), +1, MPC_RND_IM (rnd)));
#endif

   mpc_set (z, rop, MPC_RNDNN);
   if (overlap)
      mpc_clear (rop);

   return inex;
}


int
mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
{
   /* computes rop=op1*op2 by a Karatsuba algorithm, where op1 and op2
      are assumed to be finite and without zero parts                  */
  mpfr_srcptr a, b, c, d;
  int mul_i, ok, inexact, mul_a, mul_c, inex_re = 0, inex_im = 0, sign_x, sign_u;
  mpfr_t u, v, w, x;
  mpfr_prec_t prec, prec_re, prec_u, prec_v, prec_w;
  mpfr_rnd_t rnd_re, rnd_u;
  int overlap;
     /* true if rop == op1 or rop == op2 */
  mpc_t result;
     /* overlap is quite difficult to handle, because we have to tentatively
        round the variable u in the end to either the real or the imaginary
        part of rop (it is not possible to tell now whether the real or
        imaginary part is used). If this fails, we have to start again and
        need the correct values of op1 and op2.
        So we just create a new variable for the result in this case. */
  int loop;
  const int MAX_MUL_LOOP = 1;

  overlap = (rop == op1) || (rop == op2);
  if (overlap)
     mpc_init3 (result, MPC_PREC_RE (rop), MPC_PREC_IM (rop));
  else
     result [0] = rop [0];

  a = mpc_realref(op1);
  b = mpc_imagref(op1);
  c = mpc_realref(op2);
  d = mpc_imagref(op2);

  /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */

  mul_i = 0; /* number of multiplications by i */
  mul_a = 1; /* implicit factor for a */
  mul_c = 1; /* implicit factor for c */

  if (mpfr_cmp_abs (a, b) < 0)
    {
      MPFR_SWAP (a, b);
      mul_i ++;
      mul_a = -1; /* consider i * (a+i*b) = -b + i*a */
    }

  if (mpfr_cmp_abs (c, d) < 0)
    {
      MPFR_SWAP (c, d);
      mul_i ++;
      mul_c = -1; /* consider -d + i*c instead of c + i*d */
    }

  /* find the precision and rounding mode for the new real part */
  if (mul_i % 2)
    {
      prec_re = MPC_PREC_IM(rop);
      rnd_re = MPC_RND_IM(rnd);
    }
  else /* mul_i = 0 or 2 */
    {
      prec_re = MPC_PREC_RE(rop);
      rnd_re = MPC_RND_RE(rnd);
    }

  if (mul_i)
    rnd_re = INV_RND(rnd_re);

  /* now |a| >= |b| and |c| >= |d| */
  prec = MPC_MAX_PREC(rop);

  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, MPFR_RNDN);
  if (inexact) {
     /* over- or underflow */
    ok = 0;
    goto clear;
  }
  if (mul_a == -1)
    mpfr_neg (v, v, MPFR_RNDN);

  inexact = mpfr_mul (w, b, c, MPFR_RNDN);
  if (inexact) {
     /* over- or underflow */
     ok = 0;
     goto clear;
  }
  if (mul_c == -1)
    mpfr_neg (w, w, MPFR_RNDN);

  /* compute sign(v-w) */
  sign_x = mpfr_cmp_abs (v, w);
  if (sign_x > 0)
    sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w);
  else if (sign_x == 0)
    sign_x = mpfr_sgn (v) - mpfr_sgn (w);
  else
    sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w);

   sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c);

   if (sign_x * sign_u < 0)
    {  /* swap inputs */
      MPFR_SWAP (a, c);
      MPFR_SWAP (b, d);
      mpfr_swap (v, w);
      { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; }
      sign_x = - sign_x;
    }

   /* now sign_x * sign_u >= 0 */
   loop = 0;
   do
     {
        loop++;
         /* the following should give failures with prob. <= 1/prec */
         prec += mpc_ceil_log2 (prec) + 3;

         mpfr_set_prec (u, prec_u = prec);
         mpfr_set_prec (x, prec);

         /* first compute away(b +/- a) and store it in u */
         inexact = (mul_a == -1 ?
                    mpfr_sub (u, b, a, MPFR_RNDA) :
                    mpfr_add (u, b, a, MPFR_RNDA));

         /* then compute away(+/-c - d) and store it in x */
         inexact |= (mul_c == -1 ?
                     mpfr_add (x, c, d, MPFR_RNDA) :
                     mpfr_sub (x, c, d, MPFR_RNDA));
         if (mul_c == -1)
           mpfr_neg (x, x, MPFR_RNDN);

         if (inexact == 0)
            mpfr_prec_round (u, prec_u = 2 * prec, MPFR_RNDN);

         /* compute away(u*x) and store it in u */
         inexact |= mpfr_mul (u, u, x, MPFR_RNDA);
            /* (a+b)*(c-d) */

	 /* if all computations are exact up to here, it may be that
	    the real part is exact, thus we need if possible to
	    compute v - w exactly */
	 if (inexact == 0)
	   {
	     mpfr_prec_t prec_x;
             /* v and w are different from 0, so mpfr_get_exp is safe to use */
             prec_x = SAFE_ABS (mpfr_exp_t, mpfr_get_exp (v) - mpfr_get_exp (w))
                      + MPC_MAX (prec_v, prec_w) + 1;
                      /* +1 is necessary for a potential carry */
	     /* ensure we do not use a too large precision */
	     if (prec_x > prec_u)
               prec_x = prec_u;
	     if (prec_x > prec)
	       mpfr_prec_round (x, prec_x, MPFR_RNDN);
	   }

         rnd_u = (sign_u > 0) ? MPFR_RNDU : MPFR_RNDD;
         inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */

         /* in case u=0, ensure that rnd_u rounds x away from zero */
         if (mpfr_sgn (u) == 0)
           rnd_u = (mpfr_sgn (x) > 0) ? MPFR_RNDU : MPFR_RNDD;
         inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */

         ok = inexact == 0 ||
           mpfr_can_round (u, prec_u - 3, rnd_u, MPFR_RNDZ,
                           prec_re + (rnd_re == MPFR_RNDN));
         /* this ensures both we can round correctly and determine the correct
            inexact flag (for rounding to nearest) */
     }
   while (!ok && loop <= MAX_MUL_LOOP);
      /* after MAX_MUL_LOOP rounds, use mpc_naive instead */

   if (ok) {
      /* if inexact is zero, then u is exactly ac-bd, otherwise fix the sign
         of the inexact flag for u, which was rounded away from ac-bd */
      if (inexact != 0)
      inexact = mpfr_sgn (u);

      if (mul_i == 0)
      {
         inex_re = mpfr_set (mpc_realref(result), u, MPC_RND_RE(rnd));
         if (inex_re == 0)
            {
            inex_re = inexact; /* u is rounded away from 0 */
            inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd));
            }
         else
            inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd));
      }
      else if (mul_i == 1) /* (x+i*y)/i = y - i*x */
      {
         inex_im = mpfr_neg (mpc_imagref(result), u, MPC_RND_IM(rnd));
         if (inex_im == 0)
            {
            inex_im = -inexact; /* u is rounded away from 0 */
            inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd));
            }
         else
            inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd));
      }
      else /* mul_i = 2, z/i^2 = -z */
      {
         inex_re = mpfr_neg (mpc_realref(result), u, MPC_RND_RE(rnd));
         if (inex_re == 0)
            {
            inex_re = -inexact; /* u is rounded away from 0 */
            inex_im = -mpfr_add (mpc_imagref(result), v, w,
                                 INV_RND(MPC_RND_IM(rnd)));
            mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd));
            }
         else
            {
            inex_im = -mpfr_add (mpc_imagref(result), v, w,
                                 INV_RND(MPC_RND_IM(rnd)));
            mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd));
            }
      }

      mpc_set (rop, result, MPC_RNDNN);
   }

clear:
   mpfr_clear (u);
   mpfr_clear (v);
   mpfr_clear (w);
   mpfr_clear (x);
   if (overlap)
      mpc_clear (result);

   if (ok)
      return MPC_INEX(inex_re, inex_im);
   else
      return mpc_mul_naive (rop, op1, op2, rnd);
}


int
mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
{
   /* Conforming to ISO C99 standard (G.5.1 multiplicative operators),
      infinities are treated specially if both parts are NaN when computed
      naively. */
   if (mpc_inf_p (b))
      return mul_infinite (a, b, c);
   if (mpc_inf_p (c))
      return mul_infinite (a, c, b);

   /* NaN contamination of both parts in result */
   if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b))
       || mpfr_nan_p (mpc_realref (c)) || mpfr_nan_p (mpc_imagref (c))) {
      mpfr_set_nan (mpc_realref (a));
      mpfr_set_nan (mpc_imagref (a));
      return MPC_INEX (0, 0);
   }

   /* check for real multiplication */
   if (mpfr_zero_p (mpc_imagref (b)))
      return mul_real (a, c, b, rnd);
   if (mpfr_zero_p (mpc_imagref (c)))
      return mul_real (a, b, c, rnd);

   /* check for purely imaginary multiplication */
   if (mpfr_zero_p (mpc_realref (b)))
      return mul_imag (a, c, b, rnd);
   if (mpfr_zero_p (mpc_realref (c)))
      return mul_imag (a, b, c, rnd);

   /* If the real and imaginary part of one argument have a very different */
   /* exponent, it is not reasonable to use Karatsuba multiplication.      */
   if (   SAFE_ABS (mpfr_exp_t,
                     mpfr_get_exp (mpc_realref (b)) - mpfr_get_exp (mpc_imagref (b)))
         > (mpfr_exp_t) MPC_MAX_PREC (b) / 2
      || SAFE_ABS (mpfr_exp_t,
                     mpfr_get_exp (mpc_realref (c)) - mpfr_get_exp (mpc_imagref (c)))
         > (mpfr_exp_t) MPC_MAX_PREC (c) / 2)
      return mpc_mul_naive (a, b, c, rnd);
   else
      return ((MPC_MAX_PREC(a)
               <= (mpfr_prec_t) MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB)
            ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd);
}