/* 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);

}