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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libstdc++-v3/] [include/] [tr1/] [exp_integral.tcc] - Blame information for rev 775

Go to most recent revision | Details | Compare with Previous | View Log


// 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
// .
 
/** @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
 
    template _Tp __expint_E1(const _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
    _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);
    }
 
  _GLIBCXX_END_NAMESPACE_VERSION
  } // namespace std::tr1::__detail
}
}
 
#endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libstdc++-v3/] [include/] [tr1/] [exp_integral.tcc] - Blame information for rev 775

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