1 /* mpfr_eint, mpfr_eint1 -- the exponential integral 2 3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire 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 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 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 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