1 /* mpfr_log -- natural logarithm of a floating-point number 2 3 Copyright 1999-2018 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 http://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_exp_t m; 125 mpfr_exp_t cancel; 126 127 /* Calculus of m (depends on p) 128 If mpfr_exp_t has N bits, then both (p + 3) / 2 and |exp_a| fit 129 on N-2 bits, so that there cannot be an overflow. */ 130 m = (p + 3) / 2 - exp_a; 131 132 /* In standard configuration (_MPFR_EXP_FORMAT <= 3), one has 133 mpfr_exp_t <= long, so that the following assertion is always 134 true. */ 135 MPFR_ASSERTN (m >= LONG_MIN && m <= LONG_MAX); 136 137 /* FIXME: Why 1 ulp and not 1/2 ulp? Ditto with some other ones 138 below. The error concerning the AGM should be explained since 139 4/s is inexact (one needs a bound on its derivative). */ 140 mpfr_mul_2si (tmp2, a, m, MPFR_RNDN); /* s=a*2^m, err<=1 ulp */ 141 MPFR_ASSERTD (MPFR_EXP (tmp2) >= (p + 3) / 2); 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, tmp2, MPFR_RNDN);/* 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