1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006, 2007, 2008, 2009, 2010
4 // Free Software Foundation, Inc.
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)
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.
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.
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/>.
26 /** @file tr1/exp_integral.tcc
27 * This is an internal header file, included by other library headers.
28 * Do not attempt to use it directly. @headername{tr1/cmath}
32 // ISO C++ 14882 TR1: 5.2 Special functions
35 // Written by Edward Smith-Rowland based on:
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.
46 #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC
47 #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1
49 #include "special_function_util.h"
51 namespace std _GLIBCXX_VISIBILITY(default)
55 // [5.2] Special functions
57 // Implementation-space details.
60 _GLIBCXX_BEGIN_NAMESPACE_VERSION
62 template<typename _Tp> _Tp __expint_E1(const _Tp);
65 * @brief Return the exponential integral @f$ E_1(x) @f$
66 * by series summation. This should be good
69 * The exponential integral is given by
71 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
74 * @param __x The argument of the exponential integral function.
75 * @return The exponential integral.
77 template<typename _Tp>
79 __expint_E1_series(const _Tp __x)
81 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
85 const unsigned int __max_iter = 100;
86 for (unsigned int __i = 1; __i < __max_iter; ++__i)
88 __term *= - __x / __i;
89 if (std::abs(__term) < __eps)
92 __esum += __term / __i;
94 __osum += __term / __i;
97 return - __esum - __osum
98 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
103 * @brief Return the exponential integral @f$ E_1(x) @f$
104 * by asymptotic expansion.
106 * The exponential integral is given by
108 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
111 * @param __x The argument of the exponential integral function.
112 * @return The exponential integral.
114 template<typename _Tp>
116 __expint_E1_asymp(const _Tp __x)
121 const unsigned int __max_iter = 1000;
122 for (unsigned int __i = 1; __i < __max_iter; ++__i)
125 __term *= - __i / __x;
126 if (std::abs(__term) > std::abs(__prev))
128 if (__term >= _Tp(0))
134 return std::exp(- __x) * (__esum + __osum) / __x;
139 * @brief Return the exponential integral @f$ E_n(x) @f$
140 * by series summation.
142 * The exponential integral is given by
144 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
147 * @param __n The order of the exponential integral function.
148 * @param __x The argument of the exponential integral function.
149 * @return The exponential integral.
151 template<typename _Tp>
153 __expint_En_series(const unsigned int __n, const _Tp __x)
155 const unsigned int __max_iter = 100;
156 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
157 const int __nm1 = __n - 1;
158 _Tp __ans = (__nm1 != 0
159 ? _Tp(1) / __nm1 : -std::log(__x)
160 - __numeric_constants<_Tp>::__gamma_e());
162 for (int __i = 1; __i <= __max_iter; ++__i)
164 __fact *= -__x / _Tp(__i);
167 __del = -__fact / _Tp(__i - __nm1);
170 _Tp __psi = -__numeric_constants<_Tp>::gamma_e();
171 for (int __ii = 1; __ii <= __nm1; ++__ii)
172 __psi += _Tp(1) / _Tp(__ii);
173 __del = __fact * (__psi - std::log(__x));
176 if (std::abs(__del) < __eps * std::abs(__ans))
179 std::__throw_runtime_error(__N("Series summation failed "
180 "in __expint_En_series."));
185 * @brief Return the exponential integral @f$ E_n(x) @f$
186 * by continued fractions.
188 * The exponential integral is given by
190 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
193 * @param __n The order of the exponential integral function.
194 * @param __x The argument of the exponential integral function.
195 * @return The exponential integral.
197 template<typename _Tp>
199 __expint_En_cont_frac(const unsigned int __n, const _Tp __x)
201 const unsigned int __max_iter = 100;
202 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
203 const _Tp __fp_min = std::numeric_limits<_Tp>::min();
204 const int __nm1 = __n - 1;
205 _Tp __b = __x + _Tp(__n);
206 _Tp __c = _Tp(1) / __fp_min;
207 _Tp __d = _Tp(1) / __b;
209 for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
211 _Tp __a = -_Tp(__i * (__nm1 + __i));
213 __d = _Tp(1) / (__a * __d + __b);
214 __c = __b + __a / __c;
215 const _Tp __del = __c * __d;
217 if (std::abs(__del - _Tp(1)) < __eps)
219 const _Tp __ans = __h * std::exp(-__x);
223 std::__throw_runtime_error(__N("Continued fraction failed "
224 "in __expint_En_cont_frac."));
229 * @brief Return the exponential integral @f$ E_n(x) @f$
230 * by recursion. Use upward recursion for @f$ x < n @f$
231 * and downward recursion (Miller's algorithm) otherwise.
233 * The exponential integral is given by
235 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
238 * @param __n The order of the exponential integral function.
239 * @param __x The argument of the exponential integral function.
240 * @return The exponential integral.
242 template<typename _Tp>
244 __expint_En_recursion(const unsigned int __n, const _Tp __x)
247 _Tp __E1 = __expint_E1(__x);
250 // Forward recursion is stable only for n < x.
252 for (unsigned int __j = 2; __j < __n; ++__j)
253 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
257 // Backward recursion is stable only for n >= x.
259 const int __N = __n + 20; // TODO: Check this starting number.
261 for (int __j = __N; __j > 0; --__j)
263 __En = (std::exp(-__x) - __j * __En) / __x;
267 _Tp __norm = __En / __E1;
275 * @brief Return the exponential integral @f$ Ei(x) @f$
276 * by series summation.
278 * The exponential integral is given by
280 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
283 * @param __x The argument of the exponential integral function.
284 * @return The exponential integral.
286 template<typename _Tp>
288 __expint_Ei_series(const _Tp __x)
292 const unsigned int __max_iter = 1000;
293 for (unsigned int __i = 1; __i < __max_iter; ++__i)
296 __sum += __term / __i;
297 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
301 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
306 * @brief Return the exponential integral @f$ Ei(x) @f$
307 * by asymptotic expansion.
309 * The exponential integral is given by
311 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
314 * @param __x The argument of the exponential integral function.
315 * @return The exponential integral.
317 template<typename _Tp>
319 __expint_Ei_asymp(const _Tp __x)
323 const unsigned int __max_iter = 1000;
324 for (unsigned int __i = 1; __i < __max_iter; ++__i)
328 if (__term < std::numeric_limits<_Tp>::epsilon())
330 if (__term >= __prev)
335 return std::exp(__x) * __sum / __x;
340 * @brief Return the exponential integral @f$ Ei(x) @f$.
342 * The exponential integral is given by
344 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
347 * @param __x The argument of the exponential integral function.
348 * @return The exponential integral.
350 template<typename _Tp>
352 __expint_Ei(const _Tp __x)
355 return -__expint_E1(-__x);
356 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
357 return __expint_Ei_series(__x);
359 return __expint_Ei_asymp(__x);
364 * @brief Return the exponential integral @f$ E_1(x) @f$.
366 * The exponential integral is given by
368 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
371 * @param __x The argument of the exponential integral function.
372 * @return The exponential integral.
374 template<typename _Tp>
376 __expint_E1(const _Tp __x)
379 return -__expint_Ei(-__x);
380 else if (__x < _Tp(1))
381 return __expint_E1_series(__x);
382 else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point.
383 return __expint_En_cont_frac(1, __x);
385 return __expint_E1_asymp(__x);
390 * @brief Return the exponential integral @f$ E_n(x) @f$
391 * for large argument.
393 * The exponential integral is given by
395 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
398 * This is something of an extension.
400 * @param __n The order of the exponential integral function.
401 * @param __x The argument of the exponential integral function.
402 * @return The exponential integral.
404 template<typename _Tp>
406 __expint_asymp(const unsigned int __n, const _Tp __x)
410 for (unsigned int __i = 1; __i <= __n; ++__i)
413 __term *= -(__n - __i + 1) / __x;
414 if (std::abs(__term) > std::abs(__prev))
419 return std::exp(-__x) * __sum / __x;
424 * @brief Return the exponential integral @f$ E_n(x) @f$
427 * The exponential integral is given by
429 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
432 * This is something of an extension.
434 * @param __n The order of the exponential integral function.
435 * @param __x The argument of the exponential integral function.
436 * @return The exponential integral.
438 template<typename _Tp>
440 __expint_large_n(const unsigned int __n, const _Tp __x)
442 const _Tp __xpn = __x + __n;
443 const _Tp __xpn2 = __xpn * __xpn;
446 for (unsigned int __i = 1; __i <= __n; ++__i)
449 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
450 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
455 return std::exp(-__x) * __sum / __xpn;
460 * @brief Return the exponential integral @f$ E_n(x) @f$.
462 * The exponential integral is given by
464 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
466 * This is something of an extension.
468 * @param __n The order of the exponential integral function.
469 * @param __x The argument of the exponential integral function.
470 * @return The exponential integral.
472 template<typename _Tp>
474 __expint(const unsigned int __n, const _Tp __x)
476 // Return NaN on NaN input.
478 return std::numeric_limits<_Tp>::quiet_NaN();
479 else if (__n <= 1 && __x == _Tp(0))
480 return std::numeric_limits<_Tp>::infinity();
483 _Tp __E0 = std::exp(__x) / __x;
487 _Tp __E1 = __expint_E1(__x);
492 return _Tp(1) / static_cast<_Tp>(__n - 1);
494 _Tp __En = __expint_En_recursion(__n, __x);
502 * @brief Return the exponential integral @f$ Ei(x) @f$.
504 * The exponential integral is given by
506 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
509 * @param __x The argument of the exponential integral function.
510 * @return The exponential integral.
512 template<typename _Tp>
514 __expint(const _Tp __x)
517 return std::numeric_limits<_Tp>::quiet_NaN();
519 return __expint_Ei(__x);
522 _GLIBCXX_END_NAMESPACE_VERSION
523 } // namespace std::tr1::__detail
527 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC