1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006-2018 Free Software Foundation, Inc. 4 // 5 // This file is part of the GNU ISO C++ Library. This library is free 6 // software; you can redistribute it and/or modify it under the 7 // terms of the GNU General Public License as published by the 8 // Free Software Foundation; either version 3, or (at your option) 9 // any later version. 10 // 11 // This library is distributed in the hope that it will be useful, 12 // but WITHOUT ANY WARRANTY; without even the implied warranty of 13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 // GNU General Public License for more details. 15 // 16 // Under Section 7 of GPL version 3, you are granted additional 17 // permissions described in the GCC Runtime Library Exception, version 18 // 3.1, as published by the Free Software Foundation. 19 20 // You should have received a copy of the GNU General Public License and 21 // a copy of the GCC Runtime Library Exception along with this program; 22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 23 // <http://www.gnu.org/licenses/>. 24 25 /** @file tr1/exp_integral.tcc 26 * This is an internal header file, included by other library headers. 27 * Do not attempt to use it directly. @headername{tr1/cmath} 28 */ 29 30 // 31 // ISO C++ 14882 TR1: 5.2 Special functions 32 // 33 34 // Written by Edward Smith-Rowland based on: 35 // 36 // (1) Handbook of Mathematical Functions, 37 // Ed. by Milton Abramowitz and Irene A. Stegun, 38 // Dover Publications, New-York, Section 5, pp. 228-251. 39 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 40 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, 41 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), 42 // 2nd ed, pp. 222-225. 43 // 44 45 #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC 46 #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1 47 48 #include "special_function_util.h" 49 50 namespace std _GLIBCXX_VISIBILITY(default) 51 { 52 _GLIBCXX_BEGIN_NAMESPACE_VERSION 53 54 #if _GLIBCXX_USE_STD_SPEC_FUNCS 55 #elif defined(_GLIBCXX_TR1_CMATH) 56 namespace tr1 57 { 58 #else 59 # error do not include this header directly, use <cmath> or <tr1/cmath> 60 #endif 61 // [5.2] Special functions 62 63 // Implementation-space details. 64 namespace __detail 65 { 66 template<typename _Tp> _Tp __expint_E1(_Tp); 67 68 /** 69 * @brief Return the exponential integral @f$ E_1(x) @f$ 70 * by series summation. This should be good 71 * for @f$ x < 1 @f$. 72 * 73 * The exponential integral is given by 74 * \f[ 75 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt 76 * \f] 77 * 78 * @param __x The argument of the exponential integral function. 79 * @return The exponential integral. 80 */ 81 template<typename _Tp> 82 _Tp 83 __expint_E1_series(_Tp __x) 84 { 85 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 86 _Tp __term = _Tp(1); 87 _Tp __esum = _Tp(0); 88 _Tp __osum = _Tp(0); 89 const unsigned int __max_iter = 1000; 90 for (unsigned int __i = 1; __i < __max_iter; ++__i) 91 { 92 __term *= - __x / __i; 93 if (std::abs(__term) < __eps) 94 break; 95 if (__term >= _Tp(0)) 96 __esum += __term / __i; 97 else 98 __osum += __term / __i; 99 } 100 101 return - __esum - __osum 102 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x); 103 } 104 105 106 /** 107 * @brief Return the exponential integral @f$ E_1(x) @f$ 108 * by asymptotic expansion. 109 * 110 * The exponential integral is given by 111 * \f[ 112 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 113 * \f] 114 * 115 * @param __x The argument of the exponential integral function. 116 * @return The exponential integral. 117 */ 118 template<typename _Tp> 119 _Tp 120 __expint_E1_asymp(_Tp __x) 121 { 122 _Tp __term = _Tp(1); 123 _Tp __esum = _Tp(1); 124 _Tp __osum = _Tp(0); 125 const unsigned int __max_iter = 1000; 126 for (unsigned int __i = 1; __i < __max_iter; ++__i) 127 { 128 _Tp __prev = __term; 129 __term *= - __i / __x; 130 if (std::abs(__term) > std::abs(__prev)) 131 break; 132 if (__term >= _Tp(0)) 133 __esum += __term; 134 else 135 __osum += __term; 136 } 137 138 return std::exp(- __x) * (__esum + __osum) / __x; 139 } 140 141 142 /** 143 * @brief Return the exponential integral @f$ E_n(x) @f$ 144 * by series summation. 145 * 146 * The exponential integral is given by 147 * \f[ 148 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 149 * \f] 150 * 151 * @param __n The order of the exponential integral function. 152 * @param __x The argument of the exponential integral function. 153 * @return The exponential integral. 154 */ 155 template<typename _Tp> 156 _Tp 157 __expint_En_series(unsigned int __n, _Tp __x) 158 { 159 const unsigned int __max_iter = 1000; 160 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 161 const int __nm1 = __n - 1; 162 _Tp __ans = (__nm1 != 0 163 ? _Tp(1) / __nm1 : -std::log(__x) 164 - __numeric_constants<_Tp>::__gamma_e()); 165 _Tp __fact = _Tp(1); 166 for (int __i = 1; __i <= __max_iter; ++__i) 167 { 168 __fact *= -__x / _Tp(__i); 169 _Tp __del; 170 if ( __i != __nm1 ) 171 __del = -__fact / _Tp(__i - __nm1); 172 else 173 { 174 _Tp __psi = -__numeric_constants<_Tp>::gamma_e(); 175 for (int __ii = 1; __ii <= __nm1; ++__ii) 176 __psi += _Tp(1) / _Tp(__ii); 177 __del = __fact * (__psi - std::log(__x)); 178 } 179 __ans += __del; 180 if (std::abs(__del) < __eps * std::abs(__ans)) 181 return __ans; 182 } 183 std::__throw_runtime_error(__N("Series summation failed " 184 "in __expint_En_series.")); 185 } 186 187 188 /** 189 * @brief Return the exponential integral @f$ E_n(x) @f$ 190 * by continued fractions. 191 * 192 * The exponential integral is given by 193 * \f[ 194 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 195 * \f] 196 * 197 * @param __n The order of the exponential integral function. 198 * @param __x The argument of the exponential integral function. 199 * @return The exponential integral. 200 */ 201 template<typename _Tp> 202 _Tp 203 __expint_En_cont_frac(unsigned int __n, _Tp __x) 204 { 205 const unsigned int __max_iter = 1000; 206 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 207 const _Tp __fp_min = std::numeric_limits<_Tp>::min(); 208 const int __nm1 = __n - 1; 209 _Tp __b = __x + _Tp(__n); 210 _Tp __c = _Tp(1) / __fp_min; 211 _Tp __d = _Tp(1) / __b; 212 _Tp __h = __d; 213 for ( unsigned int __i = 1; __i <= __max_iter; ++__i ) 214 { 215 _Tp __a = -_Tp(__i * (__nm1 + __i)); 216 __b += _Tp(2); 217 __d = _Tp(1) / (__a * __d + __b); 218 __c = __b + __a / __c; 219 const _Tp __del = __c * __d; 220 __h *= __del; 221 if (std::abs(__del - _Tp(1)) < __eps) 222 { 223 const _Tp __ans = __h * std::exp(-__x); 224 return __ans; 225 } 226 } 227 std::__throw_runtime_error(__N("Continued fraction failed " 228 "in __expint_En_cont_frac.")); 229 } 230 231 232 /** 233 * @brief Return the exponential integral @f$ E_n(x) @f$ 234 * by recursion. Use upward recursion for @f$ x < n @f$ 235 * and downward recursion (Miller's algorithm) otherwise. 236 * 237 * The exponential integral is given by 238 * \f[ 239 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 240 * \f] 241 * 242 * @param __n The order of the exponential integral function. 243 * @param __x The argument of the exponential integral function. 244 * @return The exponential integral. 245 */ 246 template<typename _Tp> 247 _Tp 248 __expint_En_recursion(unsigned int __n, _Tp __x) 249 { 250 _Tp __En; 251 _Tp __E1 = __expint_E1(__x); 252 if (__x < _Tp(__n)) 253 { 254 // Forward recursion is stable only for n < x. 255 __En = __E1; 256 for (unsigned int __j = 2; __j < __n; ++__j) 257 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1); 258 } 259 else 260 { 261 // Backward recursion is stable only for n >= x. 262 __En = _Tp(1); 263 const int __N = __n + 20; // TODO: Check this starting number. 264 _Tp __save = _Tp(0); 265 for (int __j = __N; __j > 0; --__j) 266 { 267 __En = (std::exp(-__x) - __j * __En) / __x; 268 if (__j == __n) 269 __save = __En; 270 } 271 _Tp __norm = __En / __E1; 272 __En /= __norm; 273 } 274 275 return __En; 276 } 277 278 /** 279 * @brief Return the exponential integral @f$ Ei(x) @f$ 280 * by series summation. 281 * 282 * The exponential integral is given by 283 * \f[ 284 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 285 * \f] 286 * 287 * @param __x The argument of the exponential integral function. 288 * @return The exponential integral. 289 */ 290 template<typename _Tp> 291 _Tp 292 __expint_Ei_series(_Tp __x) 293 { 294 _Tp __term = _Tp(1); 295 _Tp __sum = _Tp(0); 296 const unsigned int __max_iter = 1000; 297 for (unsigned int __i = 1; __i < __max_iter; ++__i) 298 { 299 __term *= __x / __i; 300 __sum += __term / __i; 301 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum) 302 break; 303 } 304 305 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x); 306 } 307 308 309 /** 310 * @brief Return the exponential integral @f$ Ei(x) @f$ 311 * by asymptotic expansion. 312 * 313 * The exponential integral is given by 314 * \f[ 315 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 316 * \f] 317 * 318 * @param __x The argument of the exponential integral function. 319 * @return The exponential integral. 320 */ 321 template<typename _Tp> 322 _Tp 323 __expint_Ei_asymp(_Tp __x) 324 { 325 _Tp __term = _Tp(1); 326 _Tp __sum = _Tp(1); 327 const unsigned int __max_iter = 1000; 328 for (unsigned int __i = 1; __i < __max_iter; ++__i) 329 { 330 _Tp __prev = __term; 331 __term *= __i / __x; 332 if (__term < std::numeric_limits<_Tp>::epsilon()) 333 break; 334 if (__term >= __prev) 335 break; 336 __sum += __term; 337 } 338 339 return std::exp(__x) * __sum / __x; 340 } 341 342 343 /** 344 * @brief Return the exponential integral @f$ Ei(x) @f$. 345 * 346 * The exponential integral is given by 347 * \f[ 348 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 349 * \f] 350 * 351 * @param __x The argument of the exponential integral function. 352 * @return The exponential integral. 353 */ 354 template<typename _Tp> 355 _Tp 356 __expint_Ei(_Tp __x) 357 { 358 if (__x < _Tp(0)) 359 return -__expint_E1(-__x); 360 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon())) 361 return __expint_Ei_series(__x); 362 else 363 return __expint_Ei_asymp(__x); 364 } 365 366 367 /** 368 * @brief Return the exponential integral @f$ E_1(x) @f$. 369 * 370 * The exponential integral is given by 371 * \f[ 372 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 373 * \f] 374 * 375 * @param __x The argument of the exponential integral function. 376 * @return The exponential integral. 377 */ 378 template<typename _Tp> 379 _Tp 380 __expint_E1(_Tp __x) 381 { 382 if (__x < _Tp(0)) 383 return -__expint_Ei(-__x); 384 else if (__x < _Tp(1)) 385 return __expint_E1_series(__x); 386 else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point. 387 return __expint_En_cont_frac(1, __x); 388 else 389 return __expint_E1_asymp(__x); 390 } 391 392 393 /** 394 * @brief Return the exponential integral @f$ E_n(x) @f$ 395 * for large argument. 396 * 397 * The exponential integral is given by 398 * \f[ 399 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 400 * \f] 401 * 402 * This is something of an extension. 403 * 404 * @param __n The order of the exponential integral function. 405 * @param __x The argument of the exponential integral function. 406 * @return The exponential integral. 407 */ 408 template<typename _Tp> 409 _Tp 410 __expint_asymp(unsigned int __n, _Tp __x) 411 { 412 _Tp __term = _Tp(1); 413 _Tp __sum = _Tp(1); 414 for (unsigned int __i = 1; __i <= __n; ++__i) 415 { 416 _Tp __prev = __term; 417 __term *= -(__n - __i + 1) / __x; 418 if (std::abs(__term) > std::abs(__prev)) 419 break; 420 __sum += __term; 421 } 422 423 return std::exp(-__x) * __sum / __x; 424 } 425 426 427 /** 428 * @brief Return the exponential integral @f$ E_n(x) @f$ 429 * for large order. 430 * 431 * The exponential integral is given by 432 * \f[ 433 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 434 * \f] 435 * 436 * This is something of an extension. 437 * 438 * @param __n The order of the exponential integral function. 439 * @param __x The argument of the exponential integral function. 440 * @return The exponential integral. 441 */ 442 template<typename _Tp> 443 _Tp 444 __expint_large_n(unsigned int __n, _Tp __x) 445 { 446 const _Tp __xpn = __x + __n; 447 const _Tp __xpn2 = __xpn * __xpn; 448 _Tp __term = _Tp(1); 449 _Tp __sum = _Tp(1); 450 for (unsigned int __i = 1; __i <= __n; ++__i) 451 { 452 _Tp __prev = __term; 453 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2; 454 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) 455 break; 456 __sum += __term; 457 } 458 459 return std::exp(-__x) * __sum / __xpn; 460 } 461 462 463 /** 464 * @brief Return the exponential integral @f$ E_n(x) @f$. 465 * 466 * The exponential integral is given by 467 * \f[ 468 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 469 * \f] 470 * This is something of an extension. 471 * 472 * @param __n The order of the exponential integral function. 473 * @param __x The argument of the exponential integral function. 474 * @return The exponential integral. 475 */ 476 template<typename _Tp> 477 _Tp 478 __expint(unsigned int __n, _Tp __x) 479 { 480 // Return NaN on NaN input. 481 if (__isnan(__x)) 482 return std::numeric_limits<_Tp>::quiet_NaN(); 483 else if (__n <= 1 && __x == _Tp(0)) 484 return std::numeric_limits<_Tp>::infinity(); 485 else 486 { 487 _Tp __E0 = std::exp(__x) / __x; 488 if (__n == 0) 489 return __E0; 490 491 _Tp __E1 = __expint_E1(__x); 492 if (__n == 1) 493 return __E1; 494 495 if (__x == _Tp(0)) 496 return _Tp(1) / static_cast<_Tp>(__n - 1); 497 498 _Tp __En = __expint_En_recursion(__n, __x); 499 500 return __En; 501 } 502 } 503 504 505 /** 506 * @brief Return the exponential integral @f$ Ei(x) @f$. 507 * 508 * The exponential integral is given by 509 * \f[ 510 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 511 * \f] 512 * 513 * @param __x The argument of the exponential integral function. 514 * @return The exponential integral. 515 */ 516 template<typename _Tp> 517 inline _Tp 518 __expint(_Tp __x) 519 { 520 if (__isnan(__x)) 521 return std::numeric_limits<_Tp>::quiet_NaN(); 522 else 523 return __expint_Ei(__x); 524 } 525 } // namespace __detail 526 #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) 527 } // namespace tr1 528 #endif 529 530 _GLIBCXX_END_NAMESPACE_VERSION 531 } 532 533 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC 534