src/div.c

   1 /* mpc_div -- Divide two complex numbers.
   2
   3 Copyright (C) 2002, 2003, 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
  25 mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
  26 {
  27    /* Assumes w==0, implementation according to C99 G.5.1.8 */
  28    int sign = MPFR_SIGNBIT (MPC_RE (w));
  29    mpfr_t infty;
  30    mpfr_set_inf (infty, sign);
  31    mpfr_mul (MPC_RE (a), infty, MPC_RE (z), MPC_RND_RE (rnd));
  32    mpfr_mul (MPC_IM (a), infty, MPC_IM (z), MPC_RND_IM (rnd));
  33    return MPC_INEX (0, 0); /* exact */
  34 }
  35
  36 static int
  37 mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
  38 {
  39    /* Assumes w finite and non-zero and z infinite; implementation
  40       according to C99 G.5.1.8                                     */
  41    int a, b, x, y;
  42
  43    a = (mpfr_inf_p (MPC_RE (z)) ? MPFR_SIGNBIT (MPC_RE (z)) : 0);
  44    b = (mpfr_inf_p (MPC_IM (z)) ? MPFR_SIGNBIT (MPC_IM (z)) : 0);
  45
  46    /* x = MPC_MPFR_SIGN (a * MPC_RE (w) + b * MPC_IM (w)) */
  47    /* y = MPC_MPFR_SIGN (b * MPC_RE (w) - a * MPC_IM (w)) */
  48    if (a == 0 || b == 0) {
  49       x = a * MPC_MPFR_SIGN (MPC_RE (w)) + b * MPC_MPFR_SIGN (MPC_IM (w));
  50       y = b * MPC_MPFR_SIGN (MPC_RE (w)) - a * MPC_MPFR_SIGN (MPC_IM (w));
  51    }
  52    else {
  53       /* Both parts of z are infinite; x could be determined by sign
  54          considerations and comparisons. Since operations with non-finite
  55          numbers are not considered time-critical, we let mpfr do the work. */
  56       mpfr_t sign;
  57       mpfr_init2 (sign, 2);
  58          /* This is enough to determine the sign of sums and differences. */
  59
  60       if (a == 1)
  61          if (b == 1) {
  62             mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
  63             x = MPC_MPFR_SIGN (sign);
  64             mpfr_sub (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
  65             y = MPC_MPFR_SIGN (sign);
  66          }
  67          else { /* b == -1 */
  68             mpfr_sub (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
  69             x = MPC_MPFR_SIGN (sign);
  70             mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
  71             y = -MPC_MPFR_SIGN (sign);
  72          }
  73       else /* a == -1 */
  74          if (b == 1) {
  75             mpfr_sub (sign, MPC_IM (w), MPC_RE (w), GMP_RNDN);
  76             x = MPC_MPFR_SIGN (sign);
  77             mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
  78             y = MPC_MPFR_SIGN (sign);
  79          }
  80          else { /* b == -1 */
  81             mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
  82             x = -MPC_MPFR_SIGN (sign);
  83             mpfr_sub (sign, MPC_IM (w), MPC_RE (w), GMP_RNDN);
  84             y = MPC_MPFR_SIGN (sign);
  85          }
  86       mpfr_clear (sign);
  87    }
  88
  89    if (x == 0)
  90       mpfr_set_nan (MPC_RE (rop));
  91    else
  92       mpfr_set_inf (MPC_RE (rop), x);
  93    if (y == 0)
  94       mpfr_set_nan (MPC_IM (rop));
  95    else
  96       mpfr_set_inf (MPC_IM (rop), y);
  97
  98    return MPC_INEX (0, 0); /* exact */
  99 }
 100
 101
 102 static int
 103 mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
 104 {
 105    /* Assumes z finite and w infinite; implementation according to
 106       C99 G.5.1.8                                                  */
 107    mpfr_t c, d, a, b, x, y, zero;
 108
 109    mpfr_init2 (c, 2); /* needed to hold a signed zero, +1 or -1 */
 110    mpfr_init2 (d, 2);
 111    mpfr_init2 (x, 2);
 112    mpfr_init2 (y, 2);
 113    mpfr_init2 (zero, 2);
 114    mpfr_set_ui (zero, 0ul, GMP_RNDN);
 115    mpfr_init2 (a, mpfr_get_prec (MPC_RE (z)));
 116    mpfr_init2 (b, mpfr_get_prec (MPC_IM (z)));
 117
 118    mpfr_set_ui (c, (mpfr_inf_p (MPC_RE (w)) ? 1 : 0), GMP_RNDN);
 119    MPFR_COPYSIGN (c, c, MPC_RE (w), GMP_RNDN);
 120    mpfr_set_ui (d, (mpfr_inf_p (MPC_IM (w)) ? 1 : 0), GMP_RNDN);
 121    MPFR_COPYSIGN (d, d, MPC_IM (w), GMP_RNDN);
 122
 123    mpfr_mul (a, MPC_RE (z), c, GMP_RNDN); /* exact */
 124    mpfr_mul (b, MPC_IM (z), d, GMP_RNDN);
 125    mpfr_add (x, a, b, GMP_RNDN);
 126
 127    mpfr_mul (b, MPC_IM (z), c, GMP_RNDN);
 128    mpfr_mul (a, MPC_RE (z), d, GMP_RNDN);
 129    mpfr_sub (y, b, a, GMP_RNDN);
 130
 131    MPFR_COPYSIGN (MPC_RE (rop), zero, x, GMP_RNDN);
 132    MPFR_COPYSIGN (MPC_IM (rop), zero, y, GMP_RNDN);
 133
 134    mpfr_clear (c);
 135    mpfr_clear (d);
 136    mpfr_clear (x);
 137    mpfr_clear (y);
 138    mpfr_clear (zero);
 139    mpfr_clear (a);
 140    mpfr_clear (b);
 141
 142    return MPC_INEX (0, 0); /* exact */
 143 }
 144
 145
 146 int
 147 mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
 148 {
 149    int ok_re = 0, ok_im = 0;
 150    mpc_t res, c_conj;
 151    mpfr_t q;
 152    mp_prec_t prec;
 153    int inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
 154
 155    /* save signs of operands in case there are overlaps */
 156    int brs = MPFR_SIGNBIT (MPC_RE (b));
 157    int bis = MPFR_SIGNBIT (MPC_IM (b));
 158    int crs = MPFR_SIGNBIT (MPC_RE (c));
 159    int cis = MPFR_SIGNBIT (MPC_IM (c));
 160
 161    /* According to the C standard G.3, there are three types of numbers:   */
 162    /* finite (both parts are usual real numbers; contains 0), infinite     */
 163    /* (at least one part is a real infinity) and all others; the latter    */
 164    /* are numbers containing a nan, but no infinity, and could reasonably  */
 165    /* be called nan.                                                       */
 166    /* By G.5.1.4, infinite/finite=infinite; finite/infinite=0;             */
 167    /* all other divisions that are not finite/finite return nan+i*nan.     */
 168    /* Division by 0 could be handled by the following case of division by  */
 169    /* a real; we handle it separately instead.                             */
 170    if (mpc_zero_p (c))
 171       return mpc_div_zero (a, b, c, rnd);
 172    else {
 173       if (mpc_inf_p (b) && mpc_fin_p (c))
 174          return mpc_div_inf_fin (a, b, c);
 175       else if (mpc_fin_p (b) && mpc_inf_p (c))
 176          return mpc_div_fin_inf (a, b, c);
 177       else if (!mpc_fin_p (b) || !mpc_fin_p (c)) {
 178          mpfr_set_nan (MPC_RE (a));
 179          mpfr_set_nan (MPC_IM (a));
 180          return MPC_INEX (0, 0);
 181       }
 182    }
 183
 184    /* check for real divisor */
 185    if (mpfr_zero_p(MPC_IM(c))) /* (re_b+i*im_b)/c = re_b/c + i * (im_b/c) */
 186    {
 187       /* warning: a may overlap with b,c so treat the imaginary part first */
 188       inexact_im = mpfr_div (MPC_IM(a), MPC_IM(b), MPC_RE(c), MPC_RND_IM(rnd));
 189       inexact_re = mpfr_div (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd));
 190
 191       /* correct signs of zeroes if necessary, which does not affect the
 192          inexact flags                                                    */
 193       if (mpfr_zero_p (MPC_RE (a)))
 194          mpfr_setsign (MPC_RE (a), MPC_RE (a), (brs != crs && bis != cis),
 195             GMP_RNDN); /* exact */
 196       if (mpfr_zero_p (MPC_IM (a)))
 197          mpfr_setsign (MPC_IM (a), MPC_IM (a), (bis != crs && brs == cis),
 198             GMP_RNDN);
 199
 200       return MPC_INEX(inexact_re, inexact_im);
 201    }
 202
 203    /* check for purely imaginary divisor */
 204    if (mpfr_zero_p(MPC_RE(c)))
 205    {
 206       /* (re_b+i*im_b)/(i*c) = im_b/c - i * (re_b/c) */
 207       int overlap = (a == b) || (a == c);
 208       int imag_b = mpfr_zero_p (MPC_RE (b));
 209       mpfr_t cloc;
 210       mpc_t tmpa;
 211       mpc_ptr dest = (overlap) ? tmpa : a;
 212
 213       if (overlap)
 214          mpc_init3 (tmpa, MPFR_PREC (MPC_RE (a)), MPFR_PREC (MPC_IM (a)));
 215
 216       cloc[0] = MPC_IM(c)[0]; /* copies mpfr struct IM(c) into cloc */
 217       inexact_re = mpfr_div (MPC_RE(dest), MPC_IM(b), cloc, MPC_RND_RE(rnd));
 218       mpfr_neg (cloc, cloc, GMP_RNDN);
 219       /* changes the sign only in cloc, not in c; no need to correct later */
 220       inexact_im = mpfr_div (MPC_IM(dest), MPC_RE(b), cloc, MPC_RND_IM(rnd));
 221
 222       if (overlap)
 223         {
 224           /* Note: we could use mpc_swap here, but this might cause problems
 225              if a and tmpa have been allocated using different methods, since
 226              it will swap the significands of a and tmpa. See http://
 227          lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */
 228           mpc_set (a, tmpa, MPC_RNDNN); /* exact */
 229           mpc_clear (tmpa);
 230         }
 231
 232       /* correct signs of zeroes if necessary, which does not affect the
 233          inexact flags                                                    */
 234       if (mpfr_zero_p (MPC_RE (a)))
 235          mpfr_setsign (MPC_RE (a), MPC_RE (a), (brs != crs && bis != cis),
 236             GMP_RNDN); /* exact */
 237       if (imag_b)
 238          mpfr_setsign (MPC_IM (a), MPC_IM (a), (bis != crs && brs == cis),
 239             GMP_RNDN);
 240
 241       return MPC_INEX(inexact_re, inexact_im);
 242    }
 243
 244    prec = MPC_MAX_PREC(a);
 245
 246    mpc_init2 (res, 2);
 247    mpfr_init (q);
 248
 249    /* create the conjugate of c in c_conj without allocating new memory */
 250    MPC_RE (c_conj)[0] = MPC_RE (c)[0];
 251    MPC_IM (c_conj)[0] = MPC_IM (c)[0];
 252    MPFR_CHANGE_SIGN (MPC_IM (c_conj));
 253
 254    do
 255    {
 256       loops ++;
 257       prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2;
 258
 259       mpc_set_prec (res, prec);
 260       mpfr_set_prec (q, prec);
 261
 262       /* first compute norm(c)^2 */
 263       inexact_norm = mpc_norm (q, c, GMP_RNDD);
 264
 265       /* now compute b*conjugate(c) */
 266       /* We need rounding away from zero for both the real and the imagin-  */
 267       /* ary part; then the final result is also rounded away from zero.    */
 268       /* The error is less than 1 ulp. Since this is not implemented in     */
 269       /* mpfr, we round towards zero and add 1 ulp to the absolute values   */
 270       /* if they are not exact. */
 271       inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ);
 272       inexact_re = MPC_INEX_RE (inexact_prod);
 273       inexact_im = MPC_INEX_IM (inexact_prod);
 274       if (inexact_re != 0)
 275          mpfr_add_one_ulp (MPC_RE (res), GMP_RNDN);
 276       if (inexact_im != 0)
 277          mpfr_add_one_ulp (MPC_IM (res), GMP_RNDN);
 278
 279       /* divide the product by the norm */
 280       if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0))
 281       {
 282          /* The division has good chances to be exact in at least one part.   */
 283          /* Since this can cause problems when not rounding to the nearest,   */
 284          /* we use the division code of mpfr, which handles the situation.    */
 285          if (MPFR_SIGN (MPC_RE (res)) > 0)
 286          {
 287             inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDU);
 288             ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDU,
 289                   MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
 290          }
 291          else
 292          {
 293             inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDD);
 294             ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDD,
 295                   MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
 296          }
 297
 298          if (ok_re || !inexact_re) /* compute imaginary part */
 299          {
 300             if (MPFR_SIGN (MPC_IM (res)) > 0)
 301             {
 302                inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDU);
 303                ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDU,
 304                      MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
 305             }
 306             else
 307             {
 308                inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDD);
 309                ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDD,
 310                      MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
 311             }
 312          }
 313       }
 314       else
 315       {
 316          /* The division is inexact, so for efficiency reasons we invert q */
 317          /* only once and multiply by the inverse. */
 318          /* We do not decide about the sign of the difference. */
 319          if (mpfr_ui_div (q, 1, q, GMP_RNDU) || inexact_norm)
 320            {
 321              /* if 1/q is inexact, the approximations of the real and
 322                 imaginary part below will be inexact, unless RE(res)
 323                 or IM(res) is zero */
 324              inexact_re = inexact_re || !mpfr_zero_p (MPC_RE (res));
 325              inexact_im = inexact_im || !mpfr_zero_p (MPC_IM (res));
 326            }
 327          if (MPFR_SIGN (MPC_RE (res)) > 0)
 328          {
 329            inexact_re = mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDU)
 330              || inexact_re;
 331            ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDU,
 332                                    MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
 333          }
 334          else
 335          {
 336            inexact_re = mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDD)
 337              || inexact_re;
 338            ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDD,
 339                                    MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
 340          }
 341
 342          if (ok_re) /* compute imaginary part */
 343          {
 344             if (MPFR_SIGN (MPC_IM (res)) > 0)
 345             {
 346               inexact_im = mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDU)
 347                 || inexact_im;
 348                ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDU,
 349                                        MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
 350             }
 351             else
 352             {
 353               inexact_im = mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDD)
 354                 || inexact_im;
 355               ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDD,
 356                                       MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
 357             }
 358          }
 359       }
 360    }
 361    while ((!ok_re && inexact_re) || (!ok_im && inexact_im));
 362
 363    mpc_set (a, res, rnd);
 364
 365    mpc_clear (res);
 366    mpfr_clear (q);
 367
 368    return (MPC_INEX (inexact_re, inexact_im));
 369       /* Only exactness vs. inexactness is tested, not the sign. */
 370 }