Subversion Repositories altor32

// 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/gamma.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. Milton Abramowitz and Irene A. Stegun,

//       Dover Publications,

//       Section 6, pp. 253-266

//   (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. 213-216

//   (4) Gamma, Exploring Euler's Constant, Julian Havil,

//       Princeton, 2003.

#ifndef _GLIBCXX_TR1_GAMMA_TCC

#define _GLIBCXX_TR1_GAMMA_TCC 1

#include "special_function_util.h"

namespace std _GLIBCXX_VISIBILITY(default)

{

namespace tr1

{

  // Implementation-space details.

  namespace __detail

  {

  _GLIBCXX_BEGIN_NAMESPACE_VERSION

    /**

     *   @brief This returns Bernoulli numbers from a table or by summation

     *          for larger values.

     *

     *   Recursion is unstable.

     *

     *   @param __n the order n of the Bernoulli number.

     *   @return  The Bernoulli number of order n.

     */

    template 

    _Tp __bernoulli_series(unsigned int __n)

    {

      static const _Tp __num[28] = {

        _Tp(1UL),                        -_Tp(1UL) / _Tp(2UL),

        _Tp(1UL) / _Tp(6UL),             _Tp(0UL),

        -_Tp(1UL) / _Tp(30UL),           _Tp(0UL),

        _Tp(1UL) / _Tp(42UL),            _Tp(0UL),

        -_Tp(1UL) / _Tp(30UL),           _Tp(0UL),

        _Tp(5UL) / _Tp(66UL),            _Tp(0UL),

        -_Tp(691UL) / _Tp(2730UL),       _Tp(0UL),

        _Tp(7UL) / _Tp(6UL),             _Tp(0UL),

        -_Tp(3617UL) / _Tp(510UL),       _Tp(0UL),

        _Tp(43867UL) / _Tp(798UL),       _Tp(0UL),

        -_Tp(174611) / _Tp(330UL),       _Tp(0UL),

        _Tp(854513UL) / _Tp(138UL),      _Tp(0UL),

        -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL),

        _Tp(8553103UL) / _Tp(6UL),       _Tp(0UL)

      };

      if (__n == 0)

        return _Tp(1);

      if (__n == 1)

        return -_Tp(1) / _Tp(2);

      //  Take care of the rest of the odd ones.

      if (__n % 2 == 1)

        return _Tp(0);

      //  Take care of some small evens that are painful for the series.

      if (__n < 28)

        return __num[__n];

      _Tp __fact = _Tp(1);

      if ((__n / 2) % 2 == 0)

        __fact *= _Tp(-1);

      for (unsigned int __k = 1; __k <= __n; ++__k)

        __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());

      __fact *= _Tp(2);

      _Tp __sum = _Tp(0);

      for (unsigned int __i = 1; __i < 1000; ++__i)

        {

          _Tp __term = std::pow(_Tp(__i), -_Tp(__n));

          if (__term < std::numeric_limits<_Tp>::epsilon())

            break;

          __sum += __term;

        }

      return __fact * __sum;

    }

    /**

     *   @brief This returns Bernoulli number \f$B_n\f$.

     *

     *   @param __n the order n of the Bernoulli number.

     *   @return  The Bernoulli number of order n.

     */

    template

    inline _Tp

    __bernoulli(const int __n)

    {

      return __bernoulli_series<_Tp>(__n);

    }

    /**

     *   @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion

     *          with Bernoulli number coefficients.  This is like

     *          Sterling's approximation.

     *

     *   @param __x The argument of the log of the gamma function.

     *   @return  The logarithm of the gamma function.

     */

    template

    _Tp

    __log_gamma_bernoulli(const _Tp __x)

    {

      _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x

               + _Tp(0.5L) * std::log(_Tp(2)

               * __numeric_constants<_Tp>::__pi());

      const _Tp __xx = __x * __x;

      _Tp __help = _Tp(1) / __x;

      for ( unsigned int __i = 1; __i < 20; ++__i )

        {

          const _Tp __2i = _Tp(2 * __i);

          __help /= __2i * (__2i - _Tp(1)) * __xx;

          __lg += __bernoulli<_Tp>(2 * __i) * __help;

        }

      return __lg;

    }

    /**

     *   @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method.

     *          This method dominates all others on the positive axis I think.

     *

     *   @param __x The argument of the log of the gamma function.

     *   @return  The logarithm of the gamma function.

     */

