contrib/mpc/src/mul.c

   1 /* mpc_mul -- Multiply two complex numbers
   2
   3 Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA
   4
   5 This file is part of GNU MPC.
   6
   7 GNU MPC is free software; you can redistribute it and/or modify it under
   8 the terms of the GNU Lesser General Public License as published by the
   9 Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
  13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
  15 more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with this program. If not, see http://www.gnu.org/licenses/ .
  19 */
  20
  21 #include <stdio.h>    /* for MPC_ASSERT */
  22 #include "mpc-impl.h"
  23
  24 #define mpz_add_si(z,x,y) do { \
  25    if (y >= 0) \
  26       mpz_add_ui (z, x, (long int) y); \
  27    else \
  28       mpz_sub_ui (z, x, (long int) (-y)); \
  29    } while (0);
  30
  31 /* compute z=x*y when x has an infinite part */
  32 static int
  33 mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y)
  34 {
  35    /* Let x=xr+i*xi and y=yr+i*yi; extract the signs of the operands */
  36    int xrs = mpfr_signbit (mpc_realref (x)) ? -1 : 1;
  37    int xis = mpfr_signbit (mpc_imagref (x)) ? -1 : 1;
  38    int yrs = mpfr_signbit (mpc_realref (y)) ? -1 : 1;
  39    int yis = mpfr_signbit (mpc_imagref (y)) ? -1 : 1;
  40
  41    int u, v;
  42
  43    /* compute the sign of
  44       u = xrs * yrs * xr * yr - xis * yis * xi * yi
  45       v = xrs * yis * xr * yi + xis * yrs * xi * yr
  46       +1 if positive, -1 if negatiye, 0 if NaN */
  47    if (  mpfr_nan_p (mpc_realref (x)) || mpfr_nan_p (mpc_imagref (x))
  48       || mpfr_nan_p (mpc_realref (y)) || mpfr_nan_p (mpc_imagref (y))) {
  49       u = 0;
  50       v = 0;
  51    }
  52    else if (mpfr_inf_p (mpc_realref (x))) {
  53       /* x = (+/-inf) xr + i*xi */
  54       u = (   mpfr_zero_p (mpc_realref (y))
  55            || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_imagref (y)))
  56            || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_imagref (y)))
  57            || (   (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (y)))
  58               && xrs*yrs == xis*yis)
  59            ? 0 : xrs * yrs);
  60       v = (   mpfr_zero_p (mpc_imagref (y))
  61            || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_realref (y)))
  62            || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_realref (y)))
  63            || (   (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (x)))
  64                && xrs*yis != xis*yrs)
  65            ? 0 : xrs * yis);
  66    }
  67    else {
  68       /* x = xr + i*(+/-inf) with |xr| != inf */
  69       u = (   mpfr_zero_p (mpc_imagref (y))
  70            || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_realref (y)))
  71            || (mpfr_inf_p (mpc_realref (y)) && xrs*yrs == xis*yis)
  72            ? 0 : -xis * yis);
  73       v = (   mpfr_zero_p (mpc_realref (y))
  74            || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_imagref (y)))
  75            || (mpfr_inf_p (mpc_imagref (y)) && xrs*yis != xis*yrs)
  76            ? 0 : xis * yrs);
  77    }
  78
  79    if (u == 0 && v == 0) {
  80       /* Naive result is NaN+i*NaN. Obtain an infinity using the algorithm
  81          given in Annex G.5.1 of the ISO C99 standard */
  82       int xr = (mpfr_zero_p (mpc_realref (x)) || mpfr_nan_p (mpc_realref (x)) ? 0
  83                 : (mpfr_inf_p (mpc_realref (x)) ? 1 : 0));
  84       int xi = (mpfr_zero_p (mpc_imagref (x)) || mpfr_nan_p (mpc_imagref (x)) ? 0
  85                 : (mpfr_inf_p (mpc_imagref (x)) ? 1 : 0));
  86       int yr = (mpfr_zero_p (mpc_realref (y)) || mpfr_nan_p (mpc_realref (y)) ? 0 : 1);
  87       int yi = (mpfr_zero_p (mpc_imagref (y)) || mpfr_nan_p (mpc_imagref (y)) ? 0 : 1);
  88       if (mpc_inf_p (y)) {
  89          yr = mpfr_inf_p (mpc_realref (y)) ? 1 : 0;
  90          yi = mpfr_inf_p (mpc_imagref (y)) ? 1 : 0;
  91       }
  92
  93       u = xrs * xr * yrs * yr - xis * xi * yis * yi;
  94       v = xrs * xr * yis * yi + xis * xi * yrs * yr;
  95    }
  96
  97    if (u == 0)
  98       mpfr_set_nan (mpc_realref (z));
  99    else
 100       mpfr_set_inf (mpc_realref (z), u);
 101
 102    if (v == 0)
 103       mpfr_set_nan (mpc_imagref (z));
 104    else
 105       mpfr_set_inf (mpc_imagref (z), v);
 106
 107    return MPC_INEX (0, 0); /* exact */
 108 }
 109
 110
 111 /* compute z = x*y for Im(y) == 0 */
 112 static int
 113 mul_real (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
 114 {
 115    int xrs, xis, yrs, yis;
 116    int inex;
 117
 118    /* save signs of operands */
 119    xrs = MPFR_SIGNBIT (mpc_realref (x));
 120    xis = MPFR_SIGNBIT (mpc_imagref (x));
 121    yrs = MPFR_SIGNBIT (mpc_realref (y));
 122    yis = MPFR_SIGNBIT (mpc_imagref (y));
 123
 124    inex = mpc_mul_fr (z, x, mpc_realref (y), rnd);
 125    /* Signs of zeroes may be wrong. Their correction does not change the
 126       inexact flag. */
 127    if (mpfr_zero_p (mpc_realref (z)))
 128       mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == GMP_RNDD
 129                      || (xrs != yrs && xis == yis), GMP_RNDN);
 130    if (mpfr_zero_p (mpc_imagref (z)))
 131       mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD
 132                      || (xrs != yis && xis != yrs), GMP_RNDN);
 133
 134    return inex;
 135 }
 136
 137
 138 /* compute z = x*y for Re(y) == 0, and Im(x) != 0 and Im(y) != 0 */
 139 static int
 140 mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
 141 {
 142    int sign;
 143    int inex_re, inex_im;
 144    int overlap = z == x || z == y;
 145    mpc_t rop;
 146
 147    if (overlap)
 148       mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z));
 149    else
 150       rop [0] = z[0];
 151
 152    sign =    (MPFR_SIGNBIT (mpc_realref (y)) != MPFR_SIGNBIT (mpc_imagref (x)))
 153           && (MPFR_SIGNBIT (mpc_imagref (y)) != MPFR_SIGNBIT (mpc_realref (x)));
 154
 155    inex_re = -mpfr_mul (mpc_realref (rop), mpc_imagref (x), mpc_imagref (y),
 156                         INV_RND (MPC_RND_RE (rnd)));
 157    mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); /* exact */
 158    inex_im = mpfr_mul (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y),
 159                        MPC_RND_IM (rnd));
 160    mpc_set (z, rop, MPC_RNDNN);
 161
 162    /* Sign of zeroes may be wrong (note that Re(z) cannot be zero) */
 163    if (mpfr_zero_p (mpc_imagref (z)))
 164       mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD
 165                      || sign, GMP_RNDN);
 166
 167    if (overlap)
 168       mpc_clear (rop);
 169
 170    return MPC_INEX (inex_re, inex_im);
 171 }
 172
 173
 174 static int
 175 mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
 176            mpfr_srcptr d, int sign, mpfr_rnd_t rnd)
 177 {
 178    /* Computes z = ab+cd if sign >= 0, or z = ab-cd if sign < 0.
 179       Assumes that a, b, c, d are finite and non-zero; so any multiplication
 180       of two of them yielding an infinity is an overflow, and a
 181       multiplication yielding 0 is an underflow.
 182       Assumes further that z is distinct from a, b, c, d. */
 183
 184    int inex;
 185    mpfr_t u, v;
 186
 187    /* u=a*b, v=sign*c*d exactly */
 188    mpfr_init2 (u, mpfr_get_prec (a) + mpfr_get_prec (b));
 189    mpfr_init2 (v, mpfr_get_prec (c) + mpfr_get_prec (d));
 190    mpfr_mul (u, a, b, GMP_RNDN);
 191    mpfr_mul (v, c, d, GMP_RNDN);
 192    if (sign < 0)
 193       mpfr_neg (v, v, GMP_RNDN);
 194
 195    /* tentatively compute z as u+v; here we need z to be distinct
 196       from a, b, c, d to not lose the latter */
 197    inex = mpfr_add (z, u, v, rnd);
 198
 199    if (mpfr_inf_p (z)) {
 200       /* replace by "correctly rounded overflow" */
 201       mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN);
 202       inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd);
 203    }
 204    else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) {
 205       /* exactly u underflowed, determine inexact flag */
 206       inex = (mpfr_signbit (u) ? 1 : -1);
 207    }
 208    else if (mpfr_zero_p (v) && !mpfr_zero_p (u)) {
 209       /* exactly v underflowed, determine inexact flag */
 210       inex = (mpfr_signbit (v) ? 1 : -1);
 211    }
 212    else if (mpfr_nan_p (z) || (mpfr_zero_p (u) && mpfr_zero_p (v))) {
 213       /* In the first case, u and v are infinities with opposite signs.
 214          In the second case, u and v are zeroes; their sum may be 0 or the
 215          least representable number, with a sign to be determined.
 216          Redo the computations with mpz_t exponents */
 217       mpfr_exp_t ea, eb, ec, ed;
 218       mpz_t eu, ev;
 219          /* cheat to work around the const qualifiers */
 220
 221       /* Normalise the input by shifting and keep track of the shifts in
 222          the exponents of u and v */
 223       ea = mpfr_get_exp (a);
 224       eb = mpfr_get_exp (b);
 225       ec = mpfr_get_exp (c);
 226       ed = mpfr_get_exp (d);
 227
 228       mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0);
 229       mpfr_set_exp ((mpfr_ptr) b, (mpfr_prec_t) 0);
 230       mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0);
 231       mpfr_set_exp ((mpfr_ptr) d, (mpfr_prec_t) 0);
 232
 233       mpz_init (eu);
 234       mpz_init (ev);
 235       mpz_set_si (eu, (long int) ea);
 236       mpz_add_si (eu, eu, (long int) eb);
 237       mpz_set_si (ev, (long int) ec);
 238       mpz_add_si (ev, ev, (long int) ed);
 239
 240       /* recompute u and v and move exponents to eu and ev */
 241       mpfr_mul (u, a, b, GMP_RNDN);
 242       /* exponent of u is non-positive */
 243       mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u)));
 244       mpfr_set_exp (u, (mpfr_prec_t) 0);
 245       mpfr_mul (v, c, d, GMP_RNDN);
 246       if (sign < 0)
 247          mpfr_neg (v, v, GMP_RNDN);
 248       mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v)));
 249       mpfr_set_exp (v, (mpfr_prec_t) 0);
 250
 251       if (mpfr_nan_p (z)) {
 252          mpfr_exp_t emax = mpfr_get_emax ();
 253          int overflow;
 254          /* We have a = ma * 2^ea with 1/2 <= |ma| < 1 and ea <= emax, and
 255             analogously for b. So eu <= 2*emax, and eu > emax since we have
 256             an overflow. The same holds for ev. Shift u and v by as much as
 257             possible so that one of them has exponent emax and the
 258             remaining exponents in eu and ev are the same. Then carry out
 259             the addition. Shifting u and v prevents an underflow. */
 260          if (mpz_cmp (eu, ev) >= 0) {
 261             mpfr_set_exp (u, emax);
 262             mpz_sub_ui (eu, eu, (long int) emax);
 263             mpz_sub (ev, ev, eu);
 264             mpfr_set_exp (v, (mpfr_exp_t) mpz_get_ui (ev));
 265                /* remaining common exponent is now in eu */
 266          }
 267          else {
 268             mpfr_set_exp (v, emax);
 269             mpz_sub_ui (ev, ev, (long int) emax);
 270             mpz_sub (eu, eu, ev);
 271             mpfr_set_exp (u, (mpfr_exp_t) mpz_get_ui (eu));
 272             mpz_set (eu, ev);
 273                /* remaining common exponent is now also in eu */
 274          }
 275          inex = mpfr_add (z, u, v, rnd);
 276             /* Result is finite since u and v have different signs. */
 277          overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd);
 278          if (overflow)
 279             inex = overflow;
 280       }
 281       else {
 282          int underflow;
 283          /* Addition of two zeroes with same sign. We have a = ma * 2^ea
 284             with 1/2 <= |ma| < 1 and ea >= emin and similarly for b.
 285             So 2*emin < 2*emin+1 <= eu < emin < 0, and analogously for v. */
 286          mpfr_exp_t emin = mpfr_get_emin ();
 287          if (mpz_cmp (eu, ev) <= 0) {
 288             mpfr_set_exp (u, emin);
 289             mpz_add_ui (eu, eu, (unsigned long int) (-emin));
 290             mpz_sub (ev, ev, eu);
 291             mpfr_set_exp (v, (mpfr_exp_t) mpz_get_si (ev));
 292          }
 293          else {
 294             mpfr_set_exp (v, emin);
 295             mpz_add_ui (ev, ev, (unsigned long int) (-emin));
 296             mpz_sub (eu, eu, ev);
 297             mpfr_set_exp (u, (mpfr_exp_t) mpz_get_si (eu));
 298             mpz_set (eu, ev);
 299          }
 300          inex = mpfr_add (z, u, v, rnd);
 301          mpz_neg (eu, eu);
 302          underflow = mpfr_div_2ui (z, z, mpz_get_ui (eu), rnd);
 303          if (underflow)
 304             inex = underflow;
 305       }
 306
 307       mpz_clear (eu);
 308       mpz_clear (ev);
 309
 310       mpfr_set_exp ((mpfr_ptr) a, ea);
 311       mpfr_set_exp ((mpfr_ptr) b, eb);
 312       mpfr_set_exp ((mpfr_ptr) c, ec);
 313       mpfr_set_exp ((mpfr_ptr) d, ed);
 314          /* works also when some of a, b, c, d are not all distinct */
 315    }
 316
 317    mpfr_clear (u);
 318    mpfr_clear (v);
 319
 320    return inex;
 321 }
 322
 323
 324 int
 325 mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
 326 {
 327    /* computes z=x*y by the schoolbook method, where x and y are assumed
 328       to be finite and without zero parts                                */
 329    int overlap, inex;
 330    mpc_t rop;
 331
 332    MPC_ASSERT (   mpfr_regular_p (mpc_realref (x)) && mpfr_regular_p (mpc_imagref (x))
 333                && mpfr_regular_p (mpc_realref (y)) && mpfr_regular_p (mpc_imagref (y)));
 334    overlap = (z == x) || (z == y);
 335    if (overlap)
 336       mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z));
 337    else
 338       rop [0] = z [0];
 339
 340    inex = MPC_INEX (mpfr_fmma (mpc_realref (rop), mpc_realref (x), mpc_realref (y), mpc_imagref (x),
 341                                mpc_imagref (y), -1, MPC_RND_RE (rnd)),
 342                     mpfr_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), mpc_imagref (x),
 343                                mpc_realref (y), +1, MPC_RND_IM (rnd)));
 344
 345    mpc_set (z, rop, MPC_RNDNN);
 346    if (overlap)
 347       mpc_clear (rop);
 348
 349    return inex;
 350 }
 351
 352
 353 int
 354 mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
 355 {
 356    /* computes rop=op1*op2 by a Karatsuba algorithm, where op1 and op2
 357       are assumed to be finite and without zero parts                  */
 358   mpfr_srcptr a, b, c, d;
 359   int mul_i, ok, inexact, mul_a, mul_c, inex_re = 0, inex_im = 0, sign_x, sign_u;
 360   mpfr_t u, v, w, x;
 361   mpfr_prec_t prec, prec_re, prec_u, prec_v, prec_w;
 362   mpfr_rnd_t rnd_re, rnd_u;
 363   int overlap;
 364      /* true if rop == op1 or rop == op2 */
 365   mpc_t result;
 366      /* overlap is quite difficult to handle, because we have to tentatively
 367         round the variable u in the end to either the real or the imaginary
 368         part of rop (it is not possible to tell now whether the real or
 369         imaginary part is used). If this fails, we have to start again and
 370         need the correct values of op1 and op2.
 371         So we just create a new variable for the result in this case. */
 372   int loop;
 373   const int MAX_MUL_LOOP = 1;
 374
 375   overlap = (rop == op1) || (rop == op2);
 376   if (overlap)
 377      mpc_init3 (result, MPC_PREC_RE (rop), MPC_PREC_IM (rop));
 378   else
 379      result [0] = rop [0];
 380
 381   a = mpc_realref(op1);
 382   b = mpc_imagref(op1);
 383   c = mpc_realref(op2);
 384   d = mpc_imagref(op2);
 385
 386   /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */
 387
 388   mul_i = 0; /* number of multiplications by i */
 389   mul_a = 1; /* implicit factor for a */
 390   mul_c = 1; /* implicit factor for c */
 391
 392   if (mpfr_cmp_abs (a, b) < 0)
 393     {
 394       MPFR_SWAP (a, b);
 395       mul_i ++;
 396       mul_a = -1; /* consider i * (a+i*b) = -b + i*a */
 397     }
 398
 399   if (mpfr_cmp_abs (c, d) < 0)
 400     {
 401       MPFR_SWAP (c, d);
 402       mul_i ++;
 403       mul_c = -1; /* consider -d + i*c instead of c + i*d */
 404     }
 405
 406   /* find the precision and rounding mode for the new real part */
 407   if (mul_i % 2)
 408     {
 409       prec_re = MPC_PREC_IM(rop);
 410       rnd_re = MPC_RND_IM(rnd);
 411     }
 412   else /* mul_i = 0 or 2 */
 413     {
 414       prec_re = MPC_PREC_RE(rop);
 415       rnd_re = MPC_RND_RE(rnd);
 416     }
 417
 418   if (mul_i)
 419     rnd_re = INV_RND(rnd_re);
 420
 421   /* now |a| >= |b| and |c| >= |d| */
 422   prec = MPC_MAX_PREC(rop);
 423
 424   mpfr_init2 (v, prec_v = mpfr_get_prec (a) + mpfr_get_prec (d));
 425   mpfr_init2 (w, prec_w = mpfr_get_prec (b) + mpfr_get_prec (c));
 426   mpfr_init2 (u, 2);
 427   mpfr_init2 (x, 2);
 428
 429   inexact = mpfr_mul (v, a, d, GMP_RNDN);
 430   if (inexact) {
 431      /* over- or underflow */
 432     ok = 0;
 433     goto clear;
 434   }
 435   if (mul_a == -1)
 436     mpfr_neg (v, v, GMP_RNDN);
 437
 438   inexact = mpfr_mul (w, b, c, GMP_RNDN);
 439   if (inexact) {
 440      /* over- or underflow */
 441      ok = 0;
 442      goto clear;
 443   }
 444   if (mul_c == -1)
 445     mpfr_neg (w, w, GMP_RNDN);
 446
 447   /* compute sign(v-w) */
 448   sign_x = mpfr_cmp_abs (v, w);
 449   if (sign_x > 0)
 450     sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w);
 451   else if (sign_x == 0)
 452     sign_x = mpfr_sgn (v) - mpfr_sgn (w);
 453   else
 454     sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w);
 455
 456    sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c);
 457
 458    if (sign_x * sign_u < 0)
 459     {  /* swap inputs */
 460       MPFR_SWAP (a, c);
 461       MPFR_SWAP (b, d);
 462       mpfr_swap (v, w);
 463       { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; }
 464       sign_x = - sign_x;
 465     }
 466
 467    /* now sign_x * sign_u >= 0 */
 468    loop = 0;
 469    do
 470      {
 471         loop++;
 472          /* the following should give failures with prob. <= 1/prec */
 473          prec += mpc_ceil_log2 (prec) + 3;
 474
 475          mpfr_set_prec (u, prec_u = prec);
 476          mpfr_set_prec (x, prec);
 477
 478          /* first compute away(b +/- a) and store it in u */
 479          inexact = (mul_a == -1 ?
 480                     ROUND_AWAY (mpfr_sub (u, b, a, MPFR_RNDA), u) :
 481                     ROUND_AWAY (mpfr_add (u, b, a, MPFR_RNDA), u));
 482
 483          /* then compute away(+/-c - d) and store it in x */
 484          inexact |= (mul_c == -1 ?
 485                      ROUND_AWAY (mpfr_add (x, c, d, MPFR_RNDA), x) :
 486                      ROUND_AWAY (mpfr_sub (x, c, d, MPFR_RNDA), x));
 487          if (mul_c == -1)
 488            mpfr_neg (x, x, GMP_RNDN);
 489
 490          if (inexact == 0)
 491             mpfr_prec_round (u, prec_u = 2 * prec, GMP_RNDN);
 492
 493          /* compute away(u*x) and store it in u */
 494          inexact |= ROUND_AWAY (mpfr_mul (u, u, x, MPFR_RNDA), u);
 495             /* (a+b)*(c-d) */
 496
 497          /* if all computations are exact up to here, it may be that
 498             the real part is exact, thus we need if possible to
 499             compute v - w exactly */
 500          if (inexact == 0)
 501            {
 502              mpfr_prec_t prec_x;
 503              /* v and w are different from 0, so mpfr_get_exp is safe to use */
 504              prec_x = SAFE_ABS (mpfr_exp_t, mpfr_get_exp (v) - mpfr_get_exp (w))
 505                       + MPC_MAX (prec_v, prec_w) + 1;
 506                       /* +1 is necessary for a potential carry */
 507              /* ensure we do not use a too large precision */
 508              if (prec_x > prec_u)
 509                prec_x = prec_u;
 510              if (prec_x > prec)
 511                mpfr_prec_round (x, prec_x, GMP_RNDN);
 512            }
 513
 514          rnd_u = (sign_u > 0) ? GMP_RNDU : GMP_RNDD;
 515          inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */
 516
 517          /* in case u=0, ensure that rnd_u rounds x away from zero */
 518          if (mpfr_sgn (u) == 0)
 519            rnd_u = (mpfr_sgn (x) > 0) ? GMP_RNDU : GMP_RNDD;
 520          inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */
 521
 522          ok = inexact == 0 ||
 523            mpfr_can_round (u, prec_u - 3, rnd_u, GMP_RNDZ,
 524                            prec_re + (rnd_re == GMP_RNDN));
 525          /* this ensures both we can round correctly and determine the correct
 526             inexact flag (for rounding to nearest) */
 527      }
 528    while (!ok && loop <= MAX_MUL_LOOP);
 529       /* after MAX_MUL_LOOP rounds, use mpc_naive instead */
 530
 531    if (ok) {
 532       /* if inexact is zero, then u is exactly ac-bd, otherwise fix the sign
 533          of the inexact flag for u, which was rounded away from ac-bd */
 534       if (inexact != 0)
 535       inexact = mpfr_sgn (u);
 536
 537       if (mul_i == 0)
 538       {
 539          inex_re = mpfr_set (mpc_realref(result), u, MPC_RND_RE(rnd));
 540          if (inex_re == 0)
 541             {
 542             inex_re = inexact; /* u is rounded away from 0 */
 543             inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd));
 544             }
 545          else
 546             inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd));
 547       }
 548       else if (mul_i == 1) /* (x+i*y)/i = y - i*x */
 549       {
 550          inex_im = mpfr_neg (mpc_imagref(result), u, MPC_RND_IM(rnd));
 551          if (inex_im == 0)
 552             {
 553             inex_im = -inexact; /* u is rounded away from 0 */
 554             inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd));
 555             }
 556          else
 557             inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd));
 558       }
 559       else /* mul_i = 2, z/i^2 = -z */
 560       {
 561          inex_re = mpfr_neg (mpc_realref(result), u, MPC_RND_RE(rnd));
 562          if (inex_re == 0)
 563             {
 564             inex_re = -inexact; /* u is rounded away from 0 */
 565             inex_im = -mpfr_add (mpc_imagref(result), v, w,
 566                                  INV_RND(MPC_RND_IM(rnd)));
 567             mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd));
 568             }
 569          else
 570             {
 571             inex_im = -mpfr_add (mpc_imagref(result), v, w,
 572                                  INV_RND(MPC_RND_IM(rnd)));
 573             mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd));
 574             }
 575       }
 576
 577       mpc_set (rop, result, MPC_RNDNN);
 578    }
 579
 580 clear:
 581    mpfr_clear (u);
 582    mpfr_clear (v);
 583    mpfr_clear (w);
 584    mpfr_clear (x);
 585    if (overlap)
 586       mpc_clear (result);
 587
 588    if (ok)
 589       return MPC_INEX(inex_re, inex_im);
 590    else
 591       return mpc_mul_naive (rop, op1, op2, rnd);
 592 }
 593
 594
 595 int
 596 mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
 597 {
 598    /* Conforming to ISO C99 standard (G.5.1 multiplicative operators),
 599       infinities are treated specially if both parts are NaN when computed
 600       naively. */
 601    if (mpc_inf_p (b))
 602       return mul_infinite (a, b, c);
 603    if (mpc_inf_p (c))
 604       return mul_infinite (a, c, b);
 605
 606    /* NaN contamination of both parts in result */
 607    if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b))
 608        || mpfr_nan_p (mpc_realref (c)) || mpfr_nan_p (mpc_imagref (c))) {
 609       mpfr_set_nan (mpc_realref (a));
 610       mpfr_set_nan (mpc_imagref (a));
 611       return MPC_INEX (0, 0);
 612    }
 613
 614    /* check for real multiplication */
 615    if (mpfr_zero_p (mpc_imagref (b)))
 616       return mul_real (a, c, b, rnd);
 617    if (mpfr_zero_p (mpc_imagref (c)))
 618       return mul_real (a, b, c, rnd);
 619
 620    /* check for purely imaginary multiplication */
 621    if (mpfr_zero_p (mpc_realref (b)))
 622       return mul_imag (a, c, b, rnd);
 623    if (mpfr_zero_p (mpc_realref (c)))
 624       return mul_imag (a, b, c, rnd);
 625
 626    /* If the real and imaginary part of one argument have a very different */
 627    /* exponent, it is not reasonable to use Karatsuba multiplication.      */
 628    if (   SAFE_ABS (mpfr_exp_t,
 629                      mpfr_get_exp (mpc_realref (b)) - mpfr_get_exp (mpc_imagref (b)))
 630          > (mpfr_exp_t) MPC_MAX_PREC (b) / 2
 631       || SAFE_ABS (mpfr_exp_t,
 632                      mpfr_get_exp (mpc_realref (c)) - mpfr_get_exp (mpc_imagref (c)))
 633          > (mpfr_exp_t) MPC_MAX_PREC (c) / 2)
 634       return mpc_mul_naive (a, b, c, rnd);
 635    else
 636       return ((MPC_MAX_PREC(a)
 637                <= (mpfr_prec_t) MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB)
 638             ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd);
 639 }