gamma.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 1,479 行 · 第 1/4 页

//  Copyright John Maddock 2006-7.
//  Copyright Paul A. Bristow 2007.

//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_SF_GAMMA_HPP
#define BOOST_MATH_SF_GAMMA_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/config.hpp>
#ifdef BOOST_MSVC
# pragma warning(push)
# pragma warning(disable: 4127 4701)
//  // For lexical_cast, until fixed in 1.35?
//  // conditional expression is constant &
//  // Potentially uninitialized local variable 'name' used
#endif
#include <boost/lexical_cast.hpp>
#ifdef BOOST_MSVC
# pragma warning(pop)
#endif
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/fraction.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/special_functions/powm1.hpp>
#include <boost/math/special_functions/sqrt1pm1.hpp>
#include <boost/math/special_functions/lanczos.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/detail/igamma_large.hpp>
#include <boost/math/special_functions/detail/unchecked_factorial.hpp>
#include <boost/math/special_functions/detail/lgamma_small.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/assert.hpp>
#include <boost/mpl/greater.hpp>
#include <boost/mpl/equal_to.hpp>

#include <boost/config/no_tr1/cmath.hpp>
#include <algorithm>

#ifdef BOOST_MATH_INSTRUMENT
#include <iostream>
#include <iomanip>
#include <typeinfo>
#endif

#ifdef BOOST_MSVC
# pragma warning(push)
# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
# pragma warning(disable: 4127) // conditional expression is constant.
# pragma warning(disable: 4100) // unreferenced formal parameter.
// Several variables made comments,
// but some difficulty as whether referenced on not may depend on macro values.
// So to be safe, 4100 warnings suppressed.
// TODO - revisit this?
#endif

namespace boost{ namespace math{

namespace detail{

template <class T>
inline bool is_odd(T v, const boost::true_type&)
{
   int i = static_cast<int>(v);
   return i&1;
}
template <class T>
inline bool is_odd(T v, const boost::false_type&)
{
   // Oh dear can't cast T to int!
   BOOST_MATH_STD_USING
   T modulus = v - 2 * floor(v/2);
   return static_cast<bool>(modulus != 0);
}
template <class T>
inline bool is_odd(T v)
{
   return is_odd(v, ::boost::is_convertible<T, int>());
}

template <class T>
T sinpx(T z)
{
   // Ad hoc function calculates x * sin(pi * x),
   // taking extra care near when x is near a whole number.
   BOOST_MATH_STD_USING
   int sign = 1;
   if(z < 0)
   {
      z = -z;
   }
   else
   {
      sign = -sign;
   }
   T fl = floor(z);
   T dist;
   if(is_odd(fl))
   {
      fl += 1;
      dist = fl - z;
      sign = -sign;
   }
   else
   {
      dist = z - fl;
   }
   BOOST_ASSERT(fl >= 0);
   if(dist > 0.5)
      dist = 1 - dist;
   T result = sin(dist*boost::math::constants::pi<T>());
   return sign*z*result;
} // template <class T> T sinpx(T z)
//
// tgamma(z), with Lanczos support:
//
template <class T, class Policy, class L>
T gamma_imp(T z, const Policy& pol, const L& l)
{
   BOOST_MATH_STD_USING

   T result = 1;

#ifdef BOOST_MATH_INSTRUMENT
   static bool b = false;
   if(!b)
   {
      std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
      b = true;
   }
#endif
   static const char* function = "boost::math::tgamma<%1%>(%1%)";

   if(z <= 0)
   {
      if(floor(z) == z)
         return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
      if(z <= -20)
      {
         result = gamma_imp(-z, pol, l) * sinpx(z);
         if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
            return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
         result = -boost::math::constants::pi<T>() / result;
         if(result == 0)
            return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
         if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
            return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
         return result;
      }

      // shift z to > 1:
      while(z < 0)
      {
         result /= z;
         z += 1;
      }
   }
   if((floor(z) == z) && (z < max_factorial<T>::value))
   {
      result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
   }
   else
   {
      result *= L::lanczos_sum(z);
      if(z * log(z) > tools::log_max_value<T>())
      {
         // we're going to overflow unless this is done with care:
         T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
         if(log(zgh) * z / 2 > tools::log_max_value<T>())
            return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
         T hp = pow(zgh, (z / 2) - T(0.25));
         result *= hp / exp(zgh);
         if(tools::max_value<T>() / hp < result)
            return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
         result *= hp;
      }
      else
      {
         T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
         result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
      }
   }
   return result;
}
//
// lgamma(z) with Lanczos support:
//
template <class T, class Policy, class L>
T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
{
#ifdef BOOST_MATH_INSTRUMENT
   static bool b = false;
   if(!b)
   {
      std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
      b = true;
   }
#endif

   BOOST_MATH_STD_USING

   static const char* function = "boost::math::lgamma<%1%>(%1%)";

   T result = 0;
   int sresult = 1;
   if(z <= 0)
   {
      // reflection formula:
      if(floor(z) == z)
         return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);

      T t = sinpx(z);
      z = -z;
      if(t < 0)
      {
         t = -t;
      }
      else
      {
         sresult = -sresult;
      }
      result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
   }
   else if(z < 15)
   {
      typedef typename policies::precision<T, Policy>::type precision_type;
      typedef typename mpl::if_<
         mpl::less_equal<precision_type, mpl::int_<64> >,
         mpl::int_<64>,
         typename mpl::if_<
            mpl::less_equal<precision_type, mpl::int_<113> >,
            mpl::int_<113>, mpl::int_<0> >::type
          >::type tag_type;
      result = lgamma_small_imp(z, z - 1, z - 2, tag_type(), pol, l);
   }
   else if((z >= 3) && (z < 100))
   {
      // taking the log of tgamma reduces the error, no danger of overflow here:
      result = log(gamma_imp(z, pol, l));
   }
   else
   {
      // regular evaluation:
      T zgh = static_cast<T>(z + L::g() - boost::math::constants::half<T>());
      result = log(zgh) - 1;
      result *= z - 0.5f;
      result += log(L::lanczos_sum_expG_scaled(z));
   }

