1/* mpfr_eint, mpfr_eint1 -- the exponential integral 2 3Copyright 2005-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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 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*/ 38static mpfr_exp_t 39mpfr_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*/ 176static mpfr_exp_t 177mpfr_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 */ 224int 225mpfr_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