1/* mpc_log10 -- Take the base-10 logarithm of a complex number. 2 3Copyright (C) 2012 INRIA 4 5This file is part of GNU MPC. 6 7GNU MPC is free software; you can redistribute it and/or modify it under 8the terms of the GNU Lesser General Public License as published by the 9Free Software Foundation; either version 3 of the License, or (at your 10option) any later version. 11 12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along 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. */ 28static int 29mpc_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 100int 101mpc_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} 298