src/log10.c

   1 /* mpc_log10 -- Take the base-10 logarithm of a complex number.
   2
   3 Copyright (C) 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 <limits.h> /* for CHAR_BIT */
  22 #include "mpc-impl.h"
  23
  24 /* Auxiliary functions which implement Ziv's strategy for special cases.
  25    if flag = 0: compute only real part
  26    if flag = 1: compute only imaginary
  27    Exact cases should be dealt with separately. */
  28 static int
  29 mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb)
  30 {
  31   mp_prec_t prec = (MPFR_PREC_MIN > 4) ? MPFR_PREC_MIN : 4;
  32   mpc_t tmp;
  33   mpfr_t log10;
  34   int ok = 0, ret;
  35
  36   prec = mpfr_get_prec ((flag == 0) ? mpc_realref (rop) : mpc_imagref (rop));
  37   prec += 10;
  38   mpc_init2 (tmp, prec);
  39   mpfr_init2 (log10, prec);
  40   while (ok == 0)
  41     {
  42       mpfr_set_ui (log10, 10, GMP_RNDN); /* exact since prec >= 4 */
  43       mpfr_log (log10, log10, GMP_RNDN);
  44       /* In each case we have two roundings, thus the final value is
  45          x * (1+u)^2 where x is the exact value, and |u| <= 2^(-prec-1).
  46          Thus the error is always less than 3 ulps. */
  47       switch (nb)
  48         {
  49         case 0: /* imag <- atan2(y/x) */
  50           mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op),
  51                       MPC_RND_IM (rnd));
  52           mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN);
  53           ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDN,
  54                                GMP_RNDZ, MPC_PREC_IM(rop) +
  55                                (MPC_RND_IM (rnd) == GMP_RNDN));
  56           if (ok)
  57             ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp),
  58                             MPC_RND_IM (rnd));
  59           break;
  60         case 1: /* real <- log(x) */
  61           mpfr_log (mpc_realref (tmp), mpc_realref (op), MPC_RND_RE (rnd));
  62           mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN);
  63           ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDN,
  64                                GMP_RNDZ, MPC_PREC_RE(rop) +
  65                                (MPC_RND_RE (rnd) == GMP_RNDN));
  66           if (ok)
  67             ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp),
  68                             MPC_RND_RE (rnd));
  69           break;
  70         case 2: /* imag <- pi */
  71           mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd));
  72           mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN);
  73           ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDN,
  74                                GMP_RNDZ, MPC_PREC_IM(rop) +
  75                                (MPC_RND_IM (rnd) == GMP_RNDN));
  76           if (ok)
  77             ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp),
  78                             MPC_RND_IM (rnd));
  79           break;
  80         case 3: /* real <- log(y) */
  81           mpfr_log (mpc_realref (tmp), mpc_imagref (op), MPC_RND_RE (rnd));
  82           mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN);
  83           ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDN,
  84                                GMP_RNDZ, MPC_PREC_RE(rop) +
  85                                (MPC_RND_RE (rnd) == GMP_RNDN));
  86           if (ok)
  87             ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp),
  88                             MPC_RND_RE (rnd));
  89           break;
  90         }
  91       prec += prec / 2;
  92       mpc_set_prec (tmp, prec);
  93       mpfr_set_prec (log10, prec);
  94     }
  95   mpc_clear (tmp);
  96   mpfr_clear (log10);
  97   return ret;
  98 }
  99
 100 int
 101 mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 102 {
 103   int ok = 0, loops = 0, re_cmp, im_cmp, inex_re, inex_im, negative_zero;
 104   mpfr_t w;
 105   mpfr_prec_t prec;
 106   mpfr_rnd_t rnd_im;
 107   mpc_t ww;
 108   mpc_rnd_t invrnd;
 109
 110   /* special values: NaN and infinities: same as mpc_log */
 111   if (!mpc_fin_p (op)) /* real or imaginary parts are NaN or Inf */
 112     {
 113       if (mpfr_nan_p (mpc_realref (op)))
 114         {
 115           if (mpfr_inf_p (mpc_imagref (op)))
 116             /* (NaN, Inf) -> (+Inf, NaN) */
 117             mpfr_set_inf (mpc_realref (rop), +1);
 118           else
 119             /* (NaN, xxx) -> (NaN, NaN) */
 120             mpfr_set_nan (mpc_realref (rop));
 121          mpfr_set_nan (mpc_imagref (rop));
 122          inex_im = 0; /* Inf/NaN is exact */
 123         }
 124       else if (mpfr_nan_p (mpc_imagref (op)))
 125         {
 126           if (mpfr_inf_p (mpc_realref (op)))
 127             /* (Inf, NaN) -> (+Inf, NaN) */
 128             mpfr_set_inf (mpc_realref (rop), +1);
 129           else
 130             /* (xxx, NaN) -> (NaN, NaN) */
 131             mpfr_set_nan (mpc_realref (rop));
 132           mpfr_set_nan (mpc_imagref (rop));
 133           inex_im = 0; /* Inf/NaN is exact */
 134         }
 135       else /* We have an infinity in at least one part. */
 136         {
 137           /* (+Inf, y) -> (+Inf, 0) for finite positive-signed y */
 138           if (mpfr_inf_p (mpc_realref (op)) && mpfr_signbit (mpc_realref (op))
 139               == 0 && mpfr_number_p (mpc_imagref (op)))
 140             inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op),
 141                                   mpc_realref (op), MPC_RND_IM (rnd));
 142           else
 143             /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx))
 144                (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */
 145             inex_im = mpc_log10_aux (rop, op, rnd, 1, 0);
 146           mpfr_set_inf (mpc_realref (rop), +1);
 147         }
 148       return MPC_INEX(0, inex_im);
 149     }
 150
 151    /* special cases: real and purely imaginary numbers */
 152   re_cmp = mpfr_cmp_ui (mpc_realref (op), 0);
 153   im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0);
 154   if (im_cmp == 0) /* Im(op) = 0 */
 155     {
 156       if (re_cmp == 0) /* Re(op) = 0 */
 157         {
 158           if (mpfr_signbit (mpc_realref (op)) == 0)
 159             inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op),
 160                                   mpc_realref (op), MPC_RND_IM (rnd));
 161           else
 162             inex_im = mpc_log10_aux (rop, op, rnd, 1, 0);
 163           mpfr_set_inf (mpc_realref (rop), -1);
 164           inex_re = 0; /* -Inf is exact */
 165         }
 166       else if (re_cmp > 0)
 167         {
 168           inex_re = mpfr_log10 (mpc_realref (rop), mpc_realref (op),
 169                                 MPC_RND_RE (rnd));
 170           inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op),
 171                               MPC_RND_IM (rnd));
 172         }
 173       else /* log10(x + 0*i) for negative x */
 174         { /* op = x + 0*i; let w = -x = |x| */
 175           negative_zero = mpfr_signbit (mpc_imagref (op));
 176           if (negative_zero)
 177             rnd_im = INV_RND (MPC_RND_IM (rnd));
 178           else
 179             rnd_im = MPC_RND_IM (rnd);
 180           ww->re[0] = *mpc_realref (op);
 181           MPFR_CHANGE_SIGN (ww->re);
 182           ww->im[0] = *mpc_imagref (op);
 183           if (mpfr_cmp_ui (ww->re, 1) == 0)
 184             inex_re = mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
 185           else
 186             inex_re = mpc_log10_aux (rop, ww, rnd, 0, 1);
 187           inex_im = mpc_log10_aux (rop, op, MPC_RND (0,rnd_im), 1, 2);
 188           if (negative_zero)
 189             {
 190               mpc_conj (rop, rop, MPC_RNDNN);
 191               inex_im = -inex_im;
 192             }
 193         }
 194       return MPC_INEX(inex_re, inex_im);
 195     }
 196    else if (re_cmp == 0)
 197      {
 198        if (im_cmp > 0)
 199          {
 200            inex_re = mpc_log10_aux (rop, op, rnd, 0, 3);
 201            inex_im = mpc_log10_aux (rop, op, rnd, 1, 2);
 202            /* division by 2 does not change the ternary flag */
 203            mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
 204          }
 205        else
 206          {
 207            ww->re[0] = *mpc_realref (op);
 208            ww->im[0] = *mpc_imagref (op);
 209            MPFR_CHANGE_SIGN (ww->im);
 210            inex_re = mpc_log10_aux (rop, ww, rnd, 0, 3);
 211            invrnd = MPC_RND (0, INV_RND (MPC_RND_IM (rnd)));
 212            inex_im = mpc_log10_aux (rop, op, invrnd, 1, 2);
 213            /* division by 2 does not change the ternary flag */
 214            mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN);
 215            mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN);
 216            inex_im = -inex_im; /* negate the ternary flag */
 217          }
 218        return MPC_INEX(inex_re, inex_im);
 219      }
 220
 221   /* generic case: neither Re(op) nor Im(op) is NaN, Inf or zero */
 222   prec = MPC_PREC_RE(rop);
 223   mpfr_init2 (w, prec);
 224   mpc_init2 (ww, prec);
 225   /* let op = x + iy; compute log(op)/log(10) */
 226   while (ok == 0)
 227     {
 228       loops ++;
 229       prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
 230       mpfr_set_prec (w, prec);
 231       mpc_set_prec (ww, prec);
 232
 233       mpc_log (ww, op, MPC_RNDNN);
 234       mpfr_set_ui (w, 10, GMP_RNDN); /* exact since prec >= 4 */
 235       mpfr_log (w, w, GMP_RNDN);
 236       mpc_div_fr (ww, ww, w, MPC_RNDNN);
 237
 238       ok = mpfr_can_round (mpc_realref (ww), prec - 2, GMP_RNDN, GMP_RNDZ,
 239                            MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == GMP_RNDN));
 240
 241       /* Special code to deal with cases where the real part of log10(x+i*y)
 242          is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2
 243          this happens whenever x^2+y^2 is a nonnegative power of 10.
 244          Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0,
 245          since x^2+y^2 is rational, and 10^(a/2^b) is irrational.
 246          Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2
 247          is a rational with denominator a power of 2.
 248          Now let x^2+y^2 = 10^s. Without loss of generality we can assume
 249          x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e)
 250          thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily
 251          u = v = 0 mod 2^e, thus x and y are necessarily integers.
 252       */
 253       if ((ok == 0) && (loops == 1) && mpfr_integer_p (mpc_realref (op)) &&
 254           mpfr_integer_p (mpc_imagref (op)))
 255         {
 256           mpz_t x, y;
 257           unsigned long s, v;
 258
 259           mpz_init (x);
 260           mpz_init (y);
 261           mpfr_get_z (x, mpc_realref (op), GMP_RNDN); /* exact */
 262           mpfr_get_z (y, mpc_imagref (op), GMP_RNDN); /* exact */
 263           mpz_mul (x, x, x);
 264           mpz_mul (y, y, y);
 265           mpz_add (x, x, y); /* x^2+y^2 */
 266           v = mpz_scan1 (x, 0);
 267           /* if x = 10^s then necessarily s = v */
 268           s = mpz_sizeinbase (x, 10);
 269           /* since s is either the number of digits of x or one more,
 270              then x = 10^(s-1) or 10^(s-2) */
 271           if (s == v + 1 || s == v + 2)
 272             {
 273               mpz_div_2exp (x, x, v);
 274               mpz_ui_pow_ui (y, 5, v);
 275               if (mpz_cmp (y, x) == 0) /* Re(log10(x+i*y)) is exactly v/2 */
 276                 {
 277                   /* we reset the precision of Re(ww) so that v can be
 278                      represented exactly */
 279                   mpfr_set_prec (mpc_realref (ww), sizeof(unsigned long)*CHAR_BIT);
 280                   mpfr_set_ui_2exp (mpc_realref (ww), v, -1, GMP_RNDN); /* exact */
 281                   ok = 1;
 282                 }
 283             }
 284           mpz_clear (x);
 285           mpz_clear (y);
 286         }
 287
 288       ok = ok && mpfr_can_round (mpc_imagref (ww), prec-2, GMP_RNDN, GMP_RNDZ,
 289                             MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == GMP_RNDN));
 290     }
 291
 292   inex_re = mpfr_set (mpc_realref(rop), mpc_realref (ww), MPC_RND_RE (rnd));
 293   inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (ww), MPC_RND_IM (rnd));
 294   mpfr_clear (w);
 295   mpc_clear (ww);
 296   return MPC_INEX(inex_re, inex_im);
 297 }