1 /* mpfr_eint, mpfr_eint1 -- the exponential integral
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
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
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.
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) */
36 mpfr_eint_aux (mpfr_t y
, mpfr_srcptr x
)
38 mpfr_t eps
; /* dynamic (absolute) error bound on t */
41 mpfr_exp_t e
, sizeinbase
;
42 mpfr_prec_t w
= MPFR_PREC(y
);
44 MPFR_GROUP_DECL (group
);
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) */
51 if (MPFR_GET_EXP(x
) <= - (mpfr_exp_t
) w
)
53 mpfr_set (y
, x
, MPFR_RNDN
);
57 mpz_init (s
); /* initializes to 0 */
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
)
66 e
+= MPFR_PREC (x
) - w
;
67 mpz_tdiv_q_2exp (m
, m
, MPFR_PREC (x
) - w
);
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 */
72 mpz_tdiv_q_2exp (m
, m
, k
);
74 /* initialize t to 2^w */
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
);
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
);
94 mpz_tdiv_q_2exp (t
, t
, -e
);
96 mpz_mul_2exp (t
, t
, e
);
97 mpz_tdiv_q_ui (t
, t
, k
);
98 mpz_tdiv_q_ui (u
, t
, k
);
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)
109 MPFR_MPZ_SIZEINBASE2 (sizeinbase
, t
);
110 if (sizeinbase
< MPFR_GET_EXP (errs
))
113 /* the truncation error is bounded by (|t|+eps)/k*(|x|/k + |x|^2/k^2 + ...)
114 <= (|t|+eps)/k*|x|/(k-|x|) */
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
))
123 /* the truncated series does not converge, return fail */
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
);
138 MPFR_GROUP_CLEAR (group
);
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).
152 mpfr_eint_asympt (mpfr_ptr y
, mpfr_srcptr x
)
154 mpfr_prec_t p
= MPFR_PREC(y
);
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
++)
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
);
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
);
194 mpfr_eint (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd
)
200 MPFR_SAVE_EXPO_DECL (expo
);
201 MPFR_ZIV_DECL (loop
);
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
));
207 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
209 /* exp(NaN) = exp(-Inf) = NaN */
210 if (MPFR_IS_NAN (x
) || (MPFR_IS_INF (x
) && MPFR_IS_NEG(x
)))
215 /* eint(+inf) = +inf */
216 else if (MPFR_IS_INF (x
))
222 else /* eint(+/-0) = -Inf */
231 /* eint(x) = NaN for x < 0 */
238 MPFR_SAVE_EXPO_MARK (expo
);
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)
255 MPFR_SAVE_EXPO_FREE (expo
);
256 return mpfr_overflow (y
, rnd
, 1);
260 prec
= MPFR_PREC (y
) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y
)) + 6;
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 */
268 d
= mpfr_get_d (x
, MPFR_RNDN
) - 0.37250741078136663;
269 d
= (d
== 0.0) ? -53 : __gmpfr_ceil_log2 (d
);
273 mpfr_set_prec (tmp
, prec
);
274 mpfr_set_prec (ump
, prec
);
276 MPFR_ZIV_INIT (loop
, prec
); /* Initialize the ZivLoop controler */
277 for (;;) /* Infinite loop */
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
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
);
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
);
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 */
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
));
305 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp
, prec
- err
, MPFR_PREC (y
), rnd
)))
307 MPFR_ZIV_NEXT (loop
, prec
); /* Increase used precision */
308 mpfr_set_prec (tmp
, prec
);
309 mpfr_set_prec (ump
, prec
);
311 MPFR_ZIV_FREE (loop
); /* Free the ZivLoop Controler */
313 inex
= mpfr_set (y
, tmp
, rnd
); /* Set y to the computed value */
317 MPFR_SAVE_EXPO_FREE (expo
);
318 return mpfr_check_range (y
, inex
, rnd
);