1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006, 2007, 2008, 2009 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 * You should not attempt to use it directly. 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 52 { 53 namespace tr1 54 { 55 56 // [5.2] Special functions 57 58 // Implementation-space details. 59 namespace __detail 60 { 61 62 /** 63 * @brief Return the exponential integral @f$ E_1(x) @f$ 64 * by series summation. This should be good 65 * for @f$ x < 1 @f$. 66 * 67 * The exponential integral is given by 68 * \f[ 69 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt 70 * \f] 71 * 72 * @param __x The argument of the exponential integral function. 73 * @return The exponential integral. 74 */ 75 template<typename _Tp> 76 _Tp 77 __expint_E1_series(const _Tp __x) 78 { 79 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 80 _Tp __term = _Tp(1); 81 _Tp __esum = _Tp(0); 82 _Tp __osum = _Tp(0); 83 const unsigned int __max_iter = 100; 84 for (unsigned int __i = 1; __i < __max_iter; ++__i) 85 { 86 __term *= - __x / __i; 87 if (std::abs(__term) < __eps) 88 break; 89 if (__term >= _Tp(0)) 90 __esum += __term / __i; 91 else 92 __osum += __term / __i; 93 } 94 95 return - __esum - __osum 96 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x); 97 } 98 99 100 /** 101 * @brief Return the exponential integral @f$ E_1(x) @f$ 102 * by asymptotic expansion. 103 * 104 * The exponential integral is given by 105 * \f[ 106 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 107 * \f] 108 * 109 * @param __x The argument of the exponential integral function. 110 * @return The exponential integral. 111 */ 112 template<typename _Tp> 113 _Tp 114 __expint_E1_asymp(const _Tp __x) 115 { 116 _Tp __term = _Tp(1); 117 _Tp __esum = _Tp(1); 118 _Tp __osum = _Tp(0); 119 const unsigned int __max_iter = 1000; 120 for (unsigned int __i = 1; __i < __max_iter; ++__i) 121 { 122 _Tp __prev = __term; 123 __term *= - __i / __x; 124 if (std::abs(__term) > std::abs(__prev)) 125 break; 126 if (__term >= _Tp(0)) 127 __esum += __term; 128 else 129 __osum += __term; 130 } 131 132 return std::exp(- __x) * (__esum + __osum) / __x; 133 } 134 135 136 /** 137 * @brief Return the exponential integral @f$ E_n(x) @f$ 138 * by series summation. 139 * 140 * The exponential integral is given by 141 * \f[ 142 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 143 * \f] 144 * 145 * @param __n The order of the exponential integral function. 146 * @param __x The argument of the exponential integral function. 147 * @return The exponential integral. 148 */ 149 template<typename _Tp> 150 _Tp 151 __expint_En_series(const unsigned int __n, const _Tp __x) 152 { 153 const unsigned int __max_iter = 100; 154 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 155 const int __nm1 = __n - 1; 156 _Tp __ans = (__nm1 != 0 157 ? _Tp(1) / __nm1 : -std::log(__x) 158 - __numeric_constants<_Tp>::__gamma_e()); 159 _Tp __fact = _Tp(1); 160 for (int __i = 1; __i <= __max_iter; ++__i) 161 { 162 __fact *= -__x / _Tp(__i); 163 _Tp __del; 164 if ( __i != __nm1 ) 165 __del = -__fact / _Tp(__i - __nm1); 166 else 167 { 168 _Tp __psi = -_TR1_GAMMA_TCC; 169 for (int __ii = 1; __ii <= __nm1; ++__ii) 170 __psi += _Tp(1) / _Tp(__ii); 171 __del = __fact * (__psi - std::log(__x)); 172 } 173 __ans += __del; 174 if (std::abs(__del) < __eps * std::abs(__ans)) 175 return __ans; 176 } 177 std::__throw_runtime_error(__N("Series summation failed " 178 "in __expint_En_series.")); 179 } 180 181 182 /** 183 * @brief Return the exponential integral @f$ E_n(x) @f$ 184 * by continued fractions. 185 * 186 * The exponential integral is given by 187 * \f[ 188 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 189 * \f] 190 * 191 * @param __n The order of the exponential integral function. 192 * @param __x The argument of the exponential integral function. 193 * @return The exponential integral. 194 */ 195 template<typename _Tp> 196 _Tp 197 __expint_En_cont_frac(const unsigned int __n, const _Tp __x) 198 { 199 const unsigned int __max_iter = 100; 200 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 201 const _Tp __fp_min = std::numeric_limits<_Tp>::min(); 202 const int __nm1 = __n - 1; 203 _Tp __b = __x + _Tp(__n); 204 _Tp __c = _Tp(1) / __fp_min; 205 _Tp __d = _Tp(1) / __b; 206 _Tp __h = __d; 207 for ( unsigned int __i = 1; __i <= __max_iter; ++__i ) 208 { 209 _Tp __a = -_Tp(__i * (__nm1 + __i)); 210 __b += _Tp(2); 211 __d = _Tp(1) / (__a * __d + __b); 212 __c = __b + __a / __c; 213 const _Tp __del = __c * __d; 214 __h *= __del; 215 if (std::abs(__del - _Tp(1)) < __eps) 216 { 217 const _Tp __ans = __h * std::exp(-__x); 218 return __ans; 219 } 220 } 221 std::__throw_runtime_error(__N("Continued fraction failed " 222 "in __expint_En_cont_frac.")); 223 } 224 225 226 /** 227 * @brief Return the exponential integral @f$ E_n(x) @f$ 228 * by recursion. Use upward recursion for @f$ x < n @f$ 229 * and downward recursion (Miller's algorithm) otherwise. 230 * 231 * The exponential integral is given by 232 * \f[ 233 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 234 * \f] 235 * 236 * @param __n The order of the exponential integral function. 237 * @param __x The argument of the exponential integral function. 238 * @return The exponential integral. 239 */ 240 template<typename _Tp> 241 _Tp 242 __expint_En_recursion(const unsigned int __n, const _Tp __x) 243 { 244 _Tp __En; 245 _Tp __E1 = __expint_E1(__x); 246 if (__x < _Tp(__n)) 247 { 248 // Forward recursion is stable only for n < x. 249 __En = __E1; 250 for (unsigned int __j = 2; __j < __n; ++__j) 251 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1); 252 } 253 else 254 { 255 // Backward recursion is stable only for n >= x. 256 __En = _Tp(1); 257 const int __N = __n + 20; // TODO: Check this starting number. 258 _Tp __save = _Tp(0); 259 for (int __j = __N; __j > 0; --__j) 260 { 261 __En = (std::exp(-__x) - __j * __En) / __x; 262 if (__j == __n) 263 __save = __En; 264 } 265 _Tp __norm = __En / __E1; 266 __En /= __norm; 267 } 268 269 return __En; 270 } 271 272 /** 273 * @brief Return the exponential integral @f$ Ei(x) @f$ 274 * by series summation. 275 * 276 * The exponential integral is given by 277 * \f[ 278 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 279 * \f] 280 * 281 * @param __x The argument of the exponential integral function. 282 * @return The exponential integral. 283 */ 284 template<typename _Tp> 285 _Tp 286 __expint_Ei_series(const _Tp __x) 287 { 288 _Tp __term = _Tp(1); 289 _Tp __sum = _Tp(0); 290 const unsigned int __max_iter = 1000; 291 for (unsigned int __i = 1; __i < __max_iter; ++__i) 292 { 293 __term *= __x / __i; 294 __sum += __term / __i; 295 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum) 296 break; 297 } 298 299 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x); 300 } 301 302 303 /** 304 * @brief Return the exponential integral @f$ Ei(x) @f$ 305 * by asymptotic expansion. 306 * 307 * The exponential integral is given by 308 * \f[ 309 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 310 * \f] 311 * 312 * @param __x The argument of the exponential integral function. 313 * @return The exponential integral. 314 */ 315 template<typename _Tp> 316 _Tp 317 __expint_Ei_asymp(const _Tp __x) 318 { 319 _Tp __term = _Tp(1); 320 _Tp __sum = _Tp(1); 321 const unsigned int __max_iter = 1000; 322 for (unsigned int __i = 1; __i < __max_iter; ++__i) 323 { 324 _Tp __prev = __term; 325 __term *= __i / __x; 326 if (__term < std::numeric_limits<_Tp>::epsilon()) 327 break; 328 if (__term >= __prev) 329 break; 330 __sum += __term; 331 } 332 333 return std::exp(__x) * __sum / __x; 334 } 335 336 337 /** 338 * @brief Return the exponential integral @f$ Ei(x) @f$. 339 * 340 * The exponential integral is given by 341 * \f[ 342 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 343 * \f] 344 * 345 * @param __x The argument of the exponential integral function. 346 * @return The exponential integral. 347 */ 348 template<typename _Tp> 349 _Tp 350 __expint_Ei(const _Tp __x) 351 { 352 if (__x < _Tp(0)) 353 return -__expint_E1(-__x); 354 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon())) 355 return __expint_Ei_series(__x); 356 else 357 return __expint_Ei_asymp(__x); 358 } 359 360 361 /** 362 * @brief Return the exponential integral @f$ E_1(x) @f$. 363 * 364 * The exponential integral is given by 365 * \f[ 366 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 367 * \f] 368 * 369 * @param __x The argument of the exponential integral function. 370 * @return The exponential integral. 371 */ 372 template<typename _Tp> 373 _Tp 374 __expint_E1(const _Tp __x) 375 { 376 if (__x < _Tp(0)) 377 return -__expint_Ei(-__x); 378 else if (__x < _Tp(1)) 379 return __expint_E1_series(__x); 380 else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point. 381 return __expint_En_cont_frac(1, __x); 382 else 383 return __expint_E1_asymp(__x); 384 } 385 386 387 /** 388 * @brief Return the exponential integral @f$ E_n(x) @f$ 389 * for large argument. 390 * 391 * The exponential integral is given by 392 * \f[ 393 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 394 * \f] 395 * 396 * This is something of an extension. 397 * 398 * @param __n The order of the exponential integral function. 399 * @param __x The argument of the exponential integral function. 400 * @return The exponential integral. 401 */ 402 template<typename _Tp> 403 _Tp 404 __expint_asymp(const unsigned int __n, const _Tp __x) 405 { 406 _Tp __term = _Tp(1); 407 _Tp __sum = _Tp(1); 408 for (unsigned int __i = 1; __i <= __n; ++__i) 409 { 410 _Tp __prev = __term; 411 __term *= -(__n - __i + 1) / __x; 412 if (std::abs(__term) > std::abs(__prev)) 413 break; 414 __sum += __term; 415 } 416 417 return std::exp(-__x) * __sum / __x; 418 } 419 420 421 /** 422 * @brief Return the exponential integral @f$ E_n(x) @f$ 423 * for large order. 424 * 425 * The exponential integral is given by 426 * \f[ 427 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 428 * \f] 429 * 430 * This is something of an extension. 431 * 432 * @param __n The order of the exponential integral function. 433 * @param __x The argument of the exponential integral function. 434 * @return The exponential integral. 435 */ 436 template<typename _Tp> 437 _Tp 438 __expint_large_n(const unsigned int __n, const _Tp __x) 439 { 440 const _Tp __xpn = __x + __n; 441 const _Tp __xpn2 = __xpn * __xpn; 442 _Tp __term = _Tp(1); 443 _Tp __sum = _Tp(1); 444 for (unsigned int __i = 1; __i <= __n; ++__i) 445 { 446 _Tp __prev = __term; 447 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2; 448 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) 449 break; 450 __sum += __term; 451 } 452 453 return std::exp(-__x) * __sum / __xpn; 454 } 455 456 457 /** 458 * @brief Return the exponential integral @f$ E_n(x) @f$. 459 * 460 * The exponential integral is given by 461 * \f[ 462 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 463 * \f] 464 * This is something of an extension. 465 * 466 * @param __n The order of the exponential integral function. 467 * @param __x The argument of the exponential integral function. 468 * @return The exponential integral. 469 */ 470 template<typename _Tp> 471 _Tp 472 __expint(const unsigned int __n, const _Tp __x) 473 { 474 // Return NaN on NaN input. 475 if (__isnan(__x)) 476 return std::numeric_limits<_Tp>::quiet_NaN(); 477 else if (__n <= 1 && __x == _Tp(0)) 478 return std::numeric_limits<_Tp>::infinity(); 479 else 480 { 481 _Tp __E0 = std::exp(__x) / __x; 482 if (__n == 0) 483 return __E0; 484 485 _Tp __E1 = __expint_E1(__x); 486 if (__n == 1) 487 return __E1; 488 489 if (__x == _Tp(0)) 490 return _Tp(1) / static_cast<_Tp>(__n - 1); 491 492 _Tp __En = __expint_En_recursion(__n, __x); 493 494 return __En; 495 } 496 } 497 498 499 /** 500 * @brief Return the exponential integral @f$ Ei(x) @f$. 501 * 502 * The exponential integral is given by 503 * \f[ 504 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 505 * \f] 506 * 507 * @param __x The argument of the exponential integral function. 508 * @return The exponential integral. 509 */ 510 template<typename _Tp> 511 inline _Tp 512 __expint(const _Tp __x) 513 { 514 if (__isnan(__x)) 515 return std::numeric_limits<_Tp>::quiet_NaN(); 516 else 517 return __expint_Ei(__x); 518 } 519 520 } // namespace std::tr1::__detail 521 } 522 } 523 524 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC 525