   if(sign)
      *sign = sresult;
   return result;
}

//
// Incomplete gamma functions follow:
//
template <class T>
struct upper_incomplete_gamma_fract
{
private:
   T z, a;
   int k;
public:
   typedef std::pair<T,T> result_type;

   upper_incomplete_gamma_fract(T a1, T z1)
      : z(z1-a1+1), a(a1), k(0)
   {
   }

   result_type operator()()
   {
      ++k;
      z += 2;
      return result_type(k * (a - k), z);
   }
};

template <class T>
inline T upper_gamma_fraction(T a, T z, int bits)
{
   // Multiply result by z^a * e^-z to get the full
   // upper incomplete integral.  Divide by tgamma(z)
   // to normalise.
   upper_incomplete_gamma_fract<T> f(a, z);
   return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, bits));
}

template <class T>
struct lower_incomplete_gamma_series
{
private:
   T a, z, result;
public:
   typedef T result_type;
   lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}

   T operator()()
   {
      T r = result;
      a += 1;
      result *= z/a;
      return r;
   }
};

template <class T, class Policy>
inline T lower_gamma_series(T a, T z, const Policy& pol)
{
   // Multiply result by ((z^a) * (e^-z) / a) to get the full
   // lower incomplete integral. Then divide by tgamma(a)
   // to get the normalised value.
   lower_incomplete_gamma_series<T> s(a, z);
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
   int bits = policies::digits<T, Policy>();
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
   T zero = 0;
   T result = boost::math::tools::sum_series(s, bits, max_iter, zero);
#else
   T result = boost::math::tools::sum_series(s, bits, max_iter);
#endif
   policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
   return result;
}

//
// Fully generic tgamma and lgamma use the incomplete partial
// sums added together:
//
template <class T, class Policy>
T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
{
   static const char* function = "boost::math::tgamma<%1%>(%1%)";
   BOOST_MATH_STD_USING
   if((z <= 0) && (floor(z) == z))
      return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
   if(z <= -20)
   {
      T result = gamma_imp(-z, pol, l) * sinpx(z);
      if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
         return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
      result = -boost::math::constants::pi<T>() / result;
      if(result == 0)
         return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
      if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
         return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
      return result;
   }
   //
   // The upper gamma fraction is *very* slow for z < 6, actually it's very
   // slow to converge everywhere but recursing until z > 6 gets rid of the
   // worst of it's behaviour.
   //
   T prefix = 1;
   while(z < 6)
   {
      prefix /= z;
      z += 1;