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