1 /* mpfr_eint, mpfr_eint1 -- the exponential integral
2
3 Copyright 2005-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 /* eint1(x) = -gamma - log(x) - sum((-1)^k*z^k/k/k!, k=1..infinity) for x > 0
27 = - eint(-x) for x < 0
28 where
29 eint (x) = gamma + log(x) + sum(z^k/k/k!, k=1..infinity) for x > 0
30 eint (x) is undefined for x < 0.
31 */
32
33 /* Compute in y an approximation of sum(x^k/k/k!, k=1..infinity),
34 assuming x != 0, and return e such that the absolute error is
35 bounded by 2^e ulp(y).
36 Return PREC(y) when the truncated series does not converge.
37 */
38 static mpfr_exp_t
mpfr_eint_aux(mpfr_ptr y,mpfr_srcptr x)39 mpfr_eint_aux (mpfr_ptr y, mpfr_srcptr x)
40 {
41 mpfr_t eps; /* dynamic (absolute) error bound on t */
42 mpfr_t erru, errs;
43 mpz_t m, s, t, u;
44 mpfr_exp_t e, sizeinbase;
45 mpfr_prec_t w = MPFR_PREC(y);
46 unsigned long k;
47 MPFR_GROUP_DECL (group);
48
49 MPFR_LOG_FUNC (
50 ("x[%Pd]=%.*Rg", mpfr_get_prec (x), mpfr_log_prec, x),
51 ("y[%Pd]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y));
52
53 /* for |x| <= 1, we have S := sum(x^k/k/k!, k=1..infinity) = x + R(x)
54 where |R(x)| <= (x/2)^2/(1-|x|/2) <= 2*(x/2)^2
55 thus |R(x)/x| <= |x|/2
56 thus if |x| <= 2^(-PREC(y)) we have |S - o(x)| <= ulp(y) */
57
58 if (MPFR_GET_EXP(x) <= - (mpfr_exp_t) w)
59 {
60 mpfr_set (y, x, MPFR_RNDN);
61 return 0;
62 }
63
64 mpz_init (s); /* initializes to 0 */
65 mpz_init (t);
66 mpz_init (u);
67 mpz_init (m);
68 MPFR_GROUP_INIT_3 (group, 31, eps, erru, errs);
69 e = mpfr_get_z_2exp (m, x); /* x = m * 2^e with m != 0 */
70 MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e));
71 MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x)); /* since m != 0 */
72 if (MPFR_PREC (x) > w)
73 {
74 e += MPFR_PREC (x) - w;
75 mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w); /* one still has m != 0 */
76 MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e));
77 }
78 /* Remove trailing zeroes from m: this will speed up much cases where
79 x is a small integer divided by a power of 2.
80 Note: As shown above, m != 0. This is needed for the "e += ..." below,
81 otherwise n would take the largest value of mp_bitcnt_t and could be
82 too large. */
83 {
84 mp_bitcnt_t n = mpz_scan1 (m, 0);
85 mpz_tdiv_q_2exp (m, m, n);
86 /* Since one initially has mpz_sizeinbase (m, 2) == MPFR_PREC (x)
87 and m has not increased, one can deduce that n <= MPFR_PREC (x),
88 so that the cast to mpfr_prec_t is valid. This cast is needed to
89 ensure that the operand e of the addition below is not converted
90 to an unsigned integer type, which could yield incorrect results
91 with some C implementations. */
92 MPFR_ASSERTD (n <= MPFR_PREC (x));
93 e += (mpfr_prec_t) n;
94 }
95 /* initialize t to 2^w */
96 mpz_set_ui (t, 1);
97 mpz_mul_2exp (t, t, w);
98 mpfr_set_ui (eps, 0, MPFR_RNDN); /* eps[0] = 0 */
99 mpfr_set_ui (errs, 0, MPFR_RNDN); /* maximal error on s */
100 for (k = 1;; k++)
101 {
102 /* let t[k] = x^k/k/k!, and eps[k] be the absolute error on t[k]:
103 since t[k] = trunc(t[k-1]*m*2^e/k), we have
104 eps[k+1] <= 1 + eps[k-1]*|m|*2^e/k + |t[k-1]|*|m|*2^(1-w)*2^e/k
105 = 1 + (eps[k-1] + |t[k-1]|*2^(1-w))*|m|*2^e/k
106 = 1 + (eps[k-1]*2^(w-1) + |t[k-1]|)*2^(1-w)*|m|*2^e/k */
107 mpfr_mul_2ui (eps, eps, w - 1, MPFR_RNDU);
108 if (mpz_sgn (t) >= 0)
109 mpfr_add_z (eps, eps, t, MPFR_RNDU);
110 else
111 mpfr_sub_z (eps, eps, t, MPFR_RNDU);
112 MPFR_MPZ_SIZEINBASE2 (sizeinbase, m);
113 mpfr_mul_2si (eps, eps, sizeinbase - (w - 1) + e, MPFR_RNDU);
114 mpfr_div_ui (eps, eps, k, MPFR_RNDU);
115 mpfr_add_ui (eps, eps, 1, MPFR_RNDU);
116 mpz_mul (t, t, m);
117 if (e < 0)
118 mpz_tdiv_q_2exp (t, t, -e);
119 else
120 mpz_mul_2exp (t, t, e);
121 mpz_tdiv_q_ui (t, t, k);
122 mpz_tdiv_q_ui (u, t, k);
123 mpz_add (s, s, u);
124 /* the absolute error on u is <= 1 + eps[k]/k */
125 mpfr_div_ui (erru, eps, k, MPFR_RNDU);
126 mpfr_add_ui (erru, erru, 1, MPFR_RNDU);
127 /* and that on s is the sum of all errors on u */
128 mpfr_add (errs, errs, erru, MPFR_RNDU);
129 /* we are done when t is smaller than errs */
130 if (mpz_sgn (t) == 0)
131 sizeinbase = 0;
132 else
133 MPFR_MPZ_SIZEINBASE2 (sizeinbase, t);
134 if (sizeinbase < MPFR_GET_EXP (errs))
135 break;
136 }
137 /* the truncation error is bounded by (|t|+eps)/k*(|x|/k + |x|^2/k^2 + ...)
138 <= (|t|+eps)/k*|x|/(k-|x|) */
139 mpz_abs (t, t);
140 mpfr_add_z (eps, eps, t, MPFR_RNDU);
141 mpfr_div_ui (eps, eps, k, MPFR_RNDU);
142 mpfr_abs (erru, x, MPFR_RNDU); /* |x| */
143 mpfr_mul (eps, eps, erru, MPFR_RNDU);
144 mpfr_ui_sub (erru, k, erru, MPFR_RNDD);
145 if (MPFR_IS_NEG (erru))
146 {
147 /* the truncated series does not converge, return fail */
148 e = w;
149 }
150 else
151 {
152 mpfr_div (eps, eps, erru, MPFR_RNDU);
153 mpfr_add (errs, errs, eps, MPFR_RNDU);
154 mpfr_set_z (y, s, MPFR_RNDN);
155 mpfr_div_2ui (y, y, w, MPFR_RNDN);
156 /* errs was an absolute error bound on s. We must convert it to an error
157 in terms of ulp(y). Since ulp(y) = 2^(EXP(y)-PREC(y)), we must
158 divide the error by 2^(EXP(y)-PREC(y)), but since we divided also
159 y by 2^w = 2^PREC(y), we must simply divide by 2^EXP(y). */
160 e = MPFR_GET_EXP (errs) - MPFR_GET_EXP (y);
161 }
162 MPFR_GROUP_CLEAR (group);
163 mpz_clear (s);
164 mpz_clear (t);
165 mpz_clear (u);
166 mpz_clear (m);
167 MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e));
168 return e;
169 }
170
171 /* Return in y an approximation of Ei(x) using the asymptotic expansion:
172 Ei(x) = exp(x)/x * (1 + 1/x + 2/x^2 + ... + k!/x^k + ...)
173 Assumes |x| >= PREC(y) * log(2).
174 Returns the error bound in terms of ulp(y).
175 */
176 static mpfr_exp_t
mpfr_eint_asympt(mpfr_ptr y,mpfr_srcptr x)177 mpfr_eint_asympt (mpfr_ptr y, mpfr_srcptr x)
178 {
179 mpfr_prec_t p = MPFR_PREC(y);
180 mpfr_t invx, t, err;
181 unsigned long k;
182 mpfr_exp_t err_exp;
183
184 MPFR_LOG_FUNC (
185 ("x[%Pd]=%.*Rg", mpfr_get_prec (x), mpfr_log_prec, x),
186 ("err_exp=%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) err_exp));
187
188 mpfr_init2 (t, p);
189 mpfr_init2 (invx, p);
190 mpfr_init2 (err, 31); /* error in ulps on y */
191 mpfr_ui_div (invx, 1, x, MPFR_RNDN); /* invx = 1/x*(1+u) with |u|<=2^(1-p) */
192 mpfr_set_ui (t, 1, MPFR_RNDN); /* exact */
193 mpfr_set (y, t, MPFR_RNDN);
194 mpfr_set_ui (err, 0, MPFR_RNDN);
195 for (k = 1; MPFR_GET_EXP(t) + (mpfr_exp_t) p > MPFR_GET_EXP(y); k++)
196 {
197 mpfr_mul (t, t, invx, MPFR_RNDN); /* 2 more roundings */
198 mpfr_mul_ui (t, t, k, MPFR_RNDN); /* 1 more rounding: t = k!/x^k*(1+u)^e
199 with u=2^{-p} and |e| <= 3*k */
200 /* we use the fact that |(1+u)^n-1| <= 2*|n*u| for |n*u| <= 1, thus
201 the error on t is less than 6*k*2^{-p}*t <= 6*k*ulp(t) */
202 /* err is in terms of ulp(y): transform it in terms of ulp(t) */
203 mpfr_mul_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), MPFR_RNDU);
204 mpfr_add_ui (err, err, 6 * k, MPFR_RNDU);
205 /* transform back in terms of ulp(y) */
206 mpfr_div_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), MPFR_RNDU);
207 mpfr_add (y, y, t, MPFR_RNDN);
208 }
209 /* add the truncation error bounded by ulp(y): 1 ulp */
210 mpfr_mul (y, y, invx, MPFR_RNDN); /* err <= 2*err + 3/2 */
211 mpfr_exp (t, x, MPFR_RNDN); /* err(t) <= 1/2*ulp(t) */
212 mpfr_mul (y, y, t, MPFR_RNDN); /* again: err <= 2*err + 3/2 */
213 mpfr_mul_2ui (err, err, 2, MPFR_RNDU);
214 mpfr_add_ui (err, err, 8, MPFR_RNDU);
215 err_exp = MPFR_GET_EXP(err);
216 mpfr_clear (t);
217 mpfr_clear (invx);
218 mpfr_clear (err);
219 return err_exp;
220 }
221
222 /* mpfr_eint returns Ei(x) for x >= 0,
223 and -E1(-x) for x < 0, following https://dlmf.nist.gov/6.2 */
224 int
mpfr_eint(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd)225 mpfr_eint (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
226 {
227 int inex;
228 mpfr_t tmp, ump, x_abs;
229 mpfr_exp_t err, te;
230 mpfr_prec_t prec;
231 MPFR_SAVE_EXPO_DECL (expo);
232 MPFR_ZIV_DECL (loop);
233
234 MPFR_LOG_FUNC (
235 ("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
236 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
237
238 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
239 {
240 if (MPFR_IS_NAN (x))
241 {
242 MPFR_SET_NAN (y);
243 MPFR_RET_NAN;
244 }
245 else if (MPFR_IS_INF (x))
246 {
247 /* eint(+inf) = +inf and eint(-inf) = -0 */
248 if (MPFR_IS_POS (x))
249 {
250 MPFR_SET_INF(y);
251 MPFR_SET_POS(y);
252 }
253 else
254 {
255 MPFR_SET_ZERO(y);
256 MPFR_SET_NEG(y);
257 }
258 MPFR_RET(0);
259 }
260 else /* eint(+/-0) = -Inf */
261 {
262 MPFR_SET_INF(y);
263 MPFR_SET_NEG(y);
264 MPFR_SET_DIVBY0 ();
265 MPFR_RET(0);
266 }
267 }
268
269 MPFR_TMP_INIT_ABS (x_abs, x);
270
271 MPFR_SAVE_EXPO_MARK (expo);
272
273 /* Init stuff */
274 prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6;
275 mpfr_init2 (tmp, 64);
276 mpfr_init2 (ump, 64);
277
278 /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2).
279 Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax,
280 then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */
281 if (MPFR_IS_POS(x))
282 {
283 mpfr_log (tmp, x, MPFR_RNDU);
284 mpfr_sub (ump, x, tmp, MPFR_RNDD);
285 mpfr_div (ump, ump, __gmpfr_const_log2_RNDU, MPFR_RNDD);
286 /* FIXME: We really need a mpfr_cmp_exp_t function. */
287 MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
288 if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0)
289 {
290 mpfr_clear (tmp);
291 mpfr_clear (ump);
292 MPFR_SAVE_EXPO_FREE (expo);
293 return mpfr_overflow (y, rnd, 1);
294 }
295 }
296
297 /* Since E1(x) <= exp(-x) for x >= 1, we have log2(E1(x)) <= -x/log(2).
298 Let's compute k >= -x/log(2) in a low precision. If k < emin
299 then log2(E1(x)) <= emin-1, and E1(x) <= 2^(emin-1): it underflows. */
300 if (MPFR_IS_NEG(x) && MPFR_GET_EXP(x) >= 1)
301 {
302 mpfr_div (ump, x, __gmpfr_const_log2_RNDD, MPFR_RNDU);
303 MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN);
304 if (mpfr_cmp_si (ump, __gmpfr_emin) < 0)
305 {
306 mpfr_clear (tmp);
307 mpfr_clear (ump);
308 MPFR_SAVE_EXPO_FREE (expo);
309 return mpfr_underflow (y, rnd, -1);
310 }
311 }
312
313 /* eint() has a root 0.37250741078136663446...,
314 so if x is near, already take more bits */
315 if (MPFR_IS_POS(x) && MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */
316 {
317 mpfr_t y;
318 mpfr_init2 (y, 32);
319 /* 1599907147/2^32 is a 32-bit approximation of 0.37250741078136663446 */
320 mpfr_set_ui_2exp (y, 1599907147UL, -32, MPFR_RNDN);
321 mpfr_sub (y, x, y, MPFR_RNDN);
322 prec += (mpfr_zero_p (y)) ? 32
323 : mpfr_get_exp (y) < 0 ? -mpfr_get_exp (y) : 0;
324 mpfr_clear (y);
325 }
326
327 mpfr_set_prec (tmp, prec);
328 mpfr_set_prec (ump, prec);
329
330 MPFR_ZIV_INIT (loop, prec); /* Initialize the ZivLoop controller */
331 for (;;) /* Infinite loop */
332 {
333 /* For the asymptotic expansion to work, we need that the smallest
334 value of k!/|x|^k is smaller than 2^(-p). The minimum is obtained for
335 x=k, and it is smaller than e*sqrt(x)/e^x for x>=1. */
336 if (MPFR_GET_EXP (x) > 0 &&
337 mpfr_cmp_d (x_abs, ((double) prec +
338 0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0)
339 err = mpfr_eint_asympt (tmp, x);
340 else
341 {
342 err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */
343 te = MPFR_GET_EXP(tmp);
344 mpfr_const_euler (ump, MPFR_RNDN); /* 0.577 -> EXP(ump)=0 */
345 mpfr_add (tmp, tmp, ump, MPFR_RNDN);
346 /* If tmp <> 0:
347 error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err)
348 <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp))
349 <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))).
350 If tmp = 0 we can use the same bound, replacing
351 EXP(tmp) by EXP(ump). */
352 err = MAX(1, te + err + 2);
353 te = MPFR_IS_ZERO(tmp) ? MPFR_GET_EXP(ump) : MPFR_GET_EXP(tmp);
354 err = err - te;
355 err = MAX(0, err);
356 mpfr_log (ump, x_abs, MPFR_RNDN);
357 mpfr_add (tmp, tmp, ump, MPFR_RNDN);
358 /* same formula as above, except now EXP(ump) is not 0 */
359 err += te + 1;
360 if (MPFR_LIKELY (!MPFR_IS_ZERO (ump)))
361 err = MAX (MPFR_GET_EXP (ump), err);
362 /* if tmp is zero, we surely cannot round correctly */
363 err = (MPFR_IS_ZERO(tmp)) ? prec : MAX(0, err - MPFR_GET_EXP (tmp));
364 }
365 /* Note: we assume here that MPFR_CAN_ROUND returns the same result
366 for rnd and MPFR_INVERT_RND(rnd) */
367 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
368 break;
369 MPFR_ZIV_NEXT (loop, prec); /* Increase used precision */
370 mpfr_set_prec (tmp, prec);
371 mpfr_set_prec (ump, prec);
372 }
373 MPFR_ZIV_FREE (loop); /* Free the ZivLoop Controller */
374
375 /* Set y to the computed value */
376 inex = mpfr_set (y, tmp, rnd);
377 mpfr_clear (tmp);
378 mpfr_clear (ump);
379
380 MPFR_SAVE_EXPO_FREE (expo);
381 return mpfr_check_range (y, inex, rnd);
382 }
383