1 /* mpfr_log -- natural logarithm of a floating-point number 2 3 Copyright 1999-2023 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 #define MPFR_NEED_LONGLONG_H 24 #include "mpfr-impl.h" 25 26 /* The computation of log(x) is done using the formula : 27 if we want p bits of the result, 28 29 pi 30 log(x) ~ ------------ - m log 2 31 2 AG(1,4/s) 32 33 where s = x 2^m > 2^(p/2) 34 35 More precisely, if F(x) = int(1/sqrt(1-(1-x^2)*sin(t)^2), t=0..PI/2), 36 then for s>=1.26 we have log(s) < F(4/s) < log(s)*(1+4/s^2) 37 from which we deduce pi/2/AG(1,4/s)*(1-4/s^2) < log(s) < pi/2/AG(1,4/s) 38 so the relative error 4/s^2 is < 4/2^p i.e. 4 ulps. 39 */ 40 41 int 42 mpfr_log (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode) 43 { 44 int inexact; 45 mpfr_prec_t p, q; 46 mpfr_t tmp1, tmp2; 47 mpfr_exp_t exp_a; 48 MPFR_SAVE_EXPO_DECL (expo); 49 MPFR_ZIV_DECL (loop); 50 MPFR_GROUP_DECL(group); 51 52 MPFR_LOG_FUNC 53 (("a[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode), 54 ("r[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r, 55 inexact)); 56 57 /* Special cases */ 58 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a))) 59 { 60 /* If a is NaN, the result is NaN */ 61 if (MPFR_IS_NAN (a)) 62 { 63 MPFR_SET_NAN (r); 64 MPFR_RET_NAN; 65 } 66 /* check for infinity before zero */ 67 else if (MPFR_IS_INF (a)) 68 { 69 if (MPFR_IS_NEG (a)) 70 /* log(-Inf) = NaN */ 71 { 72 MPFR_SET_NAN (r); 73 MPFR_RET_NAN; 74 } 75 else /* log(+Inf) = +Inf */ 76 { 77 MPFR_SET_INF (r); 78 MPFR_SET_POS (r); 79 MPFR_RET (0); 80 } 81 } 82 else /* a is zero */ 83 { 84 MPFR_ASSERTD (MPFR_IS_ZERO (a)); 85 MPFR_SET_INF (r); 86 MPFR_SET_NEG (r); 87 MPFR_SET_DIVBY0 (); 88 MPFR_RET (0); /* log(0) is an exact -infinity */ 89 } 90 } 91 92 /* If a is negative, the result is NaN */ 93 if (MPFR_UNLIKELY (MPFR_IS_NEG (a))) 94 { 95 MPFR_SET_NAN (r); 96 MPFR_RET_NAN; 97 } 98 99 exp_a = MPFR_GET_EXP (a); 100 101 /* If a is 1, the result is +0 */ 102 if (MPFR_UNLIKELY (exp_a == 1 && mpfr_cmp_ui (a, 1) == 0)) 103 { 104 MPFR_SET_ZERO (r); 105 MPFR_SET_POS (r); 106 MPFR_RET (0); /* only "normal" case where the result is exact */ 107 } 108 109 q = MPFR_PREC (r); 110 111 /* use initial precision about q+2*lg(q)+cte */ 112 p = q + 2 * MPFR_INT_CEIL_LOG2 (q) + 10; 113 /* % ~(mpfr_prec_t)GMP_NUMB_BITS ; 114 m=q; while (m) { p++; m >>= 1; } */ 115 /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0)) 116 p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */ 117 118 MPFR_SAVE_EXPO_MARK (expo); 119 MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2); 120 121 MPFR_ZIV_INIT (loop, p); 122 for (;;) 123 { 124 mpfr_t scaled_a; 125 mpfr_exp_t m; 126 mpfr_exp_t cancel; 127 128 /* Calculus of m (depends on p) 129 If mpfr_exp_t has N bits, then both (p + 3) / 2 and |exp_a| fit 130 on N-2 bits, so that there cannot be an overflow. */ 131 m = (p + 3) / 2 - exp_a; 132 133 /* In standard configuration (_MPFR_EXP_FORMAT <= 3), one has 134 mpfr_exp_t <= long, so that the following assertion is always 135 true. This assertion is needed for the mpfr_mul_si below. */ 136 MPFR_ASSERTN (m >= LONG_MIN && m <= LONG_MAX); 137 138 /* FIXME: Redo the error analysis. The error concerning the AGM 139 should be explained since 4/s is inexact (one needs a bound 140 on its derivative). */ 141 MPFR_ALIAS (scaled_a, a, MPFR_SIGN_POS, (p + 3) / 2); /* s=a*2^m */ 142 /* [FIXME] and one can have the equality, even if p is even. 143 This means that if a is a power of 2 and p is even, then 144 s = (1/2) * 2^((p+2)/2) = 2^(p/2), so that the condition 145 s > 2^(p/2) from algorithms.tex is not satisfied. */ 146 mpfr_div (tmp1, __gmpfr_four, scaled_a, MPFR_RNDF); /* 4/s, err<=2 ulps */ 147 mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */ 148 mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s), err<=3 ulps */ 149 mpfr_const_pi (tmp1, MPFR_RNDN); /* compute pi, err<=1ulp */ 150 mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */ 151 mpfr_const_log2 (tmp1, MPFR_RNDN); /* compute log(2), err<=1ulp */ 152 mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps */ 153 mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a), err<=7ulps+cancel */ 154 155 if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2))) 156 { 157 cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1); 158 MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel)); 159 MPFR_LOG_VAR (tmp1); 160 if (MPFR_UNLIKELY (cancel < 0)) 161 cancel = 0; 162 163 /* we have 7 ulps of error from the above roundings, 164 4 ulps from the 4/s^2 second order term, 165 plus the canceled bits */ 166 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p - cancel - 4, q, rnd_mode))) 167 break; 168 169 /* VL: I think it is better to have an increment that it isn't 170 too low; in particular, the increment must be positive even 171 if cancel = 0 (can this occur?). */ 172 p += cancel + MPFR_INT_CEIL_LOG2 (p); 173 } 174 else 175 { 176 /* TODO: find why this case can occur and what is best to do 177 with it. */ 178 p += MPFR_INT_CEIL_LOG2 (p); 179 } 180 181 MPFR_ZIV_NEXT (loop, p); 182 MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2); 183 } 184 MPFR_ZIV_FREE (loop); 185 inexact = mpfr_set (r, tmp1, rnd_mode); 186 /* We clean */ 187 MPFR_GROUP_CLEAR (group); 188 189 MPFR_SAVE_EXPO_FREE (expo); 190 return mpfr_check_range (r, inexact, rnd_mode); 191 } 192