1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006-2015 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)
54 // [5.2] Special functions
56 // Implementation-space details.
59 _GLIBCXX_BEGIN_NAMESPACE_VERSION
61 template<typename _Tp> _Tp __expint_E1(_Tp);
64 * @brief Return the exponential integral @f$ E_1(x) @f$
65 * by series summation. This should be good
68 * The exponential integral is given by
70 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
73 * @param __x The argument of the exponential integral function.
74 * @return The exponential integral.
76 template<typename _Tp>
78 __expint_E1_series(_Tp __x)
80 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
84 const unsigned int __max_iter = 100;
85 for (unsigned int __i = 1; __i < __max_iter; ++__i)
87 __term *= - __x / __i;
88 if (std::abs(__term) < __eps)
91 __esum += __term / __i;
93 __osum += __term / __i;
96 return - __esum - __osum
97 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
102 * @brief Return the exponential integral @f$ E_1(x) @f$
103 * by asymptotic expansion.
105 * The exponential integral is given by
107 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
110 * @param __x The argument of the exponential integral function.
111 * @return The exponential integral.
113 template<typename _Tp>
115 __expint_E1_asymp(_Tp __x)
120 const unsigned int __max_iter = 1000;
121 for (unsigned int __i = 1; __i < __max_iter; ++__i)
124 __term *= - __i / __x;
125 if (std::abs(__term) > std::abs(__prev))
127 if (__term >= _Tp(0))
133 return std::exp(- __x) * (__esum + __osum) / __x;
138 * @brief Return the exponential integral @f$ E_n(x) @f$
139 * by series summation.
141 * The exponential integral is given by
143 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
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.
150 template<typename _Tp>
152 __expint_En_series(unsigned int __n, _Tp __x)
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());
161 for (int __i = 1; __i <= __max_iter; ++__i)
163 __fact *= -__x / _Tp(__i);
166 __del = -__fact / _Tp(__i - __nm1);
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));
175 if (std::abs(__del) < __eps * std::abs(__ans))
178 std::__throw_runtime_error(__N("Series summation failed "
179 "in __expint_En_series."));
184 * @brief Return the exponential integral @f$ E_n(x) @f$
185 * by continued fractions.
187 * The exponential integral is given by
189 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
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.
196 template<typename _Tp>
198 __expint_En_cont_frac(unsigned int __n, _Tp __x)
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;
208 for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
210 _Tp __a = -_Tp(__i * (__nm1 + __i));
212 __d = _Tp(1) / (__a * __d + __b);
213 __c = __b + __a / __c;
214 const _Tp __del = __c * __d;
216 if (std::abs(__del - _Tp(1)) < __eps)
218 const _Tp __ans = __h * std::exp(-__x);
222 std::__throw_runtime_error(__N("Continued fraction failed "
223 "in __expint_En_cont_frac."));
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.
232 * The exponential integral is given by
234 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
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.
241 template<typename _Tp>
243 __expint_En_recursion(unsigned int __n, _Tp __x)
246 _Tp __E1 = __expint_E1(__x);
249 // Forward recursion is stable only for n < x.
251 for (unsigned int __j = 2; __j < __n; ++__j)
252 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
256 // Backward recursion is stable only for n >= x.
258 const int __N = __n + 20; // TODO: Check this starting number.
260 for (int __j = __N; __j > 0; --__j)
262 __En = (std::exp(-__x) - __j * __En) / __x;
266 _Tp __norm = __En / __E1;
274 * @brief Return the exponential integral @f$ Ei(x) @f$
275 * by series summation.
277 * The exponential integral is given by
279 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
282 * @param __x The argument of the exponential integral function.
283 * @return The exponential integral.
285 template<typename _Tp>
287 __expint_Ei_series(_Tp __x)
291 const unsigned int __max_iter = 1000;
292 for (unsigned int __i = 1; __i < __max_iter; ++__i)
295 __sum += __term / __i;
296 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
300 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
305 * @brief Return the exponential integral @f$ Ei(x) @f$
306 * by asymptotic expansion.
308 * The exponential integral is given by
310 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
313 * @param __x The argument of the exponential integral function.
314 * @return The exponential integral.
316 template<typename _Tp>
318 __expint_Ei_asymp(_Tp __x)
322 const unsigned int __max_iter = 1000;
323 for (unsigned int __i = 1; __i < __max_iter; ++__i)
327 if (__term < std::numeric_limits<_Tp>::epsilon())
329 if (__term >= __prev)
334 return std::exp(__x) * __sum / __x;
339 * @brief Return the exponential integral @f$ Ei(x) @f$.
341 * The exponential integral is given by
343 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
346 * @param __x The argument of the exponential integral function.
347 * @return The exponential integral.
349 template<typename _Tp>
354 return -__expint_E1(-__x);
355 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
356 return __expint_Ei_series(__x);
358 return __expint_Ei_asymp(__x);
363 * @brief Return the exponential integral @f$ E_1(x) @f$.
365 * The exponential integral is given by
367 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
370 * @param __x The argument of the exponential integral function.
371 * @return The exponential integral.
373 template<typename _Tp>
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);
384 return __expint_E1_asymp(__x);
389 * @brief Return the exponential integral @f$ E_n(x) @f$
390 * for large argument.
392 * The exponential integral is given by
394 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
397 * This is something of an extension.
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.
403 template<typename _Tp>
405 __expint_asymp(unsigned int __n, _Tp __x)
409 for (unsigned int __i = 1; __i <= __n; ++__i)
412 __term *= -(__n - __i + 1) / __x;
413 if (std::abs(__term) > std::abs(__prev))
418 return std::exp(-__x) * __sum / __x;
423 * @brief Return the exponential integral @f$ E_n(x) @f$
426 * The exponential integral is given by
428 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
431 * This is something of an extension.
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.
437 template<typename _Tp>
439 __expint_large_n(unsigned int __n, _Tp __x)
441 const _Tp __xpn = __x + __n;
442 const _Tp __xpn2 = __xpn * __xpn;
445 for (unsigned int __i = 1; __i <= __n; ++__i)
448 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
449 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
454 return std::exp(-__x) * __sum / __xpn;
459 * @brief Return the exponential integral @f$ E_n(x) @f$.
461 * The exponential integral is given by
463 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
465 * This is something of an extension.
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.
471 template<typename _Tp>
473 __expint(unsigned int __n, _Tp __x)
475 // Return NaN on NaN input.
477 return std::numeric_limits<_Tp>::quiet_NaN();
478 else if (__n <= 1 && __x == _Tp(0))
479 return std::numeric_limits<_Tp>::infinity();
482 _Tp __E0 = std::exp(__x) / __x;
486 _Tp __E1 = __expint_E1(__x);
491 return _Tp(1) / static_cast<_Tp>(__n - 1);
493 _Tp __En = __expint_En_recursion(__n, __x);
501 * @brief Return the exponential integral @f$ Ei(x) @f$.
503 * The exponential integral is given by
505 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
508 * @param __x The argument of the exponential integral function.
509 * @return The exponential integral.
511 template<typename _Tp>
516 return std::numeric_limits<_Tp>::quiet_NaN();
518 return __expint_Ei(__x);
521 _GLIBCXX_END_NAMESPACE_VERSION
522 } // namespace std::tr1::__detail
526 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC