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