src/mul.c

   1 /* mpc_mul -- Multiply two complex numbers
   2
   3 Copyright (C) 2002, 2004, 2005, 2008, 2009 Andreas Enge, Paul Zimmermann, Philippe Th\'eveny
   4
   5 This file is part of the MPC Library.
   6
   7 The MPC Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 2.1 of the License, or (at your
  10 option) any later version.
  11
  12 The MPC Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the MPC Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include "mpc-impl.h"
  23
  24 static int mul_infinite (mpc_ptr a, mpc_srcptr u, mpc_srcptr v);
  25 static int mul_pure_imaginary (mpc_ptr a, mpc_srcptr u, mpfr_srcptr y,
  26                                mpc_rnd_t rnd, int overlap);
  27
  28 /* return inex such that MPC_INEX_RE(inex) = -1, 0, 1
  29                          MPC_INEX_IM(inex) = -1, 0, 1
  30    depending on the signs of the real/imaginary parts of the result
  31 */
  32 int
  33 mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
  34 {
  35   int brs, bis, crs, cis;
  36
  37   /* Conforming to ISO C99 standard (G.5.1 multiplicative operators),
  38      infinities have special treatment if both parts are NaN when computed
  39      naively. */
  40   if (mpc_inf_p (b))
  41     return mul_infinite (a, b, c);
  42   if (mpc_inf_p (c))
  43     return mul_infinite (a, c, b);
  44
  45   /* NaN contamination of both part in result */
  46   if (mpfr_nan_p (MPC_RE (b)) || mpfr_nan_p (MPC_IM (b))
  47       || mpfr_nan_p (MPC_RE (c)) || mpfr_nan_p (MPC_IM (c)))
  48     {
  49       mpfr_set_nan (MPC_RE (a));
  50       mpfr_set_nan (MPC_IM (a));
  51       return MPC_INEX (0, 0);
  52     }
  53
  54   /* save operands' sign */
  55   brs = MPFR_SIGNBIT (MPC_RE (b));
  56   bis = MPFR_SIGNBIT (MPC_IM (b));
  57   crs = MPFR_SIGNBIT (MPC_RE (c));
  58   cis = MPFR_SIGNBIT (MPC_IM (c));
  59
  60   /* first check for real multiplication */
  61   if (mpfr_zero_p (MPC_IM (b))) /* b * (x+i*y) = b*x + i *(b*y) */
  62     {
  63       int inex;
  64       inex = mpc_mul_fr (a, c, MPC_RE (b), rnd);
  65
  66       /* Sign of zeros is wrong in some cases. This correction doesn't change
  67          inexact flag. */
  68       if (mpfr_zero_p (MPC_RE (a)))
  69         mpfr_setsign (MPC_RE (a), MPC_RE (a), MPC_RND_RE(rnd) == GMP_RNDD
  70                       || (brs != crs && bis == cis), GMP_RNDN); /* exact */
  71       if (mpfr_zero_p (MPC_IM (a)))
  72         mpfr_setsign (MPC_IM (a), MPC_IM (a), MPC_RND_IM (rnd) == GMP_RNDD
  73                       || (brs != cis && bis != crs), GMP_RNDN); /* exact */
  74
  75       return inex;
  76     }
  77   if (mpfr_zero_p (MPC_IM (c)))
  78     {
  79       int inex;
  80       inex = mpc_mul_fr (a, b, MPC_RE (c), rnd);
  81
  82       /* Sign of zeros is wrong in some cases. This correction doesn't change
  83          inexact flag. */
  84       if (mpfr_zero_p (MPC_RE (a)))
  85         mpfr_setsign (MPC_RE (a), MPC_RE (a), MPC_RND_RE(rnd) == GMP_RNDD
  86                       || (brs != crs && bis == cis), GMP_RNDN);
  87       if (mpfr_zero_p (MPC_IM (a)))
  88         mpfr_setsign (MPC_IM (a), MPC_IM (a), MPC_RND_IM (rnd) == GMP_RNDD
  89                       || (brs != cis && bis != crs), GMP_RNDN);
  90
  91       return inex;
  92     }
  93
  94   /* check for purely complex multiplication */
  95   if (mpfr_zero_p (MPC_RE (b))) /* i*b * (x+i*y) = -b*y + i*(b*x) */
  96     {
  97       int inex;
  98       inex = mul_pure_imaginary (a, c, MPC_IM (b), rnd, (a == b || a == c));
  99
 100       /* Sign of zeros is wrong in some cases (note that Re(a) cannot be a
 101          zero) */
 102       if (mpfr_zero_p (MPC_IM (a)))
 103         mpfr_setsign (MPC_IM (a), MPC_IM (a), MPC_RND_IM (rnd) == GMP_RNDD
 104                       || (brs != cis && bis != crs), GMP_RNDN);
 105
 106       return inex;
 107     }
 108   if (mpfr_zero_p (MPC_RE (c)))
 109     /* note that a cannot be a zero */
 110     return mul_pure_imaginary (a, b, MPC_IM (c), rnd, (a == b || a == c));
 111
 112   /* If the real and imaginary part of one argument have a very different */
 113   /* exponent, it is not reasonable to use Karatsuba multiplication.      */
 114   if (   SAFE_ABS (mp_exp_t, MPFR_EXP (MPC_RE (b)) - MPFR_EXP (MPC_IM (b)))
 115          > (mp_exp_t) MPC_MAX_PREC (b) / 2
 116          || SAFE_ABS (mp_exp_t, MPFR_EXP (MPC_RE (c)) - MPFR_EXP (MPC_IM (c)))
 117          > (mp_exp_t) MPC_MAX_PREC (c) / 2)
 118     return mpc_mul_naive (a, b, c, rnd);
 119   else
 120     return ((MPC_MAX_PREC(a)
 121              <= (mp_prec_t) MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB)
 122             ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd);
 123 }
 124
 125 /* compute a=u*v when u has an infinite part */
 126 static int
 127 mul_infinite (mpc_ptr a, mpc_srcptr u, mpc_srcptr v)
 128 {
 129   /* let u = ur+i*ui and v =vr+i*vi */
 130   int x, y;
 131
 132   /* save operands' sign */
 133   int urs = mpfr_signbit (MPC_RE (u)) ? -1 : 1;
 134   int uis = mpfr_signbit (MPC_IM (u)) ? -1 : 1;
 135   int vrs = mpfr_signbit (MPC_RE (v)) ? -1 : 1;
 136   int vis = mpfr_signbit (MPC_IM (v)) ? -1 : 1;
 137
 138   /* compute the sign of
 139      x = urs * ur * vrs * vr - uis * ui * vis * vi
 140      y = urs * ur * vis * vi + uis * ui * vrs * vr
 141      +1 if positive, -1 if negative, 0 if zero or NaN */
 142   if (mpfr_nan_p (MPC_RE (u)) || mpfr_nan_p (MPC_IM (u))
 143       || mpfr_nan_p (MPC_RE (v)) || mpfr_nan_p (MPC_IM (v)))
 144     {
 145       x = 0;
 146       y = 0;
 147     }
 148   else if (mpfr_inf_p (MPC_RE (u)))
 149     {
 150       /* u = (+/-inf) +i*v */
 151       x = mpfr_zero_p (MPC_RE (v))
 152         || (mpfr_inf_p (MPC_IM (u)) && mpfr_zero_p (MPC_IM (v)))
 153         || (mpfr_zero_p (MPC_IM (u)) && mpfr_inf_p (MPC_IM (v)))
 154         || ((mpfr_inf_p (MPC_IM (u)) || mpfr_inf_p (MPC_IM (v)))
 155             && urs*vrs == uis*vis) ?
 156         0 : urs * vrs;
 157       y = mpfr_zero_p (MPC_IM (v))
 158         || (mpfr_inf_p (MPC_IM (u)) && mpfr_zero_p (MPC_RE (v)))
 159         || (mpfr_zero_p (MPC_IM (u)) && mpfr_inf_p (MPC_RE (v)))
 160         || ((mpfr_inf_p (MPC_IM (u)) || mpfr_inf_p (MPC_IM (u)))
 161             && urs*vis == -uis*vrs) ?
 162         0 : urs * vis;
 163     }
 164   else
 165     {
 166       /* u = u +i*(+/-inf) where |u| < inf */
 167       x = mpfr_zero_p (MPC_IM (v))
 168         || (mpfr_zero_p (MPC_RE (u)) && mpfr_inf_p (MPC_RE (v)))
 169         || (mpfr_inf_p (MPC_RE (v)) && urs*vrs == uis*vis) ?
 170         0 : -uis * vis;
 171       y = mpfr_zero_p (MPC_RE (v))
 172         || (mpfr_zero_p (MPC_RE (u)) && mpfr_inf_p (MPC_IM (v)))
 173         || (mpfr_inf_p (MPC_IM (v)) && urs*vis == -uis*vrs) ?
 174         0 : uis * vrs;
 175     }
 176
 177   /* Naive result is NaN+i*NaN. We will try to obtain an infinite using
 178      algorithm given in Annex G of ISO C99 standard */
 179   if (x == 0 && y ==0)
 180     {
 181       int ur = mpfr_zero_p (MPC_RE (u)) || mpfr_nan_p (MPC_RE (u)) ?
 182         0 : 1;
 183       int ui = mpfr_zero_p (MPC_IM (u)) || mpfr_nan_p (MPC_IM (u)) ?
 184         0 : 1;
 185       int vr = mpfr_zero_p (MPC_RE (v)) || mpfr_nan_p (MPC_RE (v)) ?
 186         0 : 1;
 187       int vi = mpfr_zero_p (MPC_IM (v)) || mpfr_nan_p (MPC_IM (v)) ?
 188         0 : 1;
 189       if (mpc_inf_p (u))
 190         {
 191           ur = mpfr_inf_p (MPC_RE (u)) ? 1 : 0;
 192           ui = mpfr_inf_p (MPC_IM (u)) ? 1 : 0;
 193         }
 194       if (mpc_inf_p (v))
 195         {
 196           vr = mpfr_inf_p (MPC_RE (v)) ? 1 : 0;
 197           vi = mpfr_inf_p (MPC_IM (v)) ? 1 : 0;
 198         }
 199
 200       x = urs * ur * vrs * vr - uis * ui * vis * vi;
 201       y = urs * ur * vis * vi + uis * ui * vrs * vr;
 202     }
 203
 204   if (x == 0)
 205     mpfr_set_nan (MPC_RE (a));
 206   else
 207     mpfr_set_inf (MPC_RE (a), x);
 208
 209   if (y == 0)
 210     mpfr_set_nan (MPC_IM (a));
 211   else
 212     mpfr_set_inf (MPC_IM (a), y);
 213
 214   return MPC_INEX (0, 0); /* exact */
 215 }
 216
 217 static int
 218 mul_pure_imaginary (mpc_ptr a, mpc_srcptr u, mpfr_srcptr y, mpc_rnd_t rnd,
 219                     int overlap)
 220 {
 221   int inex_re, inex_im;
 222   mpc_t result;
 223
 224   if (overlap)
 225     mpc_init3 (result, MPFR_PREC (MPC_RE (a)), MPFR_PREC (MPC_IM (a)));
 226   else
 227     result [0] = a [0];
 228
 229   inex_re = -mpfr_mul (MPC_RE(result), MPC_IM(u), y, INV_RND(MPC_RND_RE(rnd)));
 230   mpfr_neg (MPC_RE (result), MPC_RE (result), GMP_RNDN);
 231   inex_im = mpfr_mul (MPC_IM(result), MPC_RE(u), y, MPC_RND_IM(rnd));
 232   mpc_set (a, result, MPC_RNDNN);
 233
 234   if (overlap)
 235     mpc_clear (result);
 236
 237   return MPC_INEX(inex_re, inex_im);
 238 }
 239
 240 int
 241 mpc_mul_naive (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
 242 {
 243   int overlap, inex_re, inex_im;
 244   mpfr_t u, v, t;
 245   mp_prec_t prec;
 246
 247   overlap = (a == b) || (a == c);
 248
 249   prec = MPC_MAX_PREC(b) + MPC_MAX_PREC(c);
 250
 251   mpfr_init2 (u, prec);
 252   mpfr_init2 (v, prec);
 253
 254   /* Re(a) = Re(b)*Re(c) - Im(b)*Im(c) */
 255   /* FIXME: this code suffers undue overflows: u or v can overflow while u-v
 256      or u+v is representable */
 257   mpfr_mul (u, MPC_RE(b), MPC_RE(c), GMP_RNDN); /* exact */
 258   mpfr_mul (v, MPC_IM(b), MPC_IM(c), GMP_RNDN); /* exact */
 259
 260   if (overlap)
 261     {
 262       mpfr_init2 (t, MPFR_PREC(MPC_RE(a)));
 263       inex_re = mpfr_sub (t, u, v, MPC_RND_RE(rnd));
 264     }
 265   else
 266     inex_re = mpfr_sub (MPC_RE(a), u, v, MPC_RND_RE(rnd));
 267
 268   /* Im(a) = Re(b)*Im(c) + Im(b)*Re(c) */
 269   mpfr_mul (u, MPC_RE(b), MPC_IM(c), GMP_RNDN); /* exact */
 270   if (b == c) /* square case */
 271     inex_im = mpfr_mul_2exp (MPC_IM(a), u, 1, MPC_RND_IM(rnd));
 272   else
 273     {
 274       mpfr_mul (v, MPC_IM(b), MPC_RE(c), GMP_RNDN); /* exact */
 275       inex_im = mpfr_add (MPC_IM(a), u, v, MPC_RND_IM(rnd));
 276     }
 277
 278   mpfr_clear (u);
 279   mpfr_clear (v);
 280
 281   if (overlap)
 282     {
 283       mpfr_set (MPC_RE(a), t, GMP_RNDN); /* exact */
 284       mpfr_clear (t);
 285     }
 286
 287   return MPC_INEX(inex_re, inex_im);
 288 }
 289
 290
 291 /* Karatsuba multiplication, with 3 real multiplies */
 292 int
 293 mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
 294 {
 295   mpfr_srcptr a, b, c, d;
 296   int mul_i, ok, inexact, mul_a, mul_c, inex_re, inex_im, sign_x, sign_u;
 297   mpfr_t u, v, w, x;
 298   mp_prec_t prec, prec_re, prec_u, prec_v, prec_w;
 299   mp_rnd_t rnd_re, rnd_u, rnd_x;
 300   int overlap;
 301      /* true if rop == op1 or rop == op2 */
 302   mpc_t result;
 303      /* overlap is quite difficult to handle, because we have to tentatively
 304         round the variable u in the end to either the real or the imaginary
 305         part of rop (it is not possible to tell now whether the real or
 306         imaginary part is used). If this fails, we have to start again and
 307         need the correct values of op1 and op2.
 308         So we just create a new variable for the result in this case. */
 309
 310   overlap = (rop == op1) || (rop == op2);
 311   if (overlap)
 312      mpc_init3 (result, MPFR_PREC (MPC_RE (rop)),
 313                         MPFR_PREC (MPC_IM (rop)));
 314   else
 315      result [0] = rop [0];
 316
 317   a = MPC_RE(op1);
 318   b = MPC_IM(op1);
 319   c = MPC_RE(op2);
 320   d = MPC_IM(op2);
 321
 322   /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */
 323
 324   mul_i = 0; /* number of multiplications by i */
 325   mul_a = 1; /* implicit factor for a */
 326   mul_c = 1; /* implicit factor for c */
 327
 328   if (mpfr_cmp_abs (a, b) < 0)
 329     {
 330       MPFR_SWAP (a, b);
 331       mul_i ++;
 332       mul_a = -1; /* consider i * (a+i*b) = -b + i*a */
 333     }
 334
 335   if (mpfr_cmp_abs (c, d) < 0)
 336     {
 337       MPFR_SWAP (c, d);
 338       mul_i ++;
 339       mul_c = -1; /* consider -d + i*c instead of c + i*d */
 340     }
 341
 342   /* find the precision and rounding mode for the new real part */
 343   if (mul_i % 2)
 344     {
 345       prec_re = MPFR_PREC(MPC_IM(rop));
 346       rnd_re = MPC_RND_IM(rnd);
 347     }
 348   else /* mul_i = 0 or 2 */
 349     {
 350       prec_re = MPFR_PREC(MPC_RE(rop));
 351       rnd_re = MPC_RND_RE(rnd);
 352     }
 353
 354   if (mul_i)
 355     rnd_re = INV_RND(rnd_re);
 356
 357   /* now |a| >= |b| and |c| >= |d| */
 358   prec = MPC_MAX_PREC(rop);
 359
 360   mpfr_init2 (u, 2);
 361   mpfr_init2 (v, prec_v = MPFR_PREC(a) + MPFR_PREC(d));
 362   mpfr_init2 (w, prec_w = MPFR_PREC(b) + MPFR_PREC(c));
 363   mpfr_init2 (x, 2);
 364
 365   mpfr_mul (v, a, d, GMP_RNDN); /* exact */
 366   if (mul_a == -1)
 367     mpfr_neg (v, v, GMP_RNDN);
 368
 369   mpfr_mul (w, b, c, GMP_RNDN); /* exact */
 370   if (mul_c == -1)
 371     mpfr_neg (w, w, GMP_RNDN);
 372
 373   /* compute sign(v-w) */
 374   sign_x = mpfr_cmp_abs (v, w);
 375   if (sign_x > 0)
 376     sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w);
 377   else if (sign_x == 0)
 378     sign_x = mpfr_sgn (v) - mpfr_sgn (w);
 379   else
 380     sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w);
 381
 382    sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c);
 383
 384    if (sign_x * sign_u < 0)
 385     {  /* swap inputs */
 386       MPFR_SWAP (a, c);
 387       MPFR_SWAP (b, d);
 388       mpfr_swap (v, w);
 389       { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; }
 390       sign_x = - sign_x;
 391     }
 392
 393    /* now sign_x * sign_u >= 0 */
 394    do
 395      {
 396          /* the following should give failures with prob. <= 1/prec */
 397          prec += mpc_ceil_log2 (prec) + 3;
 398
 399          mpfr_set_prec (u, prec_u = prec);
 400          mpfr_set_prec (x, prec);
 401
 402          /* first compute away(b +/- a) and store it in u */
 403          rnd_u = (mpfr_sgn (a) > 0) ? GMP_RNDU : GMP_RNDD;
 404          if (mul_a == -1)
 405            rnd_u = INV_RND(rnd_u);
 406          inexact = ((mul_a == -1) ? mpfr_sub : mpfr_add) (u, b, a, rnd_u);
 407
 408          /* then compute away(+/-c - d) and store it in x */
 409          rnd_x = (mpfr_sgn (c) > 0) ? GMP_RNDU : GMP_RNDD;
 410          inexact |= ((mul_c == -1) ? mpfr_add : mpfr_sub) (x, c, d, rnd_x);
 411          if (mul_c == -1)
 412            mpfr_neg (x, x, GMP_RNDN);
 413
 414          if (inexact == 0)
 415            mpfr_prec_round (u, prec_u = 2 * prec, GMP_RNDN);
 416
 417          /* compute away(u*x) and store it in u */
 418          rnd_u = (sign_u > 0) ? GMP_RNDU : GMP_RNDD;
 419          inexact |= mpfr_mul (u, u, x, rnd_u); /* (a+b)*(c-d) */
 420
 421          /* if all computations are exact up to here, it may be that
 422             the real part is exact, thus we need if possible to
 423             compute v - w exactly */
 424          if (inexact == 0)
 425            {
 426              mp_prec_t prec_x;
 427              if (mpfr_zero_p(v))
 428                prec_x = prec_w;
 429              else if (mpfr_zero_p(w))
 430                prec_x = prec_v;
 431              else
 432                {
 433                  prec_x = (MPFR_EXP(v) > MPFR_EXP(w)) ? MPFR_EXP(v) - MPFR_EXP(w)
 434                    : MPFR_EXP(w) - MPFR_EXP(v);
 435                  prec_x += MPC_MAX (prec_v, prec_w) + 1;
 436                }
 437            /* +1 is necessary for a potential carry */
 438              /* ensure we do not use a too large precision */
 439              if (prec_x > prec_u)
 440                prec_x = prec_u;
 441              if (prec_x > prec)
 442                mpfr_prec_round (x, prec_x, GMP_RNDN);
 443            }
 444
 445          inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */
 446
 447          /* in case u=0, ensure that rnd_u rounds x away from zero */
 448          if (mpfr_sgn (u) == 0)
 449            rnd_u = (mpfr_sgn (x) > 0) ? GMP_RNDU : GMP_RNDD;
 450          inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */
 451
 452          ok = inexact == 0 ||
 453            mpfr_can_round (u, prec_u - 3, rnd_u, GMP_RNDZ,
 454                            prec_re + (rnd_re == GMP_RNDN));
 455          /* this ensures both we can round correctly and determine the correct
 456             inexact flag (for rounding to nearest) */
 457      }
 458    while (ok == 0);
 459
 460    /* if inexact is zero, then u is exactly ac-bd, otherwise fix the sign
 461       of the inexact flag for u, which was rounded away from ac-bd */
 462    if (inexact != 0)
 463      inexact = mpfr_sgn (u);
 464
 465    if (mul_i == 0)
 466      {
 467        inex_re = mpfr_set (MPC_RE(result), u, MPC_RND_RE(rnd));
 468        if (inex_re == 0)
 469          {
 470            inex_re = inexact; /* u is rounded away from 0 */
 471            inex_im = mpfr_add (MPC_IM(result), v, w, MPC_RND_IM(rnd));
 472          }
 473        else
 474          inex_im = mpfr_add (MPC_IM(result), v, w, MPC_RND_IM(rnd));
 475      }
 476    else if (mul_i == 1) /* (x+i*y)/i = y - i*x */
 477      {
 478        inex_im = mpfr_neg (MPC_IM(result), u, MPC_RND_IM(rnd));
 479        if (inex_im == 0)
 480          {
 481            inex_im = -inexact; /* u is rounded away from 0 */
 482            inex_re = mpfr_add (MPC_RE(result), v, w, MPC_RND_RE(rnd));
 483          }
 484        else
 485          inex_re = mpfr_add (MPC_RE(result), v, w, MPC_RND_RE(rnd));
 486      }
 487    else /* mul_i = 2, z/i^2 = -z */
 488      {
 489        inex_re = mpfr_neg (MPC_RE(result), u, MPC_RND_RE(rnd));
 490        if (inex_re == 0)
 491          {
 492            inex_re = -inexact; /* u is rounded away from 0 */
 493            inex_im = -mpfr_add (MPC_IM(result), v, w,
 494                                 INV_RND(MPC_RND_IM(rnd)));
 495            mpfr_neg (MPC_IM(result), MPC_IM(result), MPC_RND_IM(rnd));
 496          }
 497        else
 498          {
 499            inex_im = -mpfr_add (MPC_IM(result), v, w,
 500                                 INV_RND(MPC_RND_IM(rnd)));
 501            mpfr_neg (MPC_IM(result), MPC_IM(result), MPC_RND_IM(rnd));
 502          }
 503      }
 504
 505    mpc_set (rop, result, MPC_RNDNN);
 506
 507    mpfr_clear (u);
 508    mpfr_clear (v);
 509    mpfr_clear (w);
 510    mpfr_clear (x);
 511    if (overlap)
 512       mpc_clear (result);
 513
 514    return MPC_INEX(inex_re, inex_im);
 515 }