11.2.0/tr1/exp_integral.tcc

// Special functions -*- C++ -*-

// Copyright (C) 2006-2021 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 <tr1/special_function_util.h>

namespace std _GLIBCXX_VISIBILITY(default)
{
_GLIBCXX_BEGIN_NAMESPACE_VERSION

#if _GLIBCXX_USE_STD_SPEC_FUNCS
#elif defined(_GLIBCXX_TR1_CMATH)
namespace tr1
{
#else
# error do not include this header directly, use <cmath> or <tr1/cmath>
#endif
  // [5.2] Special functions

  // Implementation-space details.
  namespace __detail
  {
    template<typename _Tp> _Tp __expint_E1(_Tp);

    /**
     *   @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(_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 = 1000;
      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(_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(unsigned int __n, _Tp __x)
    {
      const unsigned int __max_iter = 1000;
      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(unsigned int __n, _Tp __x)
    {
      const unsigned int __max_iter = 1000;
      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(unsigned int __n, _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(_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(_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(_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(_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(unsigned int __n, _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(unsigned int __n, _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(unsigned int __n, _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(_Tp __x)
    {
      if (__isnan(__x))
        return std::numeric_limits<_Tp>::quiet_NaN();
      else
        return __expint_Ei(__x);
    }
  } // namespace __detail
#if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)
} // namespace tr1
#endif

_GLIBCXX_END_NAMESPACE_VERSION
}

#endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC
Revision:	1166
Committed:	Tue Oct 26 14:22:36 2021 UTC (4 years ago) by rossy
File size:	16013 byte(s)
Log Message:	Daodan: Replace MinGW build env with an up-to-date MSYS2 env
#	Content
1	// Special functions -- C++ --
2
3	// Copyright (C) 2006-2021 Free Software Foundation, Inc.
4	//
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)
9	// any later version.
10	//
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.
15	//
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.
19
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/>.
24
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}
28	*/
29
30	//
31	// ISO C++ 14882 TR1: 5.2 Special functions
32	//
33
34	// Written by Edward Smith-Rowland based on:
35	//
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.
43	//
44
45	#ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC
46	#define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1
47
48	#include <tr1/special_function_util.h>
49
50	namespace std _GLIBCXX_VISIBILITY(default)
51	{
52	_GLIBCXX_BEGIN_NAMESPACE_VERSION
53
54	#if _GLIBCXX_USE_STD_SPEC_FUNCS
55	#elif defined(_GLIBCXX_TR1_CMATH)
56	namespace tr1
57	{
58	#else
59	# error do not include this header directly, use <cmath> or <tr1/cmath>
60	#endif
61	// [5.2] Special functions
62
63	// Implementation-space details.
64	namespace __detail
65	{
66	template<typename _Tp> _Tp __expint_E1(_Tp);
67
68	/**
69	* @brief Return the exponential integral @f$ E_1(x) @f$
70	* by series summation. This should be good
71	* for @f$ x < 1 @f$.
72	*
73	* The exponential integral is given by
74	* \f[
75	* E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
76	* \f]
77	*
78	* @param __x The argument of the exponential integral function.
79	* @return The exponential integral.
80	*/
81	template<typename _Tp>
82	_Tp
83	__expint_E1_series(_Tp __x)
84	{
85	const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
86	_Tp __term = _Tp(1);
87	_Tp __esum = _Tp(0);
88	_Tp __osum = _Tp(0);
89	const unsigned int __max_iter = 1000;
90	for (unsigned int __i = 1; __i < __max_iter; ++__i)
91	{
92	__term *= - __x / __i;
93	if (std::abs(__term) < __eps)
94	break;
95	if (__term >= _Tp(0))
96	__esum += __term / __i;
97	else
98	__osum += __term / __i;
99	}
100
101	return - __esum - __osum
102	- __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
103	}
104
105
106	/**
107	* @brief Return the exponential integral @f$ E_1(x) @f$
108	* by asymptotic expansion.
109	*
110	* The exponential integral is given by
111	* \f[
112	* E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
113	* \f]
114	*
115	* @param __x The argument of the exponential integral function.
116	* @return The exponential integral.
117	*/
118	template<typename _Tp>
119	_Tp
120	__expint_E1_asymp(_Tp __x)
121	{
122	_Tp __term = _Tp(1);
123	_Tp __esum = _Tp(1);
124	_Tp __osum = _Tp(0);
125	const unsigned int __max_iter = 1000;
126	for (unsigned int __i = 1; __i < __max_iter; ++__i)
127	{
128	_Tp __prev = __term;
129	__term *= - __i / __x;
130	if (std::abs(__term) > std::abs(__prev))
131	break;
132	if (__term >= _Tp(0))
133	__esum += __term;
134	else
135	__osum += __term;
136	}
137
138	return std::exp(- __x) * (__esum + __osum) / __x;
139	}
140
141
142	/**
143	* @brief Return the exponential integral @f$ E_n(x) @f$
144	* by series summation.
145	*
146	* The exponential integral is given by
147	* \f[
148	* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
149	* \f]
150	*
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.
154	*/
155	template<typename _Tp>
156	_Tp
157	__expint_En_series(unsigned int __n, _Tp __x)
158	{
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());
165	_Tp __fact = _Tp(1);
166	for (int __i = 1; __i <= __max_iter; ++__i)
167	{
168	__fact *= -__x / _Tp(__i);
169	_Tp __del;
170	if ( __i != __nm1 )
171	__del = -__fact / _Tp(__i - __nm1);
172	else
173	{
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));
178	}
179	__ans += __del;
180	if (std::abs(__del) < __eps * std::abs(__ans))
181	return __ans;
182	}
183	std::__throw_runtime_error(__N("Series summation failed "
184	"in __expint_En_series."));
185	}
186
187
188	/**
189	* @brief Return the exponential integral @f$ E_n(x) @f$
190	* by continued fractions.
191	*
192	* The exponential integral is given by
193	* \f[
194	* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
195	* \f]
196	*
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.
200	*/
201	template<typename _Tp>
202	_Tp
203	__expint_En_cont_frac(unsigned int __n, _Tp __x)
204	{
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;
212	_Tp __h = __d;
213	for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
214	{
215	_Tp __a = -_Tp(__i * (__nm1 + __i));
216	__b += _Tp(2);
217	__d = _Tp(1) / (__a * __d + __b);
218	__c = __b + __a / __c;
219	const _Tp __del = __c * __d;
220	__h *= __del;
221	if (std::abs(__del - _Tp(1)) < __eps)
222	{
223	const _Tp __ans = __h * std::exp(-__x);
224	return __ans;
225	}
226	}
227	std::__throw_runtime_error(__N("Continued fraction failed "
228	"in __expint_En_cont_frac."));
229	}
230
231
232	/**
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.
236	*
237	* The exponential integral is given by
238	* \f[
239	* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
240	* \f]
241	*
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.
245	*/
246	template<typename _Tp>
247	_Tp
248	__expint_En_recursion(unsigned int __n, _Tp __x)
249	{
250	_Tp __En;
251	_Tp __E1 = __expint_E1(__x);
252	if (__x < _Tp(__n))
253	{
254	// Forward recursion is stable only for n < x.
255	__En = __E1;
256	for (unsigned int __j = 2; __j < __n; ++__j)
257	__En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
258	}
259	else
260	{
261	// Backward recursion is stable only for n >= x.
262	__En = _Tp(1);
263	const int __N = __n + 20; // TODO: Check this starting number.
264	_Tp __save = _Tp(0);
265	for (int __j = __N; __j > 0; --__j)
266	{
267	__En = (std::exp(-__x) - __j * __En) / __x;
268	if (__j == __n)
269	__save = __En;
270	}
271	_Tp __norm = __En / __E1;
272	__En /= __norm;
273	}
274
275	return __En;
276	}
277
278	/**
279	* @brief Return the exponential integral @f$ Ei(x) @f$
280	* by series summation.
281	*
282	* The exponential integral is given by
283	* \f[
284	* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
285	* \f]
286	*
287	* @param __x The argument of the exponential integral function.
288	* @return The exponential integral.
289	*/
290	template<typename _Tp>
291	_Tp
292	__expint_Ei_series(_Tp __x)
293	{
294	_Tp __term = _Tp(1);
295	_Tp __sum = _Tp(0);
296	const unsigned int __max_iter = 1000;
297	for (unsigned int __i = 1; __i < __max_iter; ++__i)
298	{
299	__term *= __x / __i;
300	__sum += __term / __i;
301	if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
302	break;
303	}
304
305	return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
306	}
307
308
309	/**
310	* @brief Return the exponential integral @f$ Ei(x) @f$
311	* by asymptotic expansion.
312	*
313	* The exponential integral is given by
314	* \f[
315	* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
316	* \f]
317	*
318	* @param __x The argument of the exponential integral function.
319	* @return The exponential integral.
320	*/
321	template<typename _Tp>
322	_Tp
323	__expint_Ei_asymp(_Tp __x)
324	{
325	_Tp __term = _Tp(1);
326	_Tp __sum = _Tp(1);
327	const unsigned int __max_iter = 1000;
328	for (unsigned int __i = 1; __i < __max_iter; ++__i)
329	{
330	_Tp __prev = __term;
331	__term *= __i / __x;
332	if (__term < std::numeric_limits<_Tp>::epsilon())
333	break;
334	if (__term >= __prev)
335	break;
336	__sum += __term;
337	}
338
339	return std::exp(__x) * __sum / __x;
340	}
341
342
343	/**
344	* @brief Return the exponential integral @f$ Ei(x) @f$.
345	*
346	* The exponential integral is given by
347	* \f[
348	* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
349	* \f]
350	*
351	* @param __x The argument of the exponential integral function.
352	* @return The exponential integral.
353	*/
354	template<typename _Tp>
355	_Tp
356	__expint_Ei(_Tp __x)
357	{
358	if (__x < _Tp(0))
359	return -__expint_E1(-__x);
360	else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
361	return __expint_Ei_series(__x);
362	else
363	return __expint_Ei_asymp(__x);
364	}
365
366
367	/**
368	* @brief Return the exponential integral @f$ E_1(x) @f$.
369	*
370	* The exponential integral is given by
371	* \f[
372	* E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
373	* \f]
374	*
375	* @param __x The argument of the exponential integral function.
376	* @return The exponential integral.
377	*/
378	template<typename _Tp>
379	_Tp
380	__expint_E1(_Tp __x)
381	{
382	if (__x < _Tp(0))
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);
388	else
389	return __expint_E1_asymp(__x);
390	}
391
392
393	/**
394	* @brief Return the exponential integral @f$ E_n(x) @f$
395	* for large argument.
396	*
397	* The exponential integral is given by
398	* \f[
399	* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
400	* \f]
401	*
402	* This is something of an extension.
403	*
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.
407	*/
408	template<typename _Tp>
409	_Tp
410	__expint_asymp(unsigned int __n, _Tp __x)
411	{
412	_Tp __term = _Tp(1);
413	_Tp __sum = _Tp(1);
414	for (unsigned int __i = 1; __i <= __n; ++__i)
415	{
416	_Tp __prev = __term;
417	__term *= -(__n - __i + 1) / __x;
418	if (std::abs(__term) > std::abs(__prev))
419	break;
420	__sum += __term;
421	}
422
423	return std::exp(-__x) * __sum / __x;
424	}
425
426
427	/**
428	* @brief Return the exponential integral @f$ E_n(x) @f$
429	* for large order.
430	*
431	* The exponential integral is given by
432	* \f[
433	* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
434	* \f]
435	*
436	* This is something of an extension.
437	*
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.
441	*/
442	template<typename _Tp>
443	_Tp
444	__expint_large_n(unsigned int __n, _Tp __x)
445	{
446	const _Tp __xpn = __x + __n;
447	const _Tp __xpn2 = __xpn * __xpn;
448	_Tp __term = _Tp(1);
449	_Tp __sum = _Tp(1);
450	for (unsigned int __i = 1; __i <= __n; ++__i)
451	{
452	_Tp __prev = __term;
453	__term = (__n - 2 (__i - 1) * __x) / __xpn2;
454	if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
455	break;
456	__sum += __term;
457	}
458
459	return std::exp(-__x) * __sum / __xpn;
460	}
461
462
463	/**
464	* @brief Return the exponential integral @f$ E_n(x) @f$.
465	*
466	* The exponential integral is given by
467	* \f[
468	* E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
469	* \f]
470	* This is something of an extension.
471	*
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.
475	*/
476	template<typename _Tp>
477	_Tp
478	__expint(unsigned int __n, _Tp __x)
479	{
480	// Return NaN on NaN input.
481	if (__isnan(__x))
482	return std::numeric_limits<_Tp>::quiet_NaN();
483	else if (__n <= 1 && __x == _Tp(0))
484	return std::numeric_limits<_Tp>::infinity();
485	else
486	{
487	_Tp __E0 = std::exp(__x) / __x;
488	if (__n == 0)
489	return __E0;
490
491	_Tp __E1 = __expint_E1(__x);
492	if (__n == 1)
493	return __E1;
494
495	if (__x == _Tp(0))
496	return _Tp(1) / static_cast<_Tp>(__n - 1);
497
498	_Tp __En = __expint_En_recursion(__n, __x);
499
500	return __En;
501	}
502	}
503
504
505	/**
506	* @brief Return the exponential integral @f$ Ei(x) @f$.
507	*
508	* The exponential integral is given by
509	* \f[
510	* Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
511	* \f]
512	*
513	* @param __x The argument of the exponential integral function.
514	* @return The exponential integral.
515	*/
516	template<typename _Tp>
517	inline _Tp
518	__expint(_Tp __x)
519	{
520	if (__isnan(__x))
521	return std::numeric_limits<_Tp>::quiet_NaN();
522	else
523	return __expint_Ei(__x);
524	}
525	} // namespace __detail
526	#if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)
527	} // namespace tr1
528	#endif
529
530	_GLIBCXX_END_NAMESPACE_VERSION
531	}
532
533	#endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC