1 /* mpfr_eint, mpfr_eint1 -- the exponential integral
2
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel 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 /* 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 and return e such that the absolute error is bound by 2^e ulp(y) */
35 static mpfr_exp_t
mpfr_eint_aux(mpfr_t y,mpfr_srcptr x)36 mpfr_eint_aux (mpfr_t y, mpfr_srcptr x)
37 {
38 mpfr_t eps; /* dynamic (absolute) error bound on t */
39 mpfr_t erru, errs;
40 mpz_t m, s, t, u;
41 mpfr_exp_t e, sizeinbase;
42 mpfr_prec_t w = MPFR_PREC(y);
43 unsigned long k;
44 MPFR_GROUP_DECL (group);
45
46 /* for |x| <= 1, we have S := sum(x^k/k/k!, k=1..infinity) = x + R(x)
47 where |R(x)| <= (x/2)^2/(1-x/2) <= 2*(x/2)^2
48 thus |R(x)/x| <= |x|/2
49 thus if |x| <= 2^(-PREC(y)) we have |S - o(x)| <= ulp(y) */
50
51 if (MPFR_GET_EXP(x) <= - (mpfr_exp_t) w)
52 {
53 mpfr_set (y, x, MPFR_RNDN);
54 return 0;
55 }
56
57 mpz_init (s); /* initializes to 0 */
58 mpz_init (t);
59 mpz_init (u);
60 mpz_init (m);
61 MPFR_GROUP_INIT_3 (group, 31, eps, erru, errs);
62 e = mpfr_get_z_2exp (m, x); /* x = m * 2^e */
63 MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x));
64 if (MPFR_PREC (x) > w)
65 {
66 e += MPFR_PREC (x) - w;
67 mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w);
68 }
69 /* remove trailing zeroes from m: this will speed up much cases where
70 x is a small integer divided by a power of 2 */
71 k = mpz_scan1 (m, 0);
72 mpz_tdiv_q_2exp (m, m, k);
73 e += k;
74 /* initialize t to 2^w */
75 mpz_set_ui (t, 1);
76 mpz_mul_2exp (t, t, w);
77 mpfr_set_ui (eps, 0, MPFR_RNDN); /* eps[0] = 0 */
78 mpfr_set_ui (errs, 0, MPFR_RNDN);
79 for (k = 1;; k++)
80 {
81 /* let eps[k] be the absolute error on t[k]:
82 since t[k] = trunc(t[k-1]*m*2^e/k), we have
83 eps[k+1] <= 1 + eps[k-1]*m*2^e/k + t[k-1]*m*2^(1-w)*2^e/k
84 = 1 + (eps[k-1] + t[k-1]*2^(1-w))*m*2^e/k
85 = 1 + (eps[k-1]*2^(w-1) + t[k-1])*2^(1-w)*m*2^e/k */
86 mpfr_mul_2ui (eps, eps, w - 1, MPFR_RNDU);
87 mpfr_add_z (eps, eps, t, MPFR_RNDU);
88 MPFR_MPZ_SIZEINBASE2 (sizeinbase, m);
89 mpfr_mul_2si (eps, eps, sizeinbase - (w - 1) + e, MPFR_RNDU);
90 mpfr_div_ui (eps, eps, k, MPFR_RNDU);
91 mpfr_add_ui (eps, eps, 1, MPFR_RNDU);
92 mpz_mul (t, t, m);
93 if (e < 0)
94 mpz_tdiv_q_2exp (t, t, -e);
95 else
96 mpz_mul_2exp (t, t, e);
97 mpz_tdiv_q_ui (t, t, k);
98 mpz_tdiv_q_ui (u, t, k);
99 mpz_add (s, s, u);
100 /* the absolute error on u is <= 1 + eps[k]/k */
101 mpfr_div_ui (erru, eps, k, MPFR_RNDU);
102 mpfr_add_ui (erru, erru, 1, MPFR_RNDU);
103 /* and that on s is the sum of all errors on u */
104 mpfr_add (errs, errs, erru, MPFR_RNDU);
105 /* we are done when t is smaller than errs */
106 if (mpz_sgn (t) == 0)
107 sizeinbase = 0;
108 else
109 MPFR_MPZ_SIZEINBASE2 (sizeinbase, t);
110 if (sizeinbase < MPFR_GET_EXP (errs))
111 break;
112 }
113 /* the truncation error is bounded by (|t|+eps)/k*(|x|/k + |x|^2/k^2 + ...)
114 <= (|t|+eps)/k*|x|/(k-|x|) */
115 mpz_abs (t, t);
116 mpfr_add_z (eps, eps, t, MPFR_RNDU);
117 mpfr_div_ui (eps, eps, k, MPFR_RNDU);
118 mpfr_abs (erru, x, MPFR_RNDU); /* |x| */
119 mpfr_mul (eps, eps, erru, MPFR_RNDU);
120 mpfr_ui_sub (erru, k, erru, MPFR_RNDD);
121 if (MPFR_IS_NEG (erru))
122 {
123 /* the truncated series does not converge, return fail */
124 e = w;
125 }
126 else
127 {
128 mpfr_div (eps, eps, erru, MPFR_RNDU);
129 mpfr_add (errs, errs, eps, MPFR_RNDU);
130 mpfr_set_z (y, s, MPFR_RNDN);
131 mpfr_div_2ui (y, y, w, MPFR_RNDN);
132 /* errs was an absolute error bound on s. We must convert it to an error
133 in terms of ulp(y). Since ulp(y) = 2^(EXP(y)-PREC(y)), we must
134 divide the error by 2^(EXP(y)-PREC(y)), but since we divided also
135 y by 2^w = 2^PREC(y), we must simply divide by 2^EXP(y). */
136 e = MPFR_GET_EXP (errs) - MPFR_GET_EXP (y);
137 }
138 MPFR_GROUP_CLEAR (group);
139 mpz_clear (s);
140 mpz_clear (t);
141 mpz_clear (u);
142 mpz_clear (m);
143 return e;
144 }
145
146 /* Return in y an approximation of Ei(x) using the asymptotic expansion:
147 Ei(x) = exp(x)/x * (1 + 1/x + 2/x^2 + ... + k!/x^k + ...)
148 Assumes x >= PREC(y) * log(2).
149 Returns the error bound in terms of ulp(y).
150 */
151 static mpfr_exp_t
mpfr_eint_asympt(mpfr_ptr y,mpfr_srcptr x)152 mpfr_eint_asympt (mpfr_ptr y, mpfr_srcptr x)
153 {
154 mpfr_prec_t p = MPFR_PREC(y);
155 mpfr_t invx, t, err;
156 unsigned long k;
157 mpfr_exp_t err_exp;
158
159 mpfr_init2 (t, p);
160 mpfr_init2 (invx, p);
161 mpfr_init2 (err, 31); /* error in ulps on y */
162 mpfr_ui_div (invx, 1, x, MPFR_RNDN); /* invx = 1/x*(1+u) with |u|<=2^(1-p) */
163 mpfr_set_ui (t, 1, MPFR_RNDN); /* exact */
164 mpfr_set (y, t, MPFR_RNDN);
165 mpfr_set_ui (err, 0, MPFR_RNDN);
166 for (k = 1; MPFR_GET_EXP(t) + (mpfr_exp_t) p > MPFR_GET_EXP(y); k++)
167 {
168 mpfr_mul (t, t, invx, MPFR_RNDN); /* 2 more roundings */
169 mpfr_mul_ui (t, t, k, MPFR_RNDN); /* 1 more rounding: t = k!/x^k*(1+u)^e
170 with u=2^{-p} and |e| <= 3*k */
171 /* we use the fact that |(1+u)^n-1| <= 2*|n*u| for |n*u| <= 1, thus
172 the error on t is less than 6*k*2^{-p}*t <= 6*k*ulp(t) */
173 /* err is in terms of ulp(y): transform it in terms of ulp(t) */
174 mpfr_mul_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), MPFR_RNDU);
175 mpfr_add_ui (err, err, 6 * k, MPFR_RNDU);
176 /* transform back in terms of ulp(y) */
177 mpfr_div_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), MPFR_RNDU);
178 mpfr_add (y, y, t, MPFR_RNDN);
179 }
180 /* add the truncation error bounded by ulp(y): 1 ulp */
181 mpfr_mul (y, y, invx, MPFR_RNDN); /* err <= 2*err + 3/2 */
182 mpfr_exp (t, x, MPFR_RNDN); /* err(t) <= 1/2*ulp(t) */
183 mpfr_mul (y, y, t, MPFR_RNDN); /* again: err <= 2*err + 3/2 */
184 mpfr_mul_2ui (err, err, 2, MPFR_RNDU);
185 mpfr_add_ui (err, err, 8, MPFR_RNDU);
186 err_exp = MPFR_GET_EXP(err);
187 mpfr_clear (t);
188 mpfr_clear (invx);
189 mpfr_clear (err);
190 return err_exp;
191 }
192
193 int
mpfr_eint(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd)194 mpfr_eint (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
195 {
196 int inex;
197 mpfr_t tmp, ump;
198 mpfr_exp_t err, te;
199 mpfr_prec_t prec;
200 MPFR_SAVE_EXPO_DECL (expo);
201 MPFR_ZIV_DECL (loop);
202
203 MPFR_LOG_FUNC (
204 ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
205 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
206
207 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
208 {
209 /* exp(NaN) = exp(-Inf) = NaN */
210 if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_IS_NEG(x)))
211 {
212 MPFR_SET_NAN (y);
213 MPFR_RET_NAN;
214 }
215 /* eint(+inf) = +inf */
216 else if (MPFR_IS_INF (x))
217 {
218 MPFR_SET_INF(y);
219 MPFR_SET_POS(y);
220 MPFR_RET(0);
221 }
222 else /* eint(+/-0) = -Inf */
223 {
224 MPFR_SET_INF(y);
225 MPFR_SET_NEG(y);
226 mpfr_set_divby0 ();
227 MPFR_RET(0);
228 }
229 }
230
231 /* eint(x) = NaN for x < 0 */
232 if (MPFR_IS_NEG(x))
233 {
234 MPFR_SET_NAN (y);
235 MPFR_RET_NAN;
236 }
237
238 MPFR_SAVE_EXPO_MARK (expo);
239
240 /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2).
241 Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax,
242 then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */
243 mpfr_init2 (tmp, 64);
244 mpfr_init2 (ump, 64);
245 mpfr_log (tmp, x, MPFR_RNDU);
246 mpfr_sub (ump, x, tmp, MPFR_RNDD);
247 mpfr_const_log2 (tmp, MPFR_RNDU);
248 mpfr_div (ump, ump, tmp, MPFR_RNDD);
249 /* FIXME: We really need mpfr_set_exp_t and mpfr_cmpfr_exp_t functions. */
250 MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
251 if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0)
252 {
253 mpfr_clear (tmp);
254 mpfr_clear (ump);
255 MPFR_SAVE_EXPO_FREE (expo);
256 return mpfr_overflow (y, rnd, 1);
257 }
258
259 /* Init stuff */
260 prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6;
261
262 /* eint() has a root 0.37250741078136663446..., so if x is near,
263 already take more bits */
264 /* FIXME: do not use native floating-point here. */
265 if (MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */
266 {
267 double d;
268 d = mpfr_get_d (x, MPFR_RNDN) - 0.37250741078136663;
269 d = (d == 0.0) ? -53 : __gmpfr_ceil_log2 (d);
270 prec += -d;
271 }
272
273 mpfr_set_prec (tmp, prec);
274 mpfr_set_prec (ump, prec);
275
276 MPFR_ZIV_INIT (loop, prec); /* Initialize the ZivLoop controler */
277 for (;;) /* Infinite loop */
278 {
279 /* We need that the smallest value of k!/x^k is smaller than 2^(-p).
280 The minimum is obtained for x=k, and it is smaller than e*sqrt(x)/e^x
281 for x>=1. */
282 if (MPFR_GET_EXP (x) > 0 && mpfr_cmp_d (x, ((double) prec +
283 0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0)
284 err = mpfr_eint_asympt (tmp, x);
285 else
286 {
287 err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */
288 te = MPFR_GET_EXP(tmp);
289 mpfr_const_euler (ump, MPFR_RNDN); /* 0.577 -> EXP(ump)=0 */
290 mpfr_add (tmp, tmp, ump, MPFR_RNDN);
291 /* error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err)
292 <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp))
293 <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))) */
294 err = MAX(1, te + err + 2) - MPFR_GET_EXP(tmp);
295 err = MAX(0, err);
296 te = MPFR_GET_EXP(tmp);
297 mpfr_log (ump, x, MPFR_RNDN);
298 mpfr_add (tmp, tmp, ump, MPFR_RNDN);
299 /* same formula as above, except now EXP(ump) is not 0 */
300 err += te + 1;
301 if (MPFR_LIKELY (!MPFR_IS_ZERO (ump)))
302 err = MAX (MPFR_GET_EXP (ump), err);
303 err = MAX(0, err - MPFR_GET_EXP (tmp));
304 }
305 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
306 break;
307 MPFR_ZIV_NEXT (loop, prec); /* Increase used precision */
308 mpfr_set_prec (tmp, prec);
309 mpfr_set_prec (ump, prec);
310 }
311 MPFR_ZIV_FREE (loop); /* Free the ZivLoop Controler */
312
313 inex = mpfr_set (y, tmp, rnd); /* Set y to the computed value */
314 mpfr_clear (tmp);
315 mpfr_clear (ump);
316
317 MPFR_SAVE_EXPO_FREE (expo);
318 return mpfr_check_range (y, inex, rnd);
319 }
320