contrib/mpc/src/log.c

   1 /* mpc_log -- Take the logarithm of a complex number.
   2
   3 Copyright (C) 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 int
  25 mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){
  26    int ok, underflow = 0;
  27    mpfr_srcptr x, y;
  28    mpfr_t v, w;
  29    mpfr_prec_t prec;
  30    int loops;
  31    int re_cmp, im_cmp;
  32    int inex_re, inex_im;
  33    int err;
  34    mpfr_exp_t expw;
  35    int sgnw;
  36
  37    /* special values: NaN and infinities */
  38    if (!mpc_fin_p (op)) {
  39       if (mpfr_nan_p (mpc_realref (op))) {
  40          if (mpfr_inf_p (mpc_imagref (op)))
  41             mpfr_set_inf (mpc_realref (rop), +1);
  42          else
  43             mpfr_set_nan (mpc_realref (rop));
  44          mpfr_set_nan (mpc_imagref (rop));
  45          inex_im = 0; /* Inf/NaN is exact */
  46       }
  47       else if (mpfr_nan_p (mpc_imagref (op))) {
  48          if (mpfr_inf_p (mpc_realref (op)))
  49             mpfr_set_inf (mpc_realref (rop), +1);
  50          else
  51             mpfr_set_nan (mpc_realref (rop));
  52          mpfr_set_nan (mpc_imagref (rop));
  53          inex_im = 0; /* Inf/NaN is exact */
  54       }
  55       else /* We have an infinity in at least one part. */ {
  56          inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),
  57                                MPC_RND_IM (rnd));
  58          mpfr_set_inf (mpc_realref (rop), +1);
  59       }
  60       return MPC_INEX(0, inex_im);
  61    }
  62
  63    /* special cases: real and purely imaginary numbers */
  64    re_cmp = mpfr_cmp_ui (mpc_realref (op), 0);
  65    im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0);
  66    if (im_cmp == 0) {
  67       if (re_cmp == 0) {
  68          inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),
  69                                MPC_RND_IM (rnd));
  70          mpfr_set_inf (mpc_realref (rop), -1);
  71          inex_re = 0; /* -Inf is exact */
  72       }
  73       else if (re_cmp > 0) {
  74          inex_re = mpfr_log (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
  75          inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
  76       }
  77       else {
  78          /* op = x + 0*y; let w = -x = |x| */
  79          int negative_zero;
  80          mpfr_rnd_t rnd_im;
  81
  82          negative_zero = mpfr_signbit (mpc_imagref (op));
  83          if (negative_zero)
  84             rnd_im = INV_RND (MPC_RND_IM (rnd));
  85          else
  86             rnd_im = MPC_RND_IM (rnd);
  87          w [0] = *mpc_realref (op);
  88          MPFR_CHANGE_SIGN (w);
  89          inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd));
  90          inex_im = mpfr_const_pi (mpc_imagref (rop), rnd_im);
  91          if (negative_zero) {
  92             mpc_conj (rop, rop, MPC_RNDNN);
  93             inex_im = -inex_im;
  94          }
  95       }
  96       return MPC_INEX(inex_re, inex_im);
  97    }
  98    else if (re_cmp == 0) {
  99       if (im_cmp > 0) {
 100          inex_re = mpfr_log (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd));
 101          inex_im = mpfr_const_pi (mpc_imagref (rop), MPC_RND_IM (rnd));
 102          /* division by 2 does not change the ternary flag */
 103          mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
 104       }
 105       else {
 106          w [0] = *mpc_imagref (op);
 107          MPFR_CHANGE_SIGN (w);
 108          inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd));
 109          inex_im = mpfr_const_pi (mpc_imagref (rop), INV_RND (MPC_RND_IM (rnd)));
 110          /* division by 2 does not change the ternary flag */
 111          mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
 112          mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN);
 113          inex_im = -inex_im; /* negate the ternary flag */
 114       }
 115       return MPC_INEX(inex_re, inex_im);
 116    }
 117
 118    prec = MPC_PREC_RE(rop);
 119    mpfr_init2 (w, 2);
 120    /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x)   */
 121    /* loop for the real part: 1/2 log (x^2 + y^2), fast, but unsafe */
 122    /* implementation                                                */
 123    ok = 0;
 124    for (loops = 1; !ok && loops <= 2; loops++) {
 125       prec += mpc_ceil_log2 (prec) + 4;
 126       mpfr_set_prec (w, prec);
 127
 128       mpc_abs (w, op, GMP_RNDN);
 129          /* error 0.5 ulp */
 130       if (mpfr_inf_p (w))
 131          /* intermediate overflow; the logarithm may be representable.
 132             Intermediate underflow is impossible.                      */
 133          break;
 134
 135       mpfr_log (w, w, GMP_RNDN);
 136          /* generic error of log: (2^(- exp(w)) + 0.5) ulp */
 137
 138       if (mpfr_zero_p (w))
 139          /* impossible to round, switch to second algorithm */
 140          break;
 141
 142       err = MPC_MAX (-mpfr_get_exp (w), 0) + 1;
 143          /* number of lost digits */
 144       ok = mpfr_can_round (w, prec - err, GMP_RNDN, GMP_RNDZ,
 145          mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN));
 146    }
 147
 148    if (!ok) {
 149       prec = MPC_PREC_RE(rop);
 150       mpfr_init2 (v, 2);
 151       /* compute 1/2 log (x^2 + y^2) = log |x| + 1/2 * log (1 + (y/x)^2)
 152             if |x| >= |y|; otherwise, exchange x and y                   */
 153       if (mpfr_cmpabs (mpc_realref (op), mpc_imagref (op)) >= 0) {
 154          x = mpc_realref (op);
 155          y = mpc_imagref (op);
 156       }
 157       else {
 158          x = mpc_imagref (op);
 159          y = mpc_realref (op);
 160       }
 161
 162       do {
 163          prec += mpc_ceil_log2 (prec) + 4;
 164          mpfr_set_prec (v, prec);
 165          mpfr_set_prec (w, prec);
 166
 167          mpfr_div (v, y, x, GMP_RNDD); /* error 1 ulp */
 168          mpfr_sqr (v, v, GMP_RNDD);
 169             /* generic error of multiplication:
 170                1 + 2*1*(2+1*2^(1-prec)) <= 5.0625 since prec >= 6 */
 171          mpfr_log1p (v, v, GMP_RNDD);
 172             /* error 1 + 4*5.0625 = 21.25 , see algorithms.tex */
 173          mpfr_div_2ui (v, v, 1, GMP_RNDD);
 174             /* If the result is 0, then there has been an underflow somewhere. */
 175
 176          mpfr_abs (w, x, GMP_RNDN); /* exact */
 177          mpfr_log (w, w, GMP_RNDN); /* error 0.5 ulp */
 178          expw = mpfr_get_exp (w);
 179          sgnw = mpfr_signbit (w);
 180
 181          mpfr_add (w, w, v, GMP_RNDN);
 182          if (!sgnw) /* v is positive, so no cancellation;
 183                        error 22.25 ulp; error counts lost bits */
 184             err = 5;
 185          else
 186             err =   MPC_MAX (5 + mpfr_get_exp (v),
 187                   /* 21.25 ulp (v) rewritten in ulp (result, now in w) */
 188                            -1 + expw             - mpfr_get_exp (w)
 189                   /* 0.5 ulp (previous w), rewritten in ulp (result) */
 190                   ) + 2;
 191
 192          /* handle one special case: |x|=1, and (y/x)^2 underflows;
 193             then 1/2*log(x^2+y^2) \approx 1/2*y^2 also underflows.  */
 194          if (   (mpfr_cmp_si (x, -1) == 0 || mpfr_cmp_ui (x, 1) == 0)
 195              && mpfr_zero_p (w))
 196             underflow = 1;
 197
 198       } while (!underflow &&
 199                !mpfr_can_round (w, prec - err, GMP_RNDN, GMP_RNDZ,
 200                mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN)));
 201       mpfr_clear (v);
 202    }
 203
 204    /* imaginary part */
 205    inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op),
 206                          MPC_RND_IM (rnd));
 207
 208    /* set the real part; cannot be done before if rop==op */
 209    if (underflow)
 210       /* create underflow in result */
 211       inex_re = mpfr_set_ui_2exp (mpc_realref (rop), 1,
 212                                   mpfr_get_emin_min () - 2, MPC_RND_RE (rnd));
 213    else
 214       inex_re = mpfr_set (mpc_realref (rop), w, MPC_RND_RE (rnd));
 215    mpfr_clear (w);
 216    return MPC_INEX(inex_re, inex_im);
 217 }