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 } 81 prec += prec / 2; 82 mpc_set_prec (tmp, prec); 83 mpfr_set_prec (log10, prec); 84 } 85 mpc_clear (tmp); 86 mpfr_clear (log10); 87 return ret; 88 } 89 90 int 91 mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 92 { 93 int ok = 0, loops = 0, re_cmp, im_cmp, inex_re, inex_im, negative_zero; 94 mpfr_t w; 95 mpfr_prec_t prec; 96 mpfr_rnd_t rnd_im; 97 mpc_t ww; 98 mpc_rnd_t invrnd; 99 100 /* special values: NaN and infinities: same as mpc_log */ 101 if (!mpc_fin_p (op)) /* real or imaginary parts are NaN or Inf */ 102 { 103 if (mpfr_nan_p (mpc_realref (op))) 104 { 105 if (mpfr_inf_p (mpc_imagref (op))) 106 /* (NaN, Inf) -> (+Inf, NaN) */ 107 mpfr_set_inf (mpc_realref (rop), +1); 108 else 109 /* (NaN, xxx) -> (NaN, NaN) */ 110 mpfr_set_nan (mpc_realref (rop)); 111 mpfr_set_nan (mpc_imagref (rop)); 112 inex_im = 0; /* Inf/NaN is exact */ 113 } 114 else if (mpfr_nan_p (mpc_imagref (op))) 115 { 116 if (mpfr_inf_p (mpc_realref (op))) 117 /* (Inf, NaN) -> (+Inf, NaN) */ 118 mpfr_set_inf (mpc_realref (rop), +1); 119 else 120 /* (xxx, NaN) -> (NaN, NaN) */ 121 mpfr_set_nan (mpc_realref (rop)); 122 mpfr_set_nan (mpc_imagref (rop)); 123 inex_im = 0; /* Inf/NaN is exact */ 124 } 125 else /* We have an infinity in at least one part. */ 126 { 127 /* (+Inf, y) -> (+Inf, 0) for finite positive-signed y */ 128 if (mpfr_inf_p (mpc_realref (op)) && mpfr_signbit (mpc_realref (op)) 129 == 0 && mpfr_number_p (mpc_imagref (op))) 130 inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), 131 mpc_realref (op), MPC_RND_IM (rnd)); 132 else 133 /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx)) 134 (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */ 135 inex_im = mpc_log10_aux (rop, op, rnd, 1, 0); 136 mpfr_set_inf (mpc_realref (rop), +1); 137 } 138 return MPC_INEX(0, inex_im); 139 } 140 141 /* special cases: real and purely imaginary numbers */ 142 re_cmp = mpfr_cmp_ui (mpc_realref (op), 0); 143 im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0); 144 if (im_cmp == 0) /* Im(op) = 0 */ 145 { 146 if (re_cmp == 0) /* Re(op) = 0 */ 147 { 148 if (mpfr_signbit (mpc_realref (op)) == 0) 149 inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), 150 mpc_realref (op), MPC_RND_IM (rnd)); 151 else 152 inex_im = mpc_log10_aux (rop, op, rnd, 1, 0); 153 mpfr_set_inf (mpc_realref (rop), -1); 154 inex_re = 0; /* -Inf is exact */ 155 } 156 else if (re_cmp > 0) 157 { 158 inex_re = mpfr_log10 (mpc_realref (rop), mpc_realref (op), 159 MPC_RND_RE (rnd)); 160 inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op), 161 MPC_RND_IM (rnd)); 162 } 163 else /* log10(x + 0*i) for negative x */ 164 { /* op = x + 0*i; let w = -x = |x| */ 165 negative_zero = mpfr_signbit (mpc_imagref (op)); 166 if (negative_zero) 167 rnd_im = INV_RND (MPC_RND_IM (rnd)); 168 else 169 rnd_im = MPC_RND_IM (rnd); 170 ww->re[0] = *mpc_realref (op); 171 MPFR_CHANGE_SIGN (ww->re); 172 ww->im[0] = *mpc_imagref (op); 173 if (mpfr_cmp_ui (ww->re, 1) == 0) 174 inex_re = mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); 175 else 176 inex_re = mpc_log10_aux (rop, ww, rnd, 0, 1); 177 inex_im = mpc_log10_aux (rop, op, MPC_RND (0,rnd_im), 1, 2); 178 if (negative_zero) 179 { 180 mpc_conj (rop, rop, MPC_RNDNN); 181 inex_im = -inex_im; 182 } 183 } 184 return MPC_INEX(inex_re, inex_im); 185 } 186 else if (re_cmp == 0) 187 { 188 if (im_cmp > 0) 189 { 190 inex_re = mpfr_log10 (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd)); 191 inex_im = mpc_log10_aux (rop, op, rnd, 1, 2); 192 /* division by 2 does not change the ternary flag */ 193 mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); 194 } 195 else 196 { 197 w [0] = *mpc_imagref (op); 198 MPFR_CHANGE_SIGN (w); 199 inex_re = mpfr_log10 (mpc_realref (rop), w, MPC_RND_RE (rnd)); 200 invrnd = MPC_RND (0, INV_RND (MPC_RND_IM (rnd))); 201 inex_im = mpc_log10_aux (rop, op, invrnd, 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 mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); 205 inex_im = -inex_im; /* negate the ternary flag */ 206 } 207 return MPC_INEX(inex_re, inex_im); 208 } 209 210 /* generic case: neither Re(op) nor Im(op) is NaN, Inf or zero */ 211 prec = MPC_PREC_RE(rop); 212 mpfr_init2 (w, prec); 213 mpc_init2 (ww, prec); 214 /* let op = x + iy; compute log(op)/log(10) */ 215 while (ok == 0) 216 { 217 loops ++; 218 prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; 219 mpfr_set_prec (w, prec); 220 mpc_set_prec (ww, prec); 221 222 mpc_log (ww, op, MPC_RNDNN); 223 mpfr_set_ui (w, 10, GMP_RNDN); /* exact since prec >= 4 */ 224 mpfr_log (w, w, GMP_RNDN); 225 mpc_div_fr (ww, ww, w, MPC_RNDNN); 226 227 ok = mpfr_can_round (mpc_realref (ww), prec - 2, GMP_RNDN, GMP_RNDZ, 228 MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == GMP_RNDN)); 229 230 /* Special code to deal with cases where the real part of log10(x+i*y) 231 is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2 232 this happens whenever x^2+y^2 is a nonnegative power of 10. 233 Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0, 234 since x^2+y^2 is rational, and 10^(a/2^b) is irrational. 235 Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2 236 is a rational with denominator a power of 2. 237 Now let x^2+y^2 = 10^s. Without loss of generality we can assume 238 x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e) 239 thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily 240 u = v = 0 mod 2^e, thus x and y are necessarily integers. 241 */ 242 if ((ok == 0) && (loops == 1) && mpfr_integer_p (mpc_realref (op)) && 243 mpfr_integer_p (mpc_imagref (op))) 244 { 245 mpz_t x, y; 246 unsigned long s, v; 247 248 mpz_init (x); 249 mpz_init (y); 250 mpfr_get_z (x, mpc_realref (op), GMP_RNDN); /* exact */ 251 mpfr_get_z (y, mpc_imagref (op), GMP_RNDN); /* exact */ 252 mpz_mul (x, x, x); 253 mpz_mul (y, y, y); 254 mpz_add (x, x, y); /* x^2+y^2 */ 255 v = mpz_scan1 (x, 0); 256 /* if x = 10^s then necessarily s = v */ 257 s = mpz_sizeinbase (x, 10); 258 /* since s is either the number of digits of x or one more, 259 then x = 10^(s-1) or 10^(s-2) */ 260 if (s == v + 1 || s == v + 2) 261 { 262 mpz_div_2exp (x, x, v); 263 mpz_ui_pow_ui (y, 5, v); 264 if (mpz_cmp (y, x) == 0) /* Re(log10(x+i*y)) is exactly v/2 */ 265 { 266 /* we reset the precision of Re(ww) so that v can be 267 represented exactly */ 268 mpfr_set_prec (mpc_realref (ww), sizeof(unsigned long)*CHAR_BIT); 269 mpfr_set_ui_2exp (mpc_realref (ww), v, -1, GMP_RNDN); /* exact */ 270 ok = 1; 271 } 272 } 273 mpz_clear (x); 274 mpz_clear (y); 275 } 276 277 ok = ok && mpfr_can_round (mpc_imagref (ww), prec-2, GMP_RNDN, GMP_RNDZ, 278 MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == GMP_RNDN)); 279 } 280 281 inex_re = mpfr_set (mpc_realref(rop), mpc_realref (ww), MPC_RND_RE (rnd)); 282 inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (ww), MPC_RND_IM (rnd)); 283 mpfr_clear (w); 284 mpc_clear (ww); 285 return MPC_INEX(inex_re, inex_im); 286 } 287