    template

    _Tp

    __log_gamma_lanczos(const _Tp __x)

    {

      const _Tp __xm1 = __x - _Tp(1);

      static const _Tp __lanczos_cheb_7[9] = {

       _Tp( 0.99999999999980993227684700473478L),

       _Tp( 676.520368121885098567009190444019L),

       _Tp(-1259.13921672240287047156078755283L),

       _Tp( 771.3234287776530788486528258894L),

       _Tp(-176.61502916214059906584551354L),

       _Tp( 12.507343278686904814458936853L),

       _Tp(-0.13857109526572011689554707L),

       _Tp( 9.984369578019570859563e-6L),

       _Tp( 1.50563273514931155834e-7L)

      };

      static const _Tp __LOGROOT2PI

          = _Tp(0.9189385332046727417803297364056176L);

      _Tp __sum = __lanczos_cheb_7[0];

      for(unsigned int __k = 1; __k < 9; ++__k)

        __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);

      const _Tp __term1 = (__xm1 + _Tp(0.5L))

                        * std::log((__xm1 + _Tp(7.5L))

                       / __numeric_constants<_Tp>::__euler());

      const _Tp __term2 = __LOGROOT2PI + std::log(__sum);

      const _Tp __result = __term1 + (__term2 - _Tp(7));

      return __result;

    }

    /**

     *   @brief Return \f$ log(|\Gamma(x)|) \f$.

     *          This will return values even for \f$ x < 0 \f$.

     *          To recover the sign of \f$ \Gamma(x) \f$ for

     *          any argument use @a __log_gamma_sign.

     *

     *   @param __x The argument of the log of the gamma function.

     *   @return  The logarithm of the gamma function.

     */

    template

    _Tp

    __log_gamma(const _Tp __x)

    {

      if (__x > _Tp(0.5L))

        return __log_gamma_lanczos(__x);

      else

        {

          const _Tp __sin_fact

                 = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x));

          if (__sin_fact == _Tp(0))

            std::__throw_domain_error(__N("Argument is nonpositive integer "

                                          "in __log_gamma"));

          return __numeric_constants<_Tp>::__lnpi()

                     - std::log(__sin_fact)

                     - __log_gamma_lanczos(_Tp(1) - __x);

        }

    }

    /**

     *   @brief Return the sign of \f$ \Gamma(x) \f$.

     *          At nonpositive integers zero is returned.

     *

     *   @param __x The argument of the gamma function.

     *   @return  The sign of the gamma function.

     */

    template

    _Tp

    __log_gamma_sign(const _Tp __x)

    {

      if (__x > _Tp(0))

        return _Tp(1);

      else

        {

          const _Tp __sin_fact

                  = std::sin(__numeric_constants<_Tp>::__pi() * __x);

          if (__sin_fact > _Tp(0))

            return (1);

          else if (__sin_fact < _Tp(0))

            return -_Tp(1);

          else

            return _Tp(0);

        }

    }

    /**

     *   @brief Return the logarithm of the binomial coefficient.

     *   The binomial coefficient is given by:

     *   @f[

     *   \left(  \right) = \frac{n!}{(n-k)! k!}

     *   @f]

     *

     *   @param __n The first argument of the binomial coefficient.

     *   @param __k The second argument of the binomial coefficient.

     *   @return  The binomial coefficient.

     */

    template

    _Tp

    __log_bincoef(const unsigned int __n, const unsigned int __k)

    {

      //  Max e exponent before overflow.

      static const _Tp __max_bincoeff

                      = std::numeric_limits<_Tp>::max_exponent10

                      * std::log(_Tp(10)) - _Tp(1);

#if _GLIBCXX_USE_C99_MATH_TR1

      _Tp __coeff =  std::tr1::lgamma(_Tp(1 + __n))

                  - std::tr1::lgamma(_Tp(1 + __k))

                  - std::tr1::lgamma(_Tp(1 + __n - __k));

#else

      _Tp __coeff =  __log_gamma(_Tp(1 + __n))

                  - __log_gamma(_Tp(1 + __k))

                  - __log_gamma(_Tp(1 + __n - __k));

#endif

    }

    /**

     *   @brief Return the binomial coefficient.

     *   The binomial coefficient is given by:

     *   @f[

     *   \left(  \right) = \frac{n!}{(n-k)! k!}

     *   @f]

     *

     *   @param __n The first argument of the binomial coefficient.

     *   @param __k The second argument of the binomial coefficient.

     *   @return  The binomial coefficient.

     */

    template

    _Tp

    __bincoef(const unsigned int __n, const unsigned int __k)

    {

      //  Max e exponent before overflow.

      static const _Tp __max_bincoeff

                      = std::numeric_limits<_Tp>::max_exponent10

                      * std::log(_Tp(10)) - _Tp(1);

      const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);

      if (__log_coeff > __max_bincoeff)

        return std::numeric_limits<_Tp>::quiet_NaN();

      else

        return std::exp(__log_coeff);

    }

    /**

     *   @brief Return \f$ \Gamma(x) \f$.

     *

     *   @param __x The argument of the gamma function.

     *   @return  The gamma function.

     */

    template

    inline _Tp

    __gamma(const _Tp __x)

    {

      return std::exp(__log_gamma(__x));

    }

    /**

     *   @brief  Return the digamma function by series expansion.

     *   The digamma or @f$ \psi(x) @f$ function is defined by

     *   @f[

     *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}

     *   @f]

     *

     *   The series is given by:

     *   @f[

     *     \psi(x) = -\gamma_E - \frac{1}{x}

     *              \sum_{k=1}^{\infty} \frac{x}{k(x + k)}

     *   @f]

     */

    template

    _Tp

    __psi_series(const _Tp __x)

    {

      _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;

      const unsigned int __max_iter = 100000;

      for (unsigned int __k = 1; __k < __max_iter; ++__k)

        {

          const _Tp __term = __x / (__k * (__k + __x));

          __sum += __term;

          if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())

            break;

        }

      return __sum;

    }

    /**

     *   @brief  Return the digamma function for large argument.

     *   The digamma or @f$ \psi(x) @f$ function is defined by

     *   @f[

     *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}

     *   @f]

     *

     *   The asymptotic series is given by:

     *   @f[

     *     \psi(x) = \ln(x) - \frac{1}{2x}

     *             - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}}

     *   @f]

     */

    template

    _Tp

    __psi_asymp(const _Tp __x)

    {

      _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;

      const _Tp __xx = __x * __x;

      _Tp __xp = __xx;

      const unsigned int __max_iter = 100;

      for (unsigned int __k = 1; __k < __max_iter; ++__k)

        {

          const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);

          __sum -= __term;

          if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())

            break;

          __xp *= __xx;

        }

      return __sum;

    }

    /**

     *   @brief  Return the digamma function.

     *   The digamma or @f$ \psi(x) @f$ function is defined by

     *   @f[

     *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}

     *   @f]

     *   For negative argument the reflection formula is used:

     *   @f[

     *     \psi(x) = \psi(1-x) - \pi \cot(\pi x)

     *   @f]

     */

    template

    _Tp

    __psi(const _Tp __x)

    {

      const int __n = static_cast(__x + 0.5L);

      const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();

      if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps)

        return std::numeric_limits<_Tp>::quiet_NaN();

      else if (__x < _Tp(0))

        {

          const _Tp __pi = __numeric_constants<_Tp>::__pi();

          return __psi(_Tp(1) - __x)

               - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);

        }

      else if (__x > _Tp(100))

        return __psi_asymp(__x);

      else

        return __psi_series(__x);

    }

    /**

     *   @brief  Return the polygamma function @f$ \psi^{(n)}(x) @f$.

     *

     *   The polygamma function is related to the Hurwitz zeta function:

     *   @f[

     *     \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x)

     *   @f]

     */

    template

    _Tp

    __psi(const unsigned int __n, const _Tp __x)

    {

      if (__x <= _Tp(0))

        std::__throw_domain_error(__N("Argument out of range "

                                      "in __psi"));

      else if (__n == 0)

        return __psi(__x);

      else

        {

          const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x);

#if _GLIBCXX_USE_C99_MATH_TR1

          const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));

#else

          const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));

#endif

          _Tp __result = std::exp(__ln_nfact) * __hzeta;

          if (__n % 2 == 1)

            __result = -__result;

          return __result;

        }

    }

  _GLIBCXX_END_NAMESPACE_VERSION

  } // namespace std::tr1::__detail

}

}

#endif // _GLIBCXX_TR1_GAMMA_TCC