1 /* mpfr_log1p -- Compute log(1+x)
2
3 Copyright 2001-2020 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 /* Put in y an approximation of log(1+x) for x small.
27 We assume |x| < 1, in which case:
28 |x/2| <= |log(1+x)| = |x - x^2/2 + x^3/3 - x^4/4 + ...| <= |x|.
29 Return k such that the error is bounded by 2^k*ulp(y).
30 */
31 static int
mpfr_log1p_small(mpfr_ptr y,mpfr_srcptr x)32 mpfr_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) <= 0); /* ensures |x| < 1 */
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
101 int
mpfr_log1p(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)102 mpfr_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) <= - (mpfr_exp_t) (Ny / k))
206 /* this implies EXP(x) <= 0 thus x < 1 */
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