lgamma_small.hpp

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

//  (C) Copyright John Maddock 2006.
//  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_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL
#define BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL

#ifdef _MSC_VER
#pragma once
#endif

namespace boost{ namespace math{ namespace detail{

//
// lgamma for small arguments:
//
template <class T, class Policy, class L>
T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l */, const L&)
{
   // This version uses rational approximations for small
   // values of z accurate enough for 64-bit mantissas
   // (80-bit long doubles), works well for 53-bit doubles as well.
   // L is only used to select the Lanczos function.

   BOOST_MATH_STD_USING  // for ADL of std names
   T result = 0;
   if(z < tools::epsilon<T>())
   {
      result = -log(z);
   }
   else if((zm1 == 0) || (zm2 == 0))
   {
      // nothing to do, result is zero....
   }
   else if(z > 2)
   {
      //
      // Begin by performing argument reduction until
      // z is in [2,3):
      //
      if(z >= 3)
      {
         do
         {
            z -= 1;
            zm2 -= 1;
            result += log(z);
         }while(z >= 3);
         // Update zm2, we need it below:
         zm2 = z - 2;
      }

      //
      // Use the following form:
      //
      // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
      //
      // where R(z-2) is a rational approximation optimised for
      // low absolute error - as long as it's absolute error
      // is small compared to the constant Y - then any rounding
      // error in it's computation will get wiped out.
      //
      // R(z-2) has the following properties:
      //
      // At double: Max error found:                    4.231e-18
      // At long double: Max error found:               1.987e-21
      // Maximum Deviation Found (approximation error): 5.900e-24
      //
      static const T P[] = {
         static_cast<T>(-0.180355685678449379109e-1L),
         static_cast<T>(0.25126649619989678683e-1L),
         static_cast<T>(0.494103151567532234274e-1L),
         static_cast<T>(0.172491608709613993966e-1L),
         static_cast<T>(-0.259453563205438108893e-3L),
         static_cast<T>(-0.541009869215204396339e-3L),
         static_cast<T>(-0.324588649825948492091e-4L)
      };
      static const T Q[] = {
         static_cast<T>(0.1e1),
         static_cast<T>(0.196202987197795200688e1L),
         static_cast<T>(0.148019669424231326694e1L),
         static_cast<T>(0.541391432071720958364e0L),
         static_cast<T>(0.988504251128010129477e-1L),
         static_cast<T>(0.82130967464889339326e-2L),
         static_cast<T>(0.224936291922115757597e-3L),
         static_cast<T>(-0.223352763208617092964e-6L)
      };

      static const float Y = 0.158963680267333984375e0f;

      T r = zm2 * (z + 1);
      T R = tools::evaluate_polynomial(P, zm2);
      R /= tools::evaluate_polynomial(Q, zm2);

      result +=  r * Y + r * R;
   }
   else
   {
      //
      // If z is less than 1 use recurrance to shift to
      // z in the interval [1,2]:
      //
      if(z < 1)
      {
         result += -log(z);
         zm2 = zm1;
         zm1 = z;
         z += 1;
      }
      //
      // Two approximations, on for z in [1,1.5] and
      // one for z in [1.5,2]:
      //
      if(z <= 1.5)
      {
         //
         // Use the following form:
         //
         // lgamma(z) = (z-1)(z-2)(Y + R(z-1))
         //
         // where R(z-1) is a rational approximation optimised for
         // low absolute error - as long as it's absolute error
         // is small compared to the constant Y - then any rounding
         // error in it's computation will get wiped out.
         //
         // R(z-1) has the following properties:
         //
         // At double precision: Max error found:                1.230011e-17
         // At 80-bit long double precision:   Max error found:  5.631355e-21
         // Maximum Deviation Found:                             3.139e-021
         // Expected Error Term:                                 3.139e-021

         //
         static const float Y = 0.52815341949462890625f;

         static const T P[] = {
            static_cast<T>(0.490622454069039543534e-1L),
            static_cast<T>(-0.969117530159521214579e-1L),
            static_cast<T>(-0.414983358359495381969e0L),
            static_cast<T>(-0.406567124211938417342e0L),
            static_cast<T>(-0.158413586390692192217e0L),
            static_cast<T>(-0.240149820648571559892e-1L),
            static_cast<T>(-0.100346687696279557415e-2L)
         };
         static const T Q[] = {
            static_cast<T>(0.1e1L),
            static_cast<T>(0.302349829846463038743e1L),
            static_cast<T>(0.348739585360723852576e1L),
            static_cast<T>(0.191415588274426679201e1L),
            static_cast<T>(0.507137738614363510846e0L),
            static_cast<T>(0.577039722690451849648e-1L),
            static_cast<T>(0.195768102601107189171e-2L)
         };

         T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1);
         T prefix = zm1 * zm2;

         result += prefix * Y + prefix * r;
      }
      else
      {
         //
         // Use the following form:
         //
         // lgamma(z) = (2-z)(1-z)(Y + R(2-z))
         //
         // where R(2-z) is a rational approximation optimised for
         // low absolute error - as long as it's absolute error
         // is small compared to the constant Y - then any rounding
         // error in it's computation will get wiped out.
         //
         // R(2-z) has the following properties:
         //
         // At double precision, max error found:              1.797565e-17
         // At 80-bit long double precision, max error found:  9.306419e-21
         // Maximum Deviation Found:                           2.151e-021
         // Expected Error Term:                               2.150e-021
         //
         static const float Y = 0.452017307281494140625f;

         static const T P[] = {
            static_cast<T>(-0.292329721830270012337e-1L), 
            static_cast<T>(0.144216267757192309184e0L),
            static_cast<T>(-0.142440390738631274135e0L),
            static_cast<T>(0.542809694055053558157e-1L),
            static_cast<T>(-0.850535976868336437746e-2L),
            static_cast<T>(0.431171342679297331241e-3L)
         };
         static const T Q[] = {
            static_cast<T>(0.1e1),
            static_cast<T>(-0.150169356054485044494e1L),
            static_cast<T>(0.846973248876495016101e0L),
            static_cast<T>(-0.220095151814995745555e0L),
            static_cast<T>(0.25582797155975869989e-1L),
            static_cast<T>(-0.100666795539143372762e-2L),
            static_cast<T>(-0.827193521891290553639e-6L)
         };
         T r = zm2 * zm1;
         T R = tools::evaluate_polynomial(P, -zm2) / tools::evaluate_polynomial(Q, -zm2);

         result += r * Y + r * R;
      }
   }
   return result;
}
template <class T, class Policy, class L>
T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l */, const L&)
{
   //
   // This version uses rational approximations for small
   // values of z accurate enough for 113-bit mantissas
   // (128-bit long doubles).
   //
   BOOST_MATH_STD_USING  // for ADL of std names
   T result = 0;
   if(z < tools::epsilon<T>())
   {
      result = -log(z);
      BOOST_MATH_INSTRUMENT_CODE(result);
   }
   else if((zm1 == 0) || (zm2 == 0))
   {
      // nothing to do, result is zero....
   }
   else if(z > 2)
   {
      //
      // Begin by performing argument reduction until
      // z is in [2,3):
      //
      if(z >= 3)
      {
         do
         {
            z -= 1;
            result += log(z);
         }while(z >= 3);
         zm2 = z - 2;
      }
      BOOST_MATH_INSTRUMENT_CODE(zm2);
      BOOST_MATH_INSTRUMENT_CODE(z);
      BOOST_MATH_INSTRUMENT_CODE(result);

      //
      // Use the following form:
      //
      // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
      //
      // where R(z-2) is a rational approximation optimised for
      // low absolute error - as long as it's absolute error
      // is small compared to the constant Y - then any rounding
      // error in it's computation will get wiped out.
      //
      // Maximum Deviation Found (approximation error)      3.73e-37

      static const T P[] = {