libstdc++
|
00001 // Special functions -*- C++ -*- 00002 00003 // Copyright (C) 2006, 2007, 2008, 2009 00004 // Free Software Foundation, Inc. 00005 // 00006 // This file is part of the GNU ISO C++ Library. This library is free 00007 // software; you can redistribute it and/or modify it under the 00008 // terms of the GNU General Public License as published by the 00009 // Free Software Foundation; either version 3, or (at your option) 00010 // any later version. 00011 // 00012 // This library is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU General Public License for more details. 00016 // 00017 // Under Section 7 of GPL version 3, you are granted additional 00018 // permissions described in the GCC Runtime Library Exception, version 00019 // 3.1, as published by the Free Software Foundation. 00020 00021 // You should have received a copy of the GNU General Public License and 00022 // a copy of the GCC Runtime Library Exception along with this program; 00023 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 00024 // <http://www.gnu.org/licenses/>. 00025 00026 /** @file tr1/exp_integral.tcc 00027 * This is an internal header file, included by other library headers. 00028 * You should not attempt to use it directly. 00029 */ 00030 00031 // 00032 // ISO C++ 14882 TR1: 5.2 Special functions 00033 // 00034 00035 // Written by Edward Smith-Rowland based on: 00036 // 00037 // (1) Handbook of Mathematical Functions, 00038 // Ed. by Milton Abramowitz and Irene A. Stegun, 00039 // Dover Publications, New-York, Section 5, pp. 228-251. 00040 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 00041 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, 00042 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), 00043 // 2nd ed, pp. 222-225. 00044 // 00045 00046 #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC 00047 #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1 00048 00049 #include "special_function_util.h" 00050 00051 namespace std 00052 { 00053 namespace tr1 00054 { 00055 00056 // [5.2] Special functions 00057 00058 // Implementation-space details. 00059 namespace __detail 00060 { 00061 00062 /** 00063 * @brief Return the exponential integral @f$ E_1(x) @f$ 00064 * by series summation. This should be good 00065 * for @f$ x < 1 @f$. 00066 * 00067 * The exponential integral is given by 00068 * \f[ 00069 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt 00070 * \f] 00071 * 00072 * @param __x The argument of the exponential integral function. 00073 * @return The exponential integral. 00074 */ 00075 template<typename _Tp> 00076 _Tp 00077 __expint_E1_series(const _Tp __x) 00078 { 00079 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00080 _Tp __term = _Tp(1); 00081 _Tp __esum = _Tp(0); 00082 _Tp __osum = _Tp(0); 00083 const unsigned int __max_iter = 100; 00084 for (unsigned int __i = 1; __i < __max_iter; ++__i) 00085 { 00086 __term *= - __x / __i; 00087 if (std::abs(__term) < __eps) 00088 break; 00089 if (__term >= _Tp(0)) 00090 __esum += __term / __i; 00091 else 00092 __osum += __term / __i; 00093 } 00094 00095 return - __esum - __osum 00096 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x); 00097 } 00098 00099 00100 /** 00101 * @brief Return the exponential integral @f$ E_1(x) @f$ 00102 * by asymptotic expansion. 00103 * 00104 * The exponential integral is given by 00105 * \f[ 00106 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 00107 * \f] 00108 * 00109 * @param __x The argument of the exponential integral function. 00110 * @return The exponential integral. 00111 */ 00112 template<typename _Tp> 00113 _Tp 00114 __expint_E1_asymp(const _Tp __x) 00115 { 00116 _Tp __term = _Tp(1); 00117 _Tp __esum = _Tp(1); 00118 _Tp __osum = _Tp(0); 00119 const unsigned int __max_iter = 1000; 00120 for (unsigned int __i = 1; __i < __max_iter; ++__i) 00121 { 00122 _Tp __prev = __term; 00123 __term *= - __i / __x; 00124 if (std::abs(__term) > std::abs(__prev)) 00125 break; 00126 if (__term >= _Tp(0)) 00127 __esum += __term; 00128 else 00129 __osum += __term; 00130 } 00131 00132 return std::exp(- __x) * (__esum + __osum) / __x; 00133 } 00134 00135 00136 /** 00137 * @brief Return the exponential integral @f$ E_n(x) @f$ 00138 * by series summation. 00139 * 00140 * The exponential integral is given by 00141 * \f[ 00142 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 00143 * \f] 00144 * 00145 * @param __n The order of the exponential integral function. 00146 * @param __x The argument of the exponential integral function. 00147 * @return The exponential integral. 00148 */ 00149 template<typename _Tp> 00150 _Tp 00151 __expint_En_series(const unsigned int __n, const _Tp __x) 00152 { 00153 const unsigned int __max_iter = 100; 00154 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00155 const int __nm1 = __n - 1; 00156 _Tp __ans = (__nm1 != 0 00157 ? _Tp(1) / __nm1 : -std::log(__x) 00158 - __numeric_constants<_Tp>::__gamma_e()); 00159 _Tp __fact = _Tp(1); 00160 for (int __i = 1; __i <= __max_iter; ++__i) 00161 { 00162 __fact *= -__x / _Tp(__i); 00163 _Tp __del; 00164 if ( __i != __nm1 ) 00165 __del = -__fact / _Tp(__i - __nm1); 00166 else 00167 { 00168 _Tp __psi = -_TR1_GAMMA_TCC; 00169 for (int __ii = 1; __ii <= __nm1; ++__ii) 00170 __psi += _Tp(1) / _Tp(__ii); 00171 __del = __fact * (__psi - std::log(__x)); 00172 } 00173 __ans += __del; 00174 if (std::abs(__del) < __eps * std::abs(__ans)) 00175 return __ans; 00176 } 00177 std::__throw_runtime_error(__N("Series summation failed " 00178 "in __expint_En_series.")); 00179 } 00180 00181 00182 /** 00183 * @brief Return the exponential integral @f$ E_n(x) @f$ 00184 * by continued fractions. 00185 * 00186 * The exponential integral is given by 00187 * \f[ 00188 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 00189 * \f] 00190 * 00191 * @param __n The order of the exponential integral function. 00192 * @param __x The argument of the exponential integral function. 00193 * @return The exponential integral. 00194 */ 00195 template<typename _Tp> 00196 _Tp 00197 __expint_En_cont_frac(const unsigned int __n, const _Tp __x) 00198 { 00199 const unsigned int __max_iter = 100; 00200 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00201 const _Tp __fp_min = std::numeric_limits<_Tp>::min(); 00202 const int __nm1 = __n - 1; 00203 _Tp __b = __x + _Tp(__n); 00204 _Tp __c = _Tp(1) / __fp_min; 00205 _Tp __d = _Tp(1) / __b; 00206 _Tp __h = __d; 00207 for ( unsigned int __i = 1; __i <= __max_iter; ++__i ) 00208 { 00209 _Tp __a = -_Tp(__i * (__nm1 + __i)); 00210 __b += _Tp(2); 00211 __d = _Tp(1) / (__a * __d + __b); 00212 __c = __b + __a / __c; 00213 const _Tp __del = __c * __d; 00214 __h *= __del; 00215 if (std::abs(__del - _Tp(1)) < __eps) 00216 { 00217 const _Tp __ans = __h * std::exp(-__x); 00218 return __ans; 00219 } 00220 } 00221 std::__throw_runtime_error(__N("Continued fraction failed " 00222 "in __expint_En_cont_frac.")); 00223 } 00224 00225 00226 /** 00227 * @brief Return the exponential integral @f$ E_n(x) @f$ 00228 * by recursion. Use upward recursion for @f$ x < n @f$ 00229 * and downward recursion (Miller's algorithm) otherwise. 00230 * 00231 * The exponential integral is given by 00232 * \f[ 00233 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 00234 * \f] 00235 * 00236 * @param __n The order of the exponential integral function. 00237 * @param __x The argument of the exponential integral function. 00238 * @return The exponential integral. 00239 */ 00240 template<typename _Tp> 00241 _Tp 00242 __expint_En_recursion(const unsigned int __n, const _Tp __x) 00243 { 00244 _Tp __En; 00245 _Tp __E1 = __expint_E1(__x); 00246 if (__x < _Tp(__n)) 00247 { 00248 // Forward recursion is stable only for n < x. 00249 __En = __E1; 00250 for (unsigned int __j = 2; __j < __n; ++__j) 00251 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1); 00252 } 00253 else 00254 { 00255 // Backward recursion is stable only for n >= x. 00256 __En = _Tp(1); 00257 const int __N = __n + 20; // TODO: Check this starting number. 00258 _Tp __save = _Tp(0); 00259 for (int __j = __N; __j > 0; --__j) 00260 { 00261 __En = (std::exp(-__x) - __j * __En) / __x; 00262 if (__j == __n) 00263 __save = __En; 00264 } 00265 _Tp __norm = __En / __E1; 00266 __En /= __norm; 00267 } 00268 00269 return __En; 00270 } 00271 00272 /** 00273 * @brief Return the exponential integral @f$ Ei(x) @f$ 00274 * by series summation. 00275 * 00276 * The exponential integral is given by 00277 * \f[ 00278 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 00279 * \f] 00280 * 00281 * @param __x The argument of the exponential integral function. 00282 * @return The exponential integral. 00283 */ 00284 template<typename _Tp> 00285 _Tp 00286 __expint_Ei_series(const _Tp __x) 00287 { 00288 _Tp __term = _Tp(1); 00289 _Tp __sum = _Tp(0); 00290 const unsigned int __max_iter = 1000; 00291 for (unsigned int __i = 1; __i < __max_iter; ++__i) 00292 { 00293 __term *= __x / __i; 00294 __sum += __term / __i; 00295 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum) 00296 break; 00297 } 00298 00299 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x); 00300 } 00301 00302 00303 /** 00304 * @brief Return the exponential integral @f$ Ei(x) @f$ 00305 * by asymptotic expansion. 00306 * 00307 * The exponential integral is given by 00308 * \f[ 00309 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 00310 * \f] 00311 * 00312 * @param __x The argument of the exponential integral function. 00313 * @return The exponential integral. 00314 */ 00315 template<typename _Tp> 00316 _Tp 00317 __expint_Ei_asymp(const _Tp __x) 00318 { 00319 _Tp __term = _Tp(1); 00320 _Tp __sum = _Tp(1); 00321 const unsigned int __max_iter = 1000; 00322 for (unsigned int __i = 1; __i < __max_iter; ++__i) 00323 { 00324 _Tp __prev = __term; 00325 __term *= __i / __x; 00326 if (__term < std::numeric_limits<_Tp>::epsilon()) 00327 break; 00328 if (__term >= __prev) 00329 break; 00330 __sum += __term; 00331 } 00332 00333 return std::exp(__x) * __sum / __x; 00334 } 00335 00336 00337 /** 00338 * @brief Return the exponential integral @f$ Ei(x) @f$. 00339 * 00340 * The exponential integral is given by 00341 * \f[ 00342 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 00343 * \f] 00344 * 00345 * @param __x The argument of the exponential integral function. 00346 * @return The exponential integral. 00347 */ 00348 template<typename _Tp> 00349 _Tp 00350 __expint_Ei(const _Tp __x) 00351 { 00352 if (__x < _Tp(0)) 00353 return -__expint_E1(-__x); 00354 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon())) 00355 return __expint_Ei_series(__x); 00356 else 00357 return __expint_Ei_asymp(__x); 00358 } 00359 00360 00361 /** 00362 * @brief Return the exponential integral @f$ E_1(x) @f$. 00363 * 00364 * The exponential integral is given by 00365 * \f[ 00366 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 00367 * \f] 00368 * 00369 * @param __x The argument of the exponential integral function. 00370 * @return The exponential integral. 00371 */ 00372 template<typename _Tp> 00373 _Tp 00374 __expint_E1(const _Tp __x) 00375 { 00376 if (__x < _Tp(0)) 00377 return -__expint_Ei(-__x); 00378 else if (__x < _Tp(1)) 00379 return __expint_E1_series(__x); 00380 else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point. 00381 return __expint_En_cont_frac(1, __x); 00382 else 00383 return __expint_E1_asymp(__x); 00384 } 00385 00386 00387 /** 00388 * @brief Return the exponential integral @f$ E_n(x) @f$ 00389 * for large argument. 00390 * 00391 * The exponential integral is given by 00392 * \f[ 00393 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 00394 * \f] 00395 * 00396 * This is something of an extension. 00397 * 00398 * @param __n The order of the exponential integral function. 00399 * @param __x The argument of the exponential integral function. 00400 * @return The exponential integral. 00401 */ 00402 template<typename _Tp> 00403 _Tp 00404 __expint_asymp(const unsigned int __n, const _Tp __x) 00405 { 00406 _Tp __term = _Tp(1); 00407 _Tp __sum = _Tp(1); 00408 for (unsigned int __i = 1; __i <= __n; ++__i) 00409 { 00410 _Tp __prev = __term; 00411 __term *= -(__n - __i + 1) / __x; 00412 if (std::abs(__term) > std::abs(__prev)) 00413 break; 00414 __sum += __term; 00415 } 00416 00417 return std::exp(-__x) * __sum / __x; 00418 } 00419 00420 00421 /** 00422 * @brief Return the exponential integral @f$ E_n(x) @f$ 00423 * for large order. 00424 * 00425 * The exponential integral is given by 00426 * \f[ 00427 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 00428 * \f] 00429 * 00430 * This is something of an extension. 00431 * 00432 * @param __n The order of the exponential integral function. 00433 * @param __x The argument of the exponential integral function. 00434 * @return The exponential integral. 00435 */ 00436 template<typename _Tp> 00437 _Tp 00438 __expint_large_n(const unsigned int __n, const _Tp __x) 00439 { 00440 const _Tp __xpn = __x + __n; 00441 const _Tp __xpn2 = __xpn * __xpn; 00442 _Tp __term = _Tp(1); 00443 _Tp __sum = _Tp(1); 00444 for (unsigned int __i = 1; __i <= __n; ++__i) 00445 { 00446 _Tp __prev = __term; 00447 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2; 00448 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) 00449 break; 00450 __sum += __term; 00451 } 00452 00453 return std::exp(-__x) * __sum / __xpn; 00454 } 00455 00456 00457 /** 00458 * @brief Return the exponential integral @f$ E_n(x) @f$. 00459 * 00460 * The exponential integral is given by 00461 * \f[ 00462 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 00463 * \f] 00464 * This is something of an extension. 00465 * 00466 * @param __n The order of the exponential integral function. 00467 * @param __x The argument of the exponential integral function. 00468 * @return The exponential integral. 00469 */ 00470 template<typename _Tp> 00471 _Tp 00472 __expint(const unsigned int __n, const _Tp __x) 00473 { 00474 // Return NaN on NaN input. 00475 if (__isnan(__x)) 00476 return std::numeric_limits<_Tp>::quiet_NaN(); 00477 else if (__n <= 1 && __x == _Tp(0)) 00478 return std::numeric_limits<_Tp>::infinity(); 00479 else 00480 { 00481 _Tp __E0 = std::exp(__x) / __x; 00482 if (__n == 0) 00483 return __E0; 00484 00485 _Tp __E1 = __expint_E1(__x); 00486 if (__n == 1) 00487 return __E1; 00488 00489 if (__x == _Tp(0)) 00490 return _Tp(1) / static_cast<_Tp>(__n - 1); 00491 00492 _Tp __En = __expint_En_recursion(__n, __x); 00493 00494 return __En; 00495 } 00496 } 00497 00498 00499 /** 00500 * @brief Return the exponential integral @f$ Ei(x) @f$. 00501 * 00502 * The exponential integral is given by 00503 * \f[ 00504 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 00505 * \f] 00506 * 00507 * @param __x The argument of the exponential integral function. 00508 * @return The exponential integral. 00509 */ 00510 template<typename _Tp> 00511 inline _Tp 00512 __expint(const _Tp __x) 00513 { 00514 if (__isnan(__x)) 00515 return std::numeric_limits<_Tp>::quiet_NaN(); 00516 else 00517 return __expint_Ei(__x); 00518 } 00519 00520 } // namespace std::tr1::__detail 00521 } 00522 } 00523 00524 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC