1/* mpfr_log1p -- Compute log(1+x) 2 3Copyright 2001-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26/* Put in y an approximation of log(1+x) for x small. 27 We assume |x| < 1/2, in which case: 28 |x/2| <= |log(1+x)| <= |2x|. 29 Return k such that the error is bounded by 2^k*ulp(y). 30*/ 31static int 32mpfr_log1p_small (mpfr_ptr y, mpfr_srcptr x) 33{ 34 mpfr_prec_t p = MPFR_PREC(y), err; 35 mpfr_t t, u; 36 unsigned long i; 37 int k; 38 39 MPFR_ASSERTD(MPFR_GET_EXP (x) <= -1); /* ensures |x| < 1/2 */ 40 41 /* in the following, theta represents a value with |theta| <= 2^(1-p) 42 (might be a different value each time) */ 43 44 mpfr_init2 (t, p); 45 mpfr_init2 (u, p); 46 mpfr_set (t, x, MPFR_RNDF); /* t = x * (1 + theta) */ 47 mpfr_set (y, t, MPFR_RNDF); /* exact */ 48 for (i = 2; ; i++) 49 { 50 mpfr_mul (t, t, x, MPFR_RNDF); /* t = x^i * (1 + theta)^i */ 51 mpfr_div_ui (u, t, i, MPFR_RNDF); /* u = x^i/i * (1 + theta)^(i+1) */ 52 if (MPFR_GET_EXP (u) <= MPFR_GET_EXP (y) - p) /* |u| < ulp(y) */ 53 break; 54 if (i & 1) 55 mpfr_add (y, y, u, MPFR_RNDF); /* error <= ulp(y) */ 56 else 57 mpfr_sub (y, y, u, MPFR_RNDF); /* error <= ulp(y) */ 58 } 59 /* We assume |(1 + theta)^(i+1)| <= 2. 60 The neglected part is at most |u| + |u|/2 + ... <= 2|u| < 2 ulp(y) 61 which has to be multiplied by |(1 + theta)^(i+1)| <= 2, thus at most 62 4 ulp(y). 63 The rounding error on y is bounded by: 64 * for the (i-2) add/sub, each error is bounded by ulp(y), 65 and since |y| <= |x|, this yields (i-2)*ulp(x) 66 * from Lemma 3.1 from [Higham02] (see algorithms.tex), 67 the relative error on u at step i is bounded by: 68 (i+1)*epsilon/(1-(i+1)*epsilon) where epsilon = 2^(1-p). 69 If (i+1)*epsilon <= 1/2, then the relative error on u at 70 step i is bounded by 2*(i+1)*epsilon, and since |u| <= 1/2^(i+1) 71 at step i, this gives an absolute error bound of; 72 2*epsilon*x*(3/2^3 + 4/2^4 + 5/2^5 + ...) <= 2*2^(1-p)*x = 73 4*2^(-p)*x <= 4*ulp(x). 74 75 If (i+1)*epsilon <= 1/2, then the relative error on u at step i 76 is bounded by (i+1)*epsilon/(1-(i+1)*epsilon) <= 1, thus it follows 77 |(1 + theta)^(i+1)| <= 2. 78 79 Finally the total error is bounded by 4*ulp(y) + (i-2)*ulp(x) + 4*ulp(x) 80 = 4*ulp(y) + (i+2)*ulp(x). 81 Since x/2 <= y, we have ulp(x) <= 2*ulp(y), thus the error is bounded by: 82 (2*i+8)*ulp(y). 83 */ 84 err = 2 * i + 8; 85 k = __gmpfr_int_ceil_log2 (err); 86 MPFR_ASSERTN(k < p); 87 /* if k < p, since k = ceil(log2(err)), we have err <= 2^k <= 2^(p-1), 88 thus i+4 = err/2 <= 2^(p-2), thus (i+4)*epsilon <= 1/2, which implies 89 our assumption (i+1)*epsilon <= 1/2. */ 90 mpfr_clear (t); 91 mpfr_clear (u); 92 return k; 93} 94 95/* The computation of log1p is done by 96 log1p(x) = log(1+x) 97 except when x is very small, in which case log1p(x) = x + tiny error, 98 or when x is small, where we use directly the Taylor expansion. 99*/ 100 101int 102mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 103{ 104 int comp, inexact; 105 mpfr_exp_t ex; 106 MPFR_SAVE_EXPO_DECL (expo); 107 108 MPFR_LOG_FUNC 109 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 110 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 111 inexact)); 112 113 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 114 { 115 if (MPFR_IS_NAN (x)) 116 { 117 MPFR_SET_NAN (y); 118 MPFR_RET_NAN; 119 } 120 /* check for inf or -inf (result is not defined) */ 121 else if (MPFR_IS_INF (x)) 122 { 123 if (MPFR_IS_POS (x)) 124 { 125 MPFR_SET_INF (y); 126 MPFR_SET_POS (y); 127 MPFR_RET (0); 128 } 129 else 130 { 131 MPFR_SET_NAN (y); 132 MPFR_RET_NAN; 133 } 134 } 135 else /* x is zero */ 136 { 137 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 138 MPFR_SET_ZERO (y); /* log1p(+/- 0) = +/- 0 */ 139 MPFR_SET_SAME_SIGN (y, x); 140 MPFR_RET (0); 141 } 142 } 143 144 ex = MPFR_GET_EXP (x); 145 if (ex < 0) /* -0.5 < x < 0.5 */ 146 { 147 /* For x > 0, abs(log(1+x)-x) < x^2/2. 148 For x > -0.5, abs(log(1+x)-x) < x^2. */ 149 if (MPFR_IS_POS (x)) 150 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, 0, rnd_mode, {}); 151 else 152 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {}); 153 } 154 155 comp = mpfr_cmp_si (x, -1); 156 /* log1p(x) is undefined for x < -1 */ 157 if (MPFR_UNLIKELY(comp <= 0)) 158 { 159 if (comp == 0) 160 /* x=0: log1p(-1)=-inf (divide-by-zero exception) */ 161 { 162 MPFR_SET_INF (y); 163 MPFR_SET_NEG (y); 164 MPFR_SET_DIVBY0 (); 165 MPFR_RET (0); 166 } 167 MPFR_SET_NAN (y); 168 MPFR_RET_NAN; 169 } 170 171 MPFR_SAVE_EXPO_MARK (expo); 172 173 /* General case */ 174 { 175 /* Declaration of the intermediary variable */ 176 mpfr_t t; 177 /* Declaration of the size variable */ 178 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */ 179 mpfr_prec_t Nt; /* working precision */ 180 mpfr_exp_t err; /* error */ 181 MPFR_ZIV_DECL (loop); 182 183 /* compute the precision of intermediary variable */ 184 /* the optimal number of bits : see algorithms.tex */ 185 Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6; 186 187 /* if |x| is smaller than 2^(-e), we will loose about e bits 188 in log(1+x) */ 189 if (MPFR_EXP(x) < 0) 190 Nt += -MPFR_EXP(x); 191 192 /* initialize of intermediary variable */ 193 mpfr_init2 (t, Nt); 194 195 /* First computation of log1p */ 196 MPFR_ZIV_INIT (loop, Nt); 197 for (;;) 198 { 199 int k; 200 /* small case: assuming the AGM algorithm used by mpfr_log uses 201 log2(p) steps for a precision of p bits, we try the special 202 variant whenever EXP(x) <= -p/log2(p). */ 203 k = 1 + __gmpfr_int_ceil_log2 (Ny); /* the +1 avoids a division by 0 204 when Ny=1 */ 205 if (MPFR_GET_EXP (x) + 1 <= - (mpfr_exp_t) (Ny / k)) 206 /* this implies EXP(x) <= -1 thus x < 1/2 */ 207 err = Nt - mpfr_log1p_small (t, x); 208 else 209 { 210 /* compute log1p */ 211 inexact = mpfr_add_ui (t, x, 1, MPFR_RNDN); /* 1+x */ 212 /* if inexact = 0, then t = x+1, and the result is simply log(t) */ 213 if (inexact == 0) 214 { 215 inexact = mpfr_log (y, t, rnd_mode); 216 goto end; 217 } 218 mpfr_log (t, t, MPFR_RNDN); /* log(1+x) */ 219 220 /* the error is bounded by (1/2+2^(1-EXP(t))*ulp(t) 221 (cf algorithms.tex) 222 if EXP(t)>=2, then error <= ulp(t) 223 if EXP(t)<=1, then error <= 2^(2-EXP(t))*ulp(t) */ 224 err = Nt - MAX (0, 2 - MPFR_GET_EXP (t)); 225 } 226 227 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 228 break; 229 230 /* increase the precision */ 231 MPFR_ZIV_NEXT (loop, Nt); 232 mpfr_set_prec (t, Nt); 233 } 234 inexact = mpfr_set (y, t, rnd_mode); 235 236 end: 237 MPFR_ZIV_FREE (loop); 238 mpfr_clear (t); 239 } 240 241 MPFR_SAVE_EXPO_FREE (expo); 242 return mpfr_check_range (y, inexact, rnd_mode); 243} 244