src/log.c

   1 /* mpc_log -- Take the logarithm of a complex number.
   2
   3 Copyright (C) 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 int
  25 mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){
  26    int ok=0;
  27    mpfr_t w;
  28    mp_prec_t prec;
  29    int loops = 0;
  30    int re_cmp, im_cmp;
  31    int inex_re, inex_im;
  32
  33    /* special values: NaN and infinities */
  34    if (!mpfr_number_p (MPC_RE (op)) || !mpfr_number_p (MPC_IM (op))) {
  35       if (mpfr_nan_p (MPC_RE (op))) {
  36          if (mpfr_inf_p (MPC_IM (op)))
  37             mpfr_set_inf (MPC_RE (rop), +1);
  38          else
  39             mpfr_set_nan (MPC_RE (rop));
  40          mpfr_set_nan (MPC_IM (rop));
  41          inex_im = 0; /* Inf/NaN is exact */
  42       }
  43       else if (mpfr_nan_p (MPC_IM (op))) {
  44          if (mpfr_inf_p (MPC_RE (op)))
  45             mpfr_set_inf (MPC_RE (rop), +1);
  46          else
  47             mpfr_set_nan (MPC_RE (rop));
  48          mpfr_set_nan (MPC_IM (rop));
  49          inex_im = 0; /* Inf/NaN is exact */
  50       }
  51       else /* We have an infinity in at least one part. */ {
  52          inex_im = mpfr_atan2 (MPC_IM (rop), MPC_IM (op), MPC_RE (op),
  53                                MPC_RND_IM (rnd));
  54          mpfr_set_inf (MPC_RE (rop), +1);
  55       }
  56       return MPC_INEX(0, inex_im);
  57    }
  58
  59    /* special cases: real and purely imaginary numbers */
  60    re_cmp = mpfr_cmp_ui (MPC_RE (op), 0);
  61    im_cmp = mpfr_cmp_ui (MPC_IM (op), 0);
  62    if (im_cmp == 0) {
  63       if (re_cmp == 0) {
  64          inex_im = mpfr_atan2 (MPC_IM (rop), MPC_IM (op), MPC_RE (op),
  65                                MPC_RND_IM (rnd));
  66          mpfr_set_inf (MPC_RE (rop), -1);
  67          inex_re = 0; /* -Inf is exact */
  68       }
  69       else if (re_cmp > 0) {
  70          inex_re = mpfr_log (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
  71          inex_im = mpfr_set (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd));
  72       }
  73       else {
  74          /* op = x + 0*y; let w = -x = |x| */
  75          int negative_zero;
  76          mpfr_rnd_t rnd_im;
  77
  78          negative_zero = mpfr_signbit (MPC_IM (op));
  79          if (negative_zero)
  80             rnd_im = INV_RND (MPC_RND_IM (rnd));
  81          else
  82             rnd_im = MPC_RND_IM (rnd);
  83          w [0] = *MPC_RE (op);
  84          MPFR_CHANGE_SIGN (w);
  85          inex_re = mpfr_log (MPC_RE (rop), w, MPC_RND_RE (rnd));
  86          inex_im = mpfr_const_pi (MPC_IM (rop), rnd_im);
  87          if (negative_zero) {
  88             mpc_conj (rop, rop, MPC_RNDNN);
  89             inex_im = -inex_im;
  90          }
  91       }
  92       return MPC_INEX(inex_re, inex_im);
  93    }
  94    else if (re_cmp == 0) {
  95       if (im_cmp > 0) {
  96          inex_re = mpfr_log (MPC_RE (rop), MPC_IM (op), MPC_RND_RE (rnd));
  97          inex_im = mpfr_const_pi (MPC_IM (rop), MPC_RND_IM (rnd));
  98          /* division by 2 does not change the ternary flag */
  99          mpfr_div_2ui (MPC_IM (rop), MPC_IM (rop), 1, GMP_RNDN);
 100       }
 101       else {
 102          w [0] = *MPC_IM (op);
 103          MPFR_CHANGE_SIGN (w);
 104          inex_re = mpfr_log (MPC_RE (rop), w, MPC_RND_RE (rnd));
 105          inex_im = mpfr_const_pi (MPC_IM (rop), INV_RND (MPC_RND_IM (rnd)));
 106          /* division by 2 does not change the ternary flag */
 107          mpfr_div_2ui (MPC_IM (rop), MPC_IM (rop), 1, GMP_RNDN);
 108          mpfr_neg (MPC_IM (rop), MPC_IM (rop), GMP_RNDN);
 109          inex_im = -inex_im; /* negate the ternary flag */
 110       }
 111       return MPC_INEX(inex_re, inex_im);
 112    }
 113
 114    prec = MPC_PREC_RE(rop);
 115    mpfr_init2 (w, prec);
 116    /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */
 117    /* loop for the real part: log (x^2 + y^2)                    */
 118    do {
 119       loops ++;
 120       prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
 121       mpfr_set_prec (w, prec);
 122
 123       /* w is rounded down */
 124       mpc_norm (w, op, GMP_RNDD);
 125       /* error 1 ulp */
 126
 127       if (mpfr_inf_p (w))
 128          /* FIXME
 129             return +inf, which is wrong since the logarithm is representable */
 130          ok = 1;
 131       else {
 132          mpfr_log (w, w, GMP_RNDD);
 133          /* generic error of log: (2^(2 - exp(w)) + 1) ulp */
 134
 135          if (MPFR_EXP (w) >= 2)
 136             ok = mpfr_can_round (w, prec - 2, GMP_RNDD, MPC_RND_RE(rnd), MPC_PREC_RE(rop));
 137          else
 138             ok = mpfr_can_round (w, prec - 3 + MPFR_EXP (w), GMP_RNDD, MPC_RND_RE(rnd), MPC_PREC_RE(rop));
 139       }
 140    } while (ok == 0);
 141
 142    /* imaginary part */
 143    inex_im = mpfr_atan2 (MPC_IM (rop), MPC_IM (op), MPC_RE (op),
 144                          MPC_RND_IM (rnd));
 145
 146    /* set the real part; cannot be done before when rop==op */
 147    inex_re = mpfr_div_2ui (MPC_RE(rop), w, 1ul, MPC_RND_RE (rnd));
 148    mpfr_clear (w);
 149    return MPC_INEX(inex_re, inex_im);
 150 }