526 lines
16 KiB
C++
526 lines
16 KiB
C++
// Special functions -*- C++ -*-
|
|
|
|
// Copyright (C) 2006, 2007, 2008, 2009, 2010
|
|
// Free Software Foundation, Inc.
|
|
//
|
|
// This file is part of the GNU ISO C++ Library. This library is free
|
|
// software; you can redistribute it and/or modify it under the
|
|
// terms of the GNU General Public License as published by the
|
|
// Free Software Foundation; either version 3, or (at your option)
|
|
// any later version.
|
|
//
|
|
// This library is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
// GNU General Public License for more details.
|
|
//
|
|
// Under Section 7 of GPL version 3, you are granted additional
|
|
// permissions described in the GCC Runtime Library Exception, version
|
|
// 3.1, as published by the Free Software Foundation.
|
|
|
|
// You should have received a copy of the GNU General Public License and
|
|
// a copy of the GCC Runtime Library Exception along with this program;
|
|
// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
// <http://www.gnu.org/licenses/>.
|
|
|
|
/** @file tr1/exp_integral.tcc
|
|
* This is an internal header file, included by other library headers.
|
|
* Do not attempt to use it directly. @headername{tr1/cmath}
|
|
*/
|
|
|
|
//
|
|
// ISO C++ 14882 TR1: 5.2 Special functions
|
|
//
|
|
|
|
// Written by Edward Smith-Rowland based on:
|
|
//
|
|
// (1) Handbook of Mathematical Functions,
|
|
// Ed. by Milton Abramowitz and Irene A. Stegun,
|
|
// Dover Publications, New-York, Section 5, pp. 228-251.
|
|
// (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
|
|
// (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,
|
|
// W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),
|
|
// 2nd ed, pp. 222-225.
|
|
//
|
|
|
|
#ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC
|
|
#define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1
|
|
|
|
#include "special_function_util.h"
|
|
|
|
namespace std _GLIBCXX_VISIBILITY(default)
|
|
{
|
|
namespace tr1
|
|
{
|
|
// [5.2] Special functions
|
|
|
|
// Implementation-space details.
|
|
namespace __detail
|
|
{
|
|
_GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_1(x) @f$
|
|
* by series summation. This should be good
|
|
* for @f$ x < 1 @f$.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_E1_series(const _Tp __x)
|
|
{
|
|
const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
|
|
_Tp __term = _Tp(1);
|
|
_Tp __esum = _Tp(0);
|
|
_Tp __osum = _Tp(0);
|
|
const unsigned int __max_iter = 100;
|
|
for (unsigned int __i = 1; __i < __max_iter; ++__i)
|
|
{
|
|
__term *= - __x / __i;
|
|
if (std::abs(__term) < __eps)
|
|
break;
|
|
if (__term >= _Tp(0))
|
|
__esum += __term / __i;
|
|
else
|
|
__osum += __term / __i;
|
|
}
|
|
|
|
return - __esum - __osum
|
|
- __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_1(x) @f$
|
|
* by asymptotic expansion.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_E1_asymp(const _Tp __x)
|
|
{
|
|
_Tp __term = _Tp(1);
|
|
_Tp __esum = _Tp(1);
|
|
_Tp __osum = _Tp(0);
|
|
const unsigned int __max_iter = 1000;
|
|
for (unsigned int __i = 1; __i < __max_iter; ++__i)
|
|
{
|
|
_Tp __prev = __term;
|
|
__term *= - __i / __x;
|
|
if (std::abs(__term) > std::abs(__prev))
|
|
break;
|
|
if (__term >= _Tp(0))
|
|
__esum += __term;
|
|
else
|
|
__osum += __term;
|
|
}
|
|
|
|
return std::exp(- __x) * (__esum + __osum) / __x;
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_n(x) @f$
|
|
* by series summation.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
|
|
* \f]
|
|
*
|
|
* @param __n The order of the exponential integral function.
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_En_series(const unsigned int __n, const _Tp __x)
|
|
{
|
|
const unsigned int __max_iter = 100;
|
|
const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
|
|
const int __nm1 = __n - 1;
|
|
_Tp __ans = (__nm1 != 0
|
|
? _Tp(1) / __nm1 : -std::log(__x)
|
|
- __numeric_constants<_Tp>::__gamma_e());
|
|
_Tp __fact = _Tp(1);
|
|
for (int __i = 1; __i <= __max_iter; ++__i)
|
|
{
|
|
__fact *= -__x / _Tp(__i);
|
|
_Tp __del;
|
|
if ( __i != __nm1 )
|
|
__del = -__fact / _Tp(__i - __nm1);
|
|
else
|
|
{
|
|
_Tp __psi = -__numeric_constants<_Tp>::gamma_e();
|
|
for (int __ii = 1; __ii <= __nm1; ++__ii)
|
|
__psi += _Tp(1) / _Tp(__ii);
|
|
__del = __fact * (__psi - std::log(__x));
|
|
}
|
|
__ans += __del;
|
|
if (std::abs(__del) < __eps * std::abs(__ans))
|
|
return __ans;
|
|
}
|
|
std::__throw_runtime_error(__N("Series summation failed "
|
|
"in __expint_En_series."));
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_n(x) @f$
|
|
* by continued fractions.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
|
|
* \f]
|
|
*
|
|
* @param __n The order of the exponential integral function.
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_En_cont_frac(const unsigned int __n, const _Tp __x)
|
|
{
|
|
const unsigned int __max_iter = 100;
|
|
const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
|
|
const _Tp __fp_min = std::numeric_limits<_Tp>::min();
|
|
const int __nm1 = __n - 1;
|
|
_Tp __b = __x + _Tp(__n);
|
|
_Tp __c = _Tp(1) / __fp_min;
|
|
_Tp __d = _Tp(1) / __b;
|
|
_Tp __h = __d;
|
|
for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
|
|
{
|
|
_Tp __a = -_Tp(__i * (__nm1 + __i));
|
|
__b += _Tp(2);
|
|
__d = _Tp(1) / (__a * __d + __b);
|
|
__c = __b + __a / __c;
|
|
const _Tp __del = __c * __d;
|
|
__h *= __del;
|
|
if (std::abs(__del - _Tp(1)) < __eps)
|
|
{
|
|
const _Tp __ans = __h * std::exp(-__x);
|
|
return __ans;
|
|
}
|
|
}
|
|
std::__throw_runtime_error(__N("Continued fraction failed "
|
|
"in __expint_En_cont_frac."));
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_n(x) @f$
|
|
* by recursion. Use upward recursion for @f$ x < n @f$
|
|
* and downward recursion (Miller's algorithm) otherwise.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
|
|
* \f]
|
|
*
|
|
* @param __n The order of the exponential integral function.
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_En_recursion(const unsigned int __n, const _Tp __x)
|
|
{
|
|
_Tp __En;
|
|
_Tp __E1 = __expint_E1(__x);
|
|
if (__x < _Tp(__n))
|
|
{
|
|
// Forward recursion is stable only for n < x.
|
|
__En = __E1;
|
|
for (unsigned int __j = 2; __j < __n; ++__j)
|
|
__En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
|
|
}
|
|
else
|
|
{
|
|
// Backward recursion is stable only for n >= x.
|
|
__En = _Tp(1);
|
|
const int __N = __n + 20; // TODO: Check this starting number.
|
|
_Tp __save = _Tp(0);
|
|
for (int __j = __N; __j > 0; --__j)
|
|
{
|
|
__En = (std::exp(-__x) - __j * __En) / __x;
|
|
if (__j == __n)
|
|
__save = __En;
|
|
}
|
|
_Tp __norm = __En / __E1;
|
|
__En /= __norm;
|
|
}
|
|
|
|
return __En;
|
|
}
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ Ei(x) @f$
|
|
* by series summation.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_Ei_series(const _Tp __x)
|
|
{
|
|
_Tp __term = _Tp(1);
|
|
_Tp __sum = _Tp(0);
|
|
const unsigned int __max_iter = 1000;
|
|
for (unsigned int __i = 1; __i < __max_iter; ++__i)
|
|
{
|
|
__term *= __x / __i;
|
|
__sum += __term / __i;
|
|
if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
|
|
break;
|
|
}
|
|
|
|
return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ Ei(x) @f$
|
|
* by asymptotic expansion.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_Ei_asymp(const _Tp __x)
|
|
{
|
|
_Tp __term = _Tp(1);
|
|
_Tp __sum = _Tp(1);
|
|
const unsigned int __max_iter = 1000;
|
|
for (unsigned int __i = 1; __i < __max_iter; ++__i)
|
|
{
|
|
_Tp __prev = __term;
|
|
__term *= __i / __x;
|
|
if (__term < std::numeric_limits<_Tp>::epsilon())
|
|
break;
|
|
if (__term >= __prev)
|
|
break;
|
|
__sum += __term;
|
|
}
|
|
|
|
return std::exp(__x) * __sum / __x;
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ Ei(x) @f$.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_Ei(const _Tp __x)
|
|
{
|
|
if (__x < _Tp(0))
|
|
return -__expint_E1(-__x);
|
|
else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
|
|
return __expint_Ei_series(__x);
|
|
else
|
|
return __expint_Ei_asymp(__x);
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_1(x) @f$.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_E1(const _Tp __x)
|
|
{
|
|
if (__x < _Tp(0))
|
|
return -__expint_Ei(-__x);
|
|
else if (__x < _Tp(1))
|
|
return __expint_E1_series(__x);
|
|
else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point.
|
|
return __expint_En_cont_frac(1, __x);
|
|
else
|
|
return __expint_E1_asymp(__x);
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_n(x) @f$
|
|
* for large argument.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
|
|
* \f]
|
|
*
|
|
* This is something of an extension.
|
|
*
|
|
* @param __n The order of the exponential integral function.
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_asymp(const unsigned int __n, const _Tp __x)
|
|
{
|
|
_Tp __term = _Tp(1);
|
|
_Tp __sum = _Tp(1);
|
|
for (unsigned int __i = 1; __i <= __n; ++__i)
|
|
{
|
|
_Tp __prev = __term;
|
|
__term *= -(__n - __i + 1) / __x;
|
|
if (std::abs(__term) > std::abs(__prev))
|
|
break;
|
|
__sum += __term;
|
|
}
|
|
|
|
return std::exp(-__x) * __sum / __x;
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_n(x) @f$
|
|
* for large order.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
|
|
* \f]
|
|
*
|
|
* This is something of an extension.
|
|
*
|
|
* @param __n The order of the exponential integral function.
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint_large_n(const unsigned int __n, const _Tp __x)
|
|
{
|
|
const _Tp __xpn = __x + __n;
|
|
const _Tp __xpn2 = __xpn * __xpn;
|
|
_Tp __term = _Tp(1);
|
|
_Tp __sum = _Tp(1);
|
|
for (unsigned int __i = 1; __i <= __n; ++__i)
|
|
{
|
|
_Tp __prev = __term;
|
|
__term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
|
|
if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
|
|
break;
|
|
__sum += __term;
|
|
}
|
|
|
|
return std::exp(-__x) * __sum / __xpn;
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ E_n(x) @f$.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
|
|
* \f]
|
|
* This is something of an extension.
|
|
*
|
|
* @param __n The order of the exponential integral function.
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
_Tp
|
|
__expint(const unsigned int __n, const _Tp __x)
|
|
{
|
|
// Return NaN on NaN input.
|
|
if (__isnan(__x))
|
|
return std::numeric_limits<_Tp>::quiet_NaN();
|
|
else if (__n <= 1 && __x == _Tp(0))
|
|
return std::numeric_limits<_Tp>::infinity();
|
|
else
|
|
{
|
|
_Tp __E0 = std::exp(__x) / __x;
|
|
if (__n == 0)
|
|
return __E0;
|
|
|
|
_Tp __E1 = __expint_E1(__x);
|
|
if (__n == 1)
|
|
return __E1;
|
|
|
|
if (__x == _Tp(0))
|
|
return _Tp(1) / static_cast<_Tp>(__n - 1);
|
|
|
|
_Tp __En = __expint_En_recursion(__n, __x);
|
|
|
|
return __En;
|
|
}
|
|
}
|
|
|
|
|
|
/**
|
|
* @brief Return the exponential integral @f$ Ei(x) @f$.
|
|
*
|
|
* The exponential integral is given by
|
|
* \f[
|
|
* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
|
|
* \f]
|
|
*
|
|
* @param __x The argument of the exponential integral function.
|
|
* @return The exponential integral.
|
|
*/
|
|
template<typename _Tp>
|
|
inline _Tp
|
|
__expint(const _Tp __x)
|
|
{
|
|
if (__isnan(__x))
|
|
return std::numeric_limits<_Tp>::quiet_NaN();
|
|
else
|
|
return __expint_Ei(__x);
|
|
}
|
|
|
|
_GLIBCXX_END_NAMESPACE_VERSION
|
|
} // namespace std::tr1::__detail
|
|
}
|
|
}
|
|
|
|
#endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC
|