1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006-2017 Free Software Foundation, Inc.
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)
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.
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.
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/>.
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}
31 // ISO C++ 14882 TR1: 5.2 Special functions
34 // Written by Edward Smith-Rowland based on:
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.
45 #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC
46 #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1
48 #include "special_function_util.h"
50 namespace std _GLIBCXX_VISIBILITY(default)
52 _GLIBCXX_BEGIN_NAMESPACE_VERSION
54 #if _GLIBCXX_USE_STD_SPEC_FUNCS
55 #elif defined(_GLIBCXX_TR1_CMATH)
59 # error do not include this header directly, use <cmath> or <tr1/cmath>
61 // [5.2] Special functions
63 // Implementation-space details.
66 template<typename _Tp> _Tp __expint_E1(_Tp);
69 * @brief Return the exponential integral @f$ E_1(x) @f$
70 * by series summation. This should be good
73 * The exponential integral is given by
75 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
78 * @param __x The argument of the exponential integral function.
79 * @return The exponential integral.
81 template<typename _Tp>
83 __expint_E1_series(_Tp __x)
85 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
89 const unsigned int __max_iter = 1000;
90 for (unsigned int __i = 1; __i < __max_iter; ++__i)
92 __term *= - __x / __i;
93 if (std::abs(__term) < __eps)
96 __esum += __term / __i;
98 __osum += __term / __i;
101 return - __esum - __osum
102 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
107 * @brief Return the exponential integral @f$ E_1(x) @f$
108 * by asymptotic expansion.
110 * The exponential integral is given by
112 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
115 * @param __x The argument of the exponential integral function.
116 * @return The exponential integral.
118 template<typename _Tp>
120 __expint_E1_asymp(_Tp __x)
125 const unsigned int __max_iter = 1000;
126 for (unsigned int __i = 1; __i < __max_iter; ++__i)
129 __term *= - __i / __x;
130 if (std::abs(__term) > std::abs(__prev))
132 if (__term >= _Tp(0))
138 return std::exp(- __x) * (__esum + __osum) / __x;
143 * @brief Return the exponential integral @f$ E_n(x) @f$
144 * by series summation.
146 * The exponential integral is given by
148 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
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.
155 template<typename _Tp>
157 __expint_En_series(unsigned int __n, _Tp __x)
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());
166 for (int __i = 1; __i <= __max_iter; ++__i)
168 __fact *= -__x / _Tp(__i);
171 __del = -__fact / _Tp(__i - __nm1);
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));
180 if (std::abs(__del) < __eps * std::abs(__ans))
183 std::__throw_runtime_error(__N("Series summation failed "
184 "in __expint_En_series."));
189 * @brief Return the exponential integral @f$ E_n(x) @f$
190 * by continued fractions.
192 * The exponential integral is given by
194 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
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.
201 template<typename _Tp>
203 __expint_En_cont_frac(unsigned int __n, _Tp __x)
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;
213 for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
215 _Tp __a = -_Tp(__i * (__nm1 + __i));
217 __d = _Tp(1) / (__a * __d + __b);
218 __c = __b + __a / __c;
219 const _Tp __del = __c * __d;
221 if (std::abs(__del - _Tp(1)) < __eps)
223 const _Tp __ans = __h * std::exp(-__x);
227 std::__throw_runtime_error(__N("Continued fraction failed "
228 "in __expint_En_cont_frac."));
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.
237 * The exponential integral is given by
239 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
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.
246 template<typename _Tp>
248 __expint_En_recursion(unsigned int __n, _Tp __x)
251 _Tp __E1 = __expint_E1(__x);
254 // Forward recursion is stable only for n < x.
256 for (unsigned int __j = 2; __j < __n; ++__j)
257 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
261 // Backward recursion is stable only for n >= x.
263 const int __N = __n + 20; // TODO: Check this starting number.
265 for (int __j = __N; __j > 0; --__j)
267 __En = (std::exp(-__x) - __j * __En) / __x;
271 _Tp __norm = __En / __E1;
279 * @brief Return the exponential integral @f$ Ei(x) @f$
280 * by series summation.
282 * The exponential integral is given by
284 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
287 * @param __x The argument of the exponential integral function.
288 * @return The exponential integral.
290 template<typename _Tp>
292 __expint_Ei_series(_Tp __x)
296 const unsigned int __max_iter = 1000;
297 for (unsigned int __i = 1; __i < __max_iter; ++__i)
300 __sum += __term / __i;
301 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
305 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
310 * @brief Return the exponential integral @f$ Ei(x) @f$
311 * by asymptotic expansion.
313 * The exponential integral is given by
315 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
318 * @param __x The argument of the exponential integral function.
319 * @return The exponential integral.
321 template<typename _Tp>
323 __expint_Ei_asymp(_Tp __x)
327 const unsigned int __max_iter = 1000;
328 for (unsigned int __i = 1; __i < __max_iter; ++__i)
332 if (__term < std::numeric_limits<_Tp>::epsilon())
334 if (__term >= __prev)
339 return std::exp(__x) * __sum / __x;
344 * @brief Return the exponential integral @f$ Ei(x) @f$.
346 * The exponential integral is given by
348 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
351 * @param __x The argument of the exponential integral function.
352 * @return The exponential integral.
354 template<typename _Tp>
359 return -__expint_E1(-__x);
360 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
361 return __expint_Ei_series(__x);
363 return __expint_Ei_asymp(__x);
368 * @brief Return the exponential integral @f$ E_1(x) @f$.
370 * The exponential integral is given by
372 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
375 * @param __x The argument of the exponential integral function.
376 * @return The exponential integral.
378 template<typename _Tp>
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);
389 return __expint_E1_asymp(__x);
394 * @brief Return the exponential integral @f$ E_n(x) @f$
395 * for large argument.
397 * The exponential integral is given by
399 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
402 * This is something of an extension.
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.
408 template<typename _Tp>
410 __expint_asymp(unsigned int __n, _Tp __x)
414 for (unsigned int __i = 1; __i <= __n; ++__i)
417 __term *= -(__n - __i + 1) / __x;
418 if (std::abs(__term) > std::abs(__prev))
423 return std::exp(-__x) * __sum / __x;
428 * @brief Return the exponential integral @f$ E_n(x) @f$
431 * The exponential integral is given by
433 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
436 * This is something of an extension.
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.
442 template<typename _Tp>
444 __expint_large_n(unsigned int __n, _Tp __x)
446 const _Tp __xpn = __x + __n;
447 const _Tp __xpn2 = __xpn * __xpn;
450 for (unsigned int __i = 1; __i <= __n; ++__i)
453 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
454 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
459 return std::exp(-__x) * __sum / __xpn;
464 * @brief Return the exponential integral @f$ E_n(x) @f$.
466 * The exponential integral is given by
468 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
470 * This is something of an extension.
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.
476 template<typename _Tp>
478 __expint(unsigned int __n, _Tp __x)
480 // Return NaN on NaN input.
482 return std::numeric_limits<_Tp>::quiet_NaN();
483 else if (__n <= 1 && __x == _Tp(0))
484 return std::numeric_limits<_Tp>::infinity();
487 _Tp __E0 = std::exp(__x) / __x;
491 _Tp __E1 = __expint_E1(__x);
496 return _Tp(1) / static_cast<_Tp>(__n - 1);
498 _Tp __En = __expint_En_recursion(__n, __x);
506 * @brief Return the exponential integral @f$ Ei(x) @f$.
508 * The exponential integral is given by
510 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
513 * @param __x The argument of the exponential integral function.
514 * @return The exponential integral.
516 template<typename _Tp>
521 return std::numeric_limits<_Tp>::quiet_NaN();
523 return __expint_Ei(__x);
525 } // namespace __detail
526 #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)
530 _GLIBCXX_END_NAMESPACE_VERSION
533 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC