contrib/mpc/src/div.c

   1 /* mpc_div -- Divide two complex numbers.
   2
   3 Copyright (C) 2002, 2003, 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 "mpc-impl.h"
  22
  23 /* this routine deals with the case where w is zero */
  24 static int
  25 mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
  26 /* Assumes w==0, implementation according to C99 G.5.1.8 */
  27 {
  28    int sign = MPFR_SIGNBIT (mpc_realref (w));
  29    mpfr_t infty;
  30
  31    mpfr_init2 (infty, MPFR_PREC_MIN);
  32    mpfr_set_inf (infty, sign);
  33    mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd));
  34    mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd));
  35    mpfr_clear (infty);
  36    return MPC_INEX (0, 0); /* exact */
  37 }
  38
  39 /* this routine deals with the case where z is infinite and w finite */
  40 static int
  41 mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
  42 /* Assumes w finite and non-zero and z infinite; implementation
  43    according to C99 G.5.1.8                                     */
  44 {
  45    int a, b, x, y;
  46
  47    a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0);
  48    b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0);
  49
  50    /* a is -1 if Re(z) = -Inf, 1 if Re(z) = +Inf, 0 if Re(z) is finite
  51       b is -1 if Im(z) = -Inf, 1 if Im(z) = +Inf, 0 if Im(z) is finite */
  52
  53    /* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */
  54    /* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */
  55    if (a == 0 || b == 0) {
  56      /* only one of a or b can be zero, since z is infinite */
  57       x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w));
  58       y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w));
  59    }
  60    else {
  61       /* Both parts of z are infinite; x could be determined by sign
  62          considerations and comparisons. Since operations with non-finite
  63          numbers are not considered time-critical, we let mpfr do the work. */
  64       mpfr_t sign;
  65
  66       mpfr_init2 (sign, 2);
  67       /* This is enough to determine the sign of sums and differences. */
  68
  69       if (a == 1)
  70          if (b == 1) {
  71             mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
  72             x = MPC_MPFR_SIGN (sign);
  73             mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
  74             y = MPC_MPFR_SIGN (sign);
  75          }
  76          else { /* b == -1 */
  77             mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
  78             x = MPC_MPFR_SIGN (sign);
  79             mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
  80             y = -MPC_MPFR_SIGN (sign);
  81          }
  82       else /* a == -1 */
  83          if (b == 1) {
  84             mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN);
  85             x = MPC_MPFR_SIGN (sign);
  86             mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
  87             y = MPC_MPFR_SIGN (sign);
  88          }
  89          else { /* b == -1 */
  90             mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
  91             x = -MPC_MPFR_SIGN (sign);
  92             mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN);
  93             y = MPC_MPFR_SIGN (sign);
  94          }
  95       mpfr_clear (sign);
  96    }
  97
  98    if (x == 0)
  99       mpfr_set_nan (mpc_realref (rop));
 100    else
 101       mpfr_set_inf (mpc_realref (rop), x);
 102    if (y == 0)
 103       mpfr_set_nan (mpc_imagref (rop));
 104    else
 105       mpfr_set_inf (mpc_imagref (rop), y);
 106
 107    return MPC_INEX (0, 0); /* exact */
 108 }
 109
 110
 111 /* this routine deals with the case where z if finite and w infinite */
 112 static int
 113 mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
 114 /* Assumes z finite and w infinite; implementation according to
 115    C99 G.5.1.8                                                  */
 116 {
 117    mpfr_t c, d, a, b, x, y, zero;
 118
 119    mpfr_init2 (c, 2); /* needed to hold a signed zero, +1 or -1 */
 120    mpfr_init2 (d, 2);
 121    mpfr_init2 (x, 2);
 122    mpfr_init2 (y, 2);
 123    mpfr_init2 (zero, 2);
 124    mpfr_set_ui (zero, 0ul, GMP_RNDN);
 125    mpfr_init2 (a, mpfr_get_prec (mpc_realref (z)));
 126    mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z)));
 127
 128    mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), GMP_RNDN);
 129    MPFR_COPYSIGN (c, c, mpc_realref (w), GMP_RNDN);
 130    mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), GMP_RNDN);
 131    MPFR_COPYSIGN (d, d, mpc_imagref (w), GMP_RNDN);
 132
 133    mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */
 134    mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN);
 135    mpfr_add (x, a, b, GMP_RNDN);
 136
 137    mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN);
 138    mpfr_mul (a, mpc_realref (z), d, GMP_RNDN);
 139    mpfr_sub (y, b, a, GMP_RNDN);
 140
 141    MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN);
 142    MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN);
 143
 144    mpfr_clear (c);
 145    mpfr_clear (d);
 146    mpfr_clear (x);
 147    mpfr_clear (y);
 148    mpfr_clear (zero);
 149    mpfr_clear (a);
 150    mpfr_clear (b);
 151
 152    return MPC_INEX (0, 0); /* exact */
 153 }
 154
 155
 156 static int
 157 mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
 158 /* Assumes z finite and w finite and non-zero, with imaginary part
 159    of w a signed zero.                                             */
 160 {
 161    int inex_re, inex_im;
 162    /* save signs of operands in case there are overlaps */
 163    int zrs = MPFR_SIGNBIT (mpc_realref (z));
 164    int zis = MPFR_SIGNBIT (mpc_imagref (z));
 165    int wrs = MPFR_SIGNBIT (mpc_realref (w));
 166    int wis = MPFR_SIGNBIT (mpc_imagref (w));
 167
 168    /* warning: rop may overlap with z,w so treat the imaginary part first */
 169    inex_im = mpfr_div (mpc_imagref(rop), mpc_imagref(z), mpc_realref(w), MPC_RND_IM(rnd));
 170    inex_re = mpfr_div (mpc_realref(rop), mpc_realref(z), mpc_realref(w), MPC_RND_RE(rnd));
 171
 172    /* correct signs of zeroes if necessary, which does not affect the
 173       inexact flags                                                    */
 174    if (mpfr_zero_p (mpc_realref (rop)))
 175       mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
 176          GMP_RNDN); /* exact */
 177    if (mpfr_zero_p (mpc_imagref (rop)))
 178       mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
 179          GMP_RNDN);
 180
 181    return MPC_INEX(inex_re, inex_im);
 182 }
 183
 184
 185 static int
 186 mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
 187 /* Assumes z finite and w finite and non-zero, with real part
 188    of w a signed zero.                                        */
 189 {
 190    int inex_re, inex_im;
 191    int overlap = (rop == z) || (rop == w);
 192    int imag_z = mpfr_zero_p (mpc_realref (z));
 193    mpfr_t wloc;
 194    mpc_t tmprop;
 195    mpc_ptr dest = (overlap) ? tmprop : rop;
 196    /* save signs of operands in case there are overlaps */
 197    int zrs = MPFR_SIGNBIT (mpc_realref (z));
 198    int zis = MPFR_SIGNBIT (mpc_imagref (z));
 199    int wrs = MPFR_SIGNBIT (mpc_realref (w));
 200    int wis = MPFR_SIGNBIT (mpc_imagref (w));
 201
 202    if (overlap)
 203       mpc_init3 (tmprop, MPC_PREC_RE (rop), MPC_PREC_IM (rop));
 204
 205    wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */
 206    inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd));
 207    mpfr_neg (wloc, wloc, GMP_RNDN);
 208    /* changes the sign only in wloc, not in w; no need to correct later */
 209    inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd));
 210
 211    if (overlap) {
 212       /* Note: we could use mpc_swap here, but this might cause problems
 213          if rop and tmprop have been allocated using different methods, since
 214          it will swap the significands of rop and tmprop. See
 215          http://lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */
 216       mpc_set (rop, tmprop, MPC_RNDNN); /* exact */
 217       mpc_clear (tmprop);
 218    }
 219
 220    /* correct signs of zeroes if necessary, which does not affect the
 221       inexact flags                                                    */
 222    if (mpfr_zero_p (mpc_realref (rop)))
 223       mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
 224          GMP_RNDN); /* exact */
 225    if (imag_z)
 226       mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
 227          GMP_RNDN);
 228
 229    return MPC_INEX(inex_re, inex_im);
 230 }
 231
 232
 233 int
 234 mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
 235 {
 236    int ok_re = 0, ok_im = 0;
 237    mpc_t res, c_conj;
 238    mpfr_t q;
 239    mpfr_prec_t prec;
 240    int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
 241    int underflow_norm, overflow_norm, underflow_prod, overflow_prod;
 242    int underflow_re = 0, overflow_re = 0, underflow_im = 0, overflow_im = 0;
 243    mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd);
 244    int saved_underflow, saved_overflow;
 245    int tmpsgn;
 246
 247    /* According to the C standard G.3, there are three types of numbers:   */
 248    /* finite (both parts are usual real numbers; contains 0), infinite     */
 249    /* (at least one part is a real infinity) and all others; the latter    */
 250    /* are numbers containing a nan, but no infinity, and could reasonably  */
 251    /* be called nan.                                                       */
 252    /* By G.5.1.4, infinite/finite=infinite; finite/infinite=0;             */
 253    /* all other divisions that are not finite/finite return nan+i*nan.     */
 254    /* Division by 0 could be handled by the following case of division by  */
 255    /* a real; we handle it separately instead.                             */
 256    if (mpc_zero_p (c))
 257       return mpc_div_zero (a, b, c, rnd);
 258    else if (mpc_inf_p (b) && mpc_fin_p (c))
 259          return mpc_div_inf_fin (a, b, c);
 260    else if (mpc_fin_p (b) && mpc_inf_p (c))
 261          return mpc_div_fin_inf (a, b, c);
 262    else if (!mpc_fin_p (b) || !mpc_fin_p (c)) {
 263       mpc_set_nan (a);
 264       return MPC_INEX (0, 0);
 265    }
 266    else if (mpfr_zero_p(mpc_imagref(c)))
 267       return mpc_div_real (a, b, c, rnd);
 268    else if (mpfr_zero_p(mpc_realref(c)))
 269       return mpc_div_imag (a, b, c, rnd);
 270
 271    prec = MPC_MAX_PREC(a);
 272
 273    mpc_init2 (res, 2);
 274    mpfr_init (q);
 275
 276    /* create the conjugate of c in c_conj without allocating new memory */
 277    mpc_realref (c_conj)[0] = mpc_realref (c)[0];
 278    mpc_imagref (c_conj)[0] = mpc_imagref (c)[0];
 279    MPFR_CHANGE_SIGN (mpc_imagref (c_conj));
 280
 281    /* save the underflow or overflow flags from MPFR */
 282    saved_underflow = mpfr_underflow_p ();
 283    saved_overflow = mpfr_overflow_p ();
 284
 285    do {
 286       loops ++;
 287       prec += loops <= 2 ? mpc_ceil_log2 (prec) + 5 : prec / 2;
 288
 289       mpc_set_prec (res, prec);
 290       mpfr_set_prec (q, prec);
 291
 292       /* first compute norm(c) */
 293       mpfr_clear_underflow ();
 294       mpfr_clear_overflow ();
 295       inexact_norm = mpc_norm (q, c, GMP_RNDU);
 296       underflow_norm = mpfr_underflow_p ();
 297       overflow_norm = mpfr_overflow_p ();
 298       if (underflow_norm)
 299          mpfr_set_ui (q, 0ul, GMP_RNDN);
 300          /* to obtain divisions by 0 later on */
 301
 302       /* now compute b*conjugate(c) */
 303       mpfr_clear_underflow ();
 304       mpfr_clear_overflow ();
 305       inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ);
 306       inexact_re = MPC_INEX_RE (inexact_prod);
 307       inexact_im = MPC_INEX_IM (inexact_prod);
 308       underflow_prod = mpfr_underflow_p ();
 309       overflow_prod = mpfr_overflow_p ();
 310          /* unfortunately, does not distinguish between under-/overflow
 311             in real or imaginary parts
 312             hopefully, the side-effects of mpc_mul do indeed raise the
 313             mpfr exceptions */
 314       if (overflow_prod) {
 315          int isinf = 0;
 316          tmpsgn = mpfr_sgn (mpc_realref(res));
 317          if (tmpsgn > 0)
 318            {
 319              mpfr_nextabove (mpc_realref(res));
 320              isinf = mpfr_inf_p (mpc_realref(res));
 321              mpfr_nextbelow (mpc_realref(res));
 322            }
 323          else if (tmpsgn < 0)
 324            {
 325              mpfr_nextbelow (mpc_realref(res));
 326              isinf = mpfr_inf_p (mpc_realref(res));
 327              mpfr_nextabove (mpc_realref(res));
 328            }
 329          if (isinf)
 330            {
 331              mpfr_set_inf (mpc_realref(res), tmpsgn);
 332              overflow_re = 1;
 333            }
 334          tmpsgn = mpfr_sgn (mpc_imagref(res));
 335          isinf = 0;
 336          if (tmpsgn > 0)
 337            {
 338              mpfr_nextabove (mpc_imagref(res));
 339              isinf = mpfr_inf_p (mpc_imagref(res));
 340              mpfr_nextbelow (mpc_imagref(res));
 341            }
 342          else if (tmpsgn < 0)
 343            {
 344              mpfr_nextbelow (mpc_imagref(res));
 345              isinf = mpfr_inf_p (mpc_imagref(res));
 346              mpfr_nextabove (mpc_imagref(res));
 347            }
 348          if (isinf)
 349            {
 350              mpfr_set_inf (mpc_imagref(res), tmpsgn);
 351              overflow_im = 1;
 352            }
 353          mpc_set (a, res, rnd);
 354          goto end;
 355       }
 356
 357       /* divide the product by the norm */
 358       if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) {
 359          /* The division has good chances to be exact in at least one part.  */
 360          /* Since this can cause problems when not rounding to the nearest,  */
 361          /* we use the division code of mpfr, which handles the situation.   */
 362          mpfr_clear_underflow ();
 363          mpfr_clear_overflow ();
 364          inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ);
 365          underflow_re = mpfr_underflow_p ();
 366          overflow_re = mpfr_overflow_p ();
 367          ok_re = !inexact_re || underflow_re || overflow_re
 368                  || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN,
 369                     GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN));
 370
 371          if (ok_re) /* compute imaginary part */ {
 372             mpfr_clear_underflow ();
 373             mpfr_clear_overflow ();
 374             inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ);
 375             underflow_im = mpfr_underflow_p ();
 376             overflow_im = mpfr_overflow_p ();
 377             ok_im = !inexact_im || underflow_im || overflow_im
 378                     || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN,
 379                        GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN));
 380          }
 381       }
 382       else {
 383          /* The division is inexact, so for efficiency reasons we invert q */
 384          /* only once and multiply by the inverse. */
 385          if (mpfr_ui_div (q, 1ul, q, GMP_RNDZ) || inexact_norm) {
 386              /* if 1/q is inexact, the approximations of the real and
 387                 imaginary part below will be inexact, unless RE(res)
 388                 or IM(res) is zero */
 389              inexact_re |= ~mpfr_zero_p (mpc_realref (res));
 390              inexact_im |= ~mpfr_zero_p (mpc_imagref (res));
 391          }
 392          mpfr_clear_underflow ();
 393          mpfr_clear_overflow ();
 394          inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ);
 395          underflow_re = mpfr_underflow_p ();
 396          overflow_re = mpfr_overflow_p ();
 397          ok_re = !inexact_re || underflow_re || overflow_re
 398                  || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN,
 399                     GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN));
 400
 401          if (ok_re) /* compute imaginary part */ {
 402             mpfr_clear_underflow ();
 403             mpfr_clear_overflow ();
 404             inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ);
 405             underflow_im = mpfr_underflow_p ();
 406             overflow_im = mpfr_overflow_p ();
 407             ok_im = !inexact_im || underflow_im || overflow_im
 408                     || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN,
 409                        GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN));
 410          }
 411       }
 412    } while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm
 413                                && !underflow_prod && !overflow_prod);
 414
 415    inex = mpc_set (a, res, rnd);
 416    inexact_re = MPC_INEX_RE (inex);
 417    inexact_im = MPC_INEX_IM (inex);
 418
 419  end:
 420    /* fix values and inexact flags in case of overflow/underflow */
 421    /* FIXME: heuristic, certainly does not cover all cases */
 422    if (overflow_re || (underflow_norm && !underflow_prod)) {
 423       mpfr_set_inf (mpc_realref (a), mpfr_sgn (mpc_realref (res)));
 424       inexact_re = mpfr_sgn (mpc_realref (res));
 425    }
 426    else if (underflow_re || (overflow_norm && !overflow_prod)) {
 427       inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1;
 428       mpfr_set_zero (mpc_realref (a), -inexact_re);
 429    }
 430    if (overflow_im || (underflow_norm && !underflow_prod)) {
 431       mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res)));
 432       inexact_im = mpfr_sgn (mpc_imagref (res));
 433    }
 434    else if (underflow_im || (overflow_norm && !overflow_prod)) {
 435       inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1;
 436       mpfr_set_zero (mpc_imagref (a), -inexact_im);
 437    }
 438
 439    mpc_clear (res);
 440    mpfr_clear (q);
 441
 442    /* restore underflow and overflow flags from MPFR */
 443    if (saved_underflow)
 444      mpfr_set_underflow ();
 445    if (saved_overflow)
 446      mpfr_set_overflow ();
 447
 448    return MPC_INEX (inexact_re, inexact_im);
 449 }