1 /* mpc_log10 -- Take the base-10 logarithm of a complex number. 2 3 Copyright (C) 2012, 2020 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://logw.gnu.org/licenses/ . 19 */ 20 21 #include <limits.h> /* for CHAR_BIT */ 22 #include "mpc-impl.h" 23 24 static void 25 mpfr_const_log10 (mpfr_ptr log10) 26 { 27 mpfr_set_ui (log10, 10, MPFR_RNDN); /* exact since prec >= 4 */ 28 mpfr_log (log10, log10, MPFR_RNDN); /* error <= 1/2 ulp */ 29 } 30 31 32 int 33 mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 34 { 35 int ok = 0, loops = 0, check_exact = 0, special_re, special_im, 36 inex, inex_re, inex_im; 37 mpfr_prec_t prec; 38 mpfr_t log10; 39 mpc_t log; 40 mpfr_exp_t saved_emin, saved_emax; 41 42 saved_emin = mpfr_get_emin (); 43 saved_emax = mpfr_get_emax (); 44 mpfr_set_emin (mpfr_get_emin_min ()); 45 mpfr_set_emax (mpfr_get_emax_max ()); 46 47 mpfr_init2 (log10, 2); 48 mpc_init2 (log, 2); 49 prec = MPC_MAX_PREC (rop); 50 /* compute log(op)/log(10) */ 51 while (ok == 0) { 52 loops ++; 53 prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; 54 mpfr_set_prec (log10, prec); 55 mpc_set_prec (log, prec); 56 57 inex = mpc_log (log, op, rnd); /* error <= 1 ulp */ 58 59 if (!mpfr_number_p (mpc_imagref (log)) 60 || mpfr_zero_p (mpc_imagref (log))) { 61 /* no need to divide by log(10) */ 62 special_im = 1; 63 ok = 1; 64 } 65 else { 66 special_im = 0; 67 mpfr_const_log10 (log10); 68 mpfr_div (mpc_imagref (log), mpc_imagref (log), log10, MPFR_RNDN); 69 70 ok = mpfr_can_round (mpc_imagref (log), prec - 2, 71 MPFR_RNDN, MPFR_RNDZ, 72 MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN)); 73 } 74 75 if (ok) { 76 if (!mpfr_number_p (mpc_realref (log)) 77 || mpfr_zero_p (mpc_realref (log))) 78 special_re = 1; 79 else { 80 special_re = 0; 81 if (special_im) 82 /* log10 not yet computed */ 83 mpfr_const_log10 (log10); 84 mpfr_div (mpc_realref (log), mpc_realref (log), log10, MPFR_RNDN); 85 /* error <= 24/7 ulp < 4 ulp for prec >= 4, see algorithms.tex */ 86 87 ok = mpfr_can_round (mpc_realref (log), prec - 2, 88 MPFR_RNDN, MPFR_RNDZ, 89 MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == MPFR_RNDN)); 90 } 91 92 /* Special code to deal with cases where the real part of log10(x+i*y) 93 is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2 94 this happens whenever x^2+y^2 is a nonnegative power of 10. 95 Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0, 96 since x^2+y^2 is rational, and 10^(a/2^b) is irrational. 97 Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2 98 is a rational with denominator a power of 2. 99 Now let x^2+y^2 = 10^s. Without loss of generality we can assume 100 x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e) 101 thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily 102 u = v = 0 mod 2^e, thus x and y are necessarily integers. 103 */ 104 if (!ok && !check_exact && mpfr_integer_p (mpc_realref (op)) && 105 mpfr_integer_p (mpc_imagref (op))) { 106 mpz_t x, y; 107 unsigned long s, v; 108 109 check_exact = 1; 110 mpz_init (x); 111 mpz_init (y); 112 mpfr_get_z (x, mpc_realref (op), MPFR_RNDN); /* exact */ 113 mpfr_get_z (y, mpc_imagref (op), MPFR_RNDN); /* exact */ 114 mpz_mul (x, x, x); 115 mpz_mul (y, y, y); 116 mpz_add (x, x, y); /* x^2+y^2 */ 117 v = mpz_scan1 (x, 0); 118 /* if x = 10^s then necessarily s = v */ 119 s = mpz_sizeinbase (x, 10); 120 /* since s is either the number of digits of x or one more, 121 then x = 10^(s-1) or 10^(s-2) */ 122 if (s == v + 1 || s == v + 2) { 123 mpz_div_2exp (x, x, v); 124 mpz_ui_pow_ui (y, 5, v); 125 if (mpz_cmp (y, x) == 0) { 126 /* Re(log10(x+i*y)) is exactly v/2 127 we reset the precision of Re(log) so that v can be 128 represented exactly */ 129 mpfr_set_prec (mpc_realref (log), 130 sizeof(unsigned long)*CHAR_BIT); 131 mpfr_set_ui_2exp (mpc_realref (log), v, -1, MPFR_RNDN); 132 /* exact */ 133 ok = 1; 134 } 135 } 136 mpz_clear (x); 137 mpz_clear (y); 138 } 139 } 140 } 141 142 inex_re = mpfr_set (mpc_realref(rop), mpc_realref (log), MPC_RND_RE (rnd)); 143 if (special_re) 144 inex_re = MPC_INEX_RE (inex); 145 /* recover flag from call to mpc_log above */ 146 inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (log), MPC_RND_IM (rnd)); 147 if (special_im) 148 inex_im = MPC_INEX_IM (inex); 149 mpfr_clear (log10); 150 mpc_clear (log); 151 152 /* restore the exponent range, and check the range of results */ 153 mpfr_set_emin (saved_emin); 154 mpfr_set_emax (saved_emax); 155 inex_re = mpfr_check_range (mpc_realref (rop), inex_re, MPC_RND_RE (rnd)); 156 inex_im = mpfr_check_range (mpc_imagref (rop), inex_im, MPC_RND_IM (rnd)); 157 158 return MPC_INEX(inex_re, inex_im); 159 } 160