URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [tags/] [gnu-src/] [gcc-4.5.1/] [gcc-4.5.1-or32-1.0rc4/] [libstdc++-v3/] [include/] [tr1/] [exp_integral.tcc] - Blame information for rev 519

Details | Compare with Previous | View Log


// Special functions -*- C++ -*-
 
// Copyright (C) 2006, 2007, 2008, 2009
// 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
// .
 
/** @file tr1/exp_integral.tcc
 *  This is an internal header file, included by other library headers.
 *  You should not attempt to use it directly.
 */
 
//
// 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
{
namespace tr1
{
 
  // [5.2] Special functions
 
  // Implementation-space details.
  namespace __detail
  {
 
    /**
     *   @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
    _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
    _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
    _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
    _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
    _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
    _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
    _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
    _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
    _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
    _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
    _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
    _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
    inline _Tp
    __expint(const _Tp __x)
    {
      if (__isnan(__x))
        return std::numeric_limits<_Tp>::quiet_NaN();
      else
        return __expint_Ei(__x);
    }
 
  } // namespace std::tr1::__detail
}
}
 
#endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [tags/] [gnu-src/] [gcc-4.5.1/] [gcc-4.5.1-or32-1.0rc4/] [libstdc++-v3/] [include/] [tr1/] [exp_integral.tcc] - Blame information for rev 519

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