1/* mpfr_eint, mpfr_eint1 -- the exponential integral 2 3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 20http://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/* 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) */ 35static mpfr_exp_t 36mpfr_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*/ 151static mpfr_exp_t 152mpfr_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 193int 194mpfr_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 (("x[%#R]=%R rnd=%d", x, x, rnd), 204 ("y[%#R]=%R inexact=%d", y, y, inex)); 205 206 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 207 { 208 /* exp(NaN) = exp(-Inf) = NaN */ 209 if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_IS_NEG(x))) 210 { 211 MPFR_SET_NAN (y); 212 MPFR_RET_NAN; 213 } 214 /* eint(+inf) = +inf */ 215 else if (MPFR_IS_INF (x)) 216 { 217 MPFR_SET_INF(y); 218 MPFR_SET_POS(y); 219 MPFR_RET(0); 220 } 221 else /* eint(+/-0) = -Inf */ 222 { 223 MPFR_SET_INF(y); 224 MPFR_SET_NEG(y); 225 MPFR_RET(0); 226 } 227 } 228 229 /* eint(x) = NaN for x < 0 */ 230 if (MPFR_IS_NEG(x)) 231 { 232 MPFR_SET_NAN (y); 233 MPFR_RET_NAN; 234 } 235 236 MPFR_SAVE_EXPO_MARK (expo); 237 238 /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2). 239 Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax, 240 then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */ 241 mpfr_init2 (tmp, 64); 242 mpfr_init2 (ump, 64); 243 mpfr_log (tmp, x, MPFR_RNDU); 244 mpfr_sub (ump, x, tmp, MPFR_RNDD); 245 mpfr_const_log2 (tmp, MPFR_RNDU); 246 mpfr_div (ump, ump, tmp, MPFR_RNDD); 247 /* FIXME: We really need mpfr_set_exp_t and mpfr_cmpfr_exp_t functions. */ 248 MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX); 249 if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0) 250 { 251 mpfr_clear (tmp); 252 mpfr_clear (ump); 253 MPFR_SAVE_EXPO_FREE (expo); 254 return mpfr_overflow (y, rnd, 1); 255 } 256 257 /* Init stuff */ 258 prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6; 259 260 /* eint() has a root 0.37250741078136663446..., so if x is near, 261 already take more bits */ 262 if (MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */ 263 { 264 double d; 265 d = mpfr_get_d (x, MPFR_RNDN) - 0.37250741078136663; 266 d = (d == 0.0) ? -53 : __gmpfr_ceil_log2 (d); 267 prec += -d; 268 } 269 270 mpfr_set_prec (tmp, prec); 271 mpfr_set_prec (ump, prec); 272 273 MPFR_ZIV_INIT (loop, prec); /* Initialize the ZivLoop controler */ 274 for (;;) /* Infinite loop */ 275 { 276 /* We need that the smallest value of k!/x^k is smaller than 2^(-p). 277 The minimum is obtained for x=k, and it is smaller than e*sqrt(x)/e^x 278 for x>=1. */ 279 if (MPFR_GET_EXP (x) > 0 && mpfr_cmp_d (x, ((double) prec + 280 0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0) 281 err = mpfr_eint_asympt (tmp, x); 282 else 283 { 284 err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */ 285 te = MPFR_GET_EXP(tmp); 286 mpfr_const_euler (ump, MPFR_RNDN); /* 0.577 -> EXP(ump)=0 */ 287 mpfr_add (tmp, tmp, ump, MPFR_RNDN); 288 /* error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err) 289 <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp)) 290 <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))) */ 291 err = MAX(1, te + err + 2) - MPFR_GET_EXP(tmp); 292 err = MAX(0, err); 293 te = MPFR_GET_EXP(tmp); 294 mpfr_log (ump, x, MPFR_RNDN); 295 mpfr_add (tmp, tmp, ump, MPFR_RNDN); 296 /* same formula as above, except now EXP(ump) is not 0 */ 297 err += te + 1; 298 if (MPFR_LIKELY (!MPFR_IS_ZERO (ump))) 299 err = MAX (MPFR_GET_EXP (ump), err); 300 err = MAX(0, err - MPFR_GET_EXP (tmp)); 301 } 302 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd))) 303 break; 304 MPFR_ZIV_NEXT (loop, prec); /* Increase used precision */ 305 mpfr_set_prec (tmp, prec); 306 mpfr_set_prec (ump, prec); 307 } 308 MPFR_ZIV_FREE (loop); /* Free the ZivLoop Controler */ 309 310 inex = mpfr_set (y, tmp, rnd); /* Set y to the computed value */ 311 mpfr_clear (tmp); 312 mpfr_clear (ump); 313 314 MPFR_SAVE_EXPO_FREE (expo); 315 return mpfr_check_range (y, inex, rnd); 316} 317