Subversion Repositories Kolibri OS

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

// Copyright (C) 2006-2013 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(_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(_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(_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(unsigned int __n, _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(unsigned int __n, _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(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

    _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

    _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

    _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

    _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

    _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

    _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

    _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

    inline _Tp

    __expint(_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