lanczos.hpp

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

//  (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_LANCZOS
#define BOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/config.hpp>
#include <boost/mpl/if.hpp>
#include <boost/limits.hpp>
#include <boost/cstdint.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/policies/policy.hpp>
#include <boost/mpl/less_equal.hpp>

#include <limits.h>

namespace boost{ namespace math{ namespace lanczos{

//
// Individual lanczos approximations start here.
//
// Optimal values for G for each N are taken from
// http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf,
// as are the theoretical error bounds.
//
// Constants calculated using the method described by Godfrey
// http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
// http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.
//
// Lanczos Coefficients for N=6 G=5.581
// Max experimental error (with arbitary precision arithmetic) 9.516e-12
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos6 : public mpl::int_<35>
{
   //
   // Produces slightly better than float precision when evaluated at
   // double precision:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[6] = {
         static_cast<T>(8706.349592549009182288174442774377925882L),
         static_cast<T>(8523.650341121874633477483696775067709735L),
         static_cast<T>(3338.029219476423550899999750161289306564L),
         static_cast<T>(653.6424994294008795995653541449610986791L),
         static_cast<T>(63.99951844938187085666201263218840287667L),
         static_cast<T>(2.506628274631006311133031631822390264407L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint16_t) denom[6] = {
         static_cast<boost::uint16_t>(0u),
         static_cast<boost::uint16_t>(24u),
         static_cast<boost::uint16_t>(50u),
         static_cast<boost::uint16_t>(35u),
         static_cast<boost::uint16_t>(10u),
         static_cast<boost::uint16_t>(1u)
      };
      return boost::math::tools::evaluate_rational(num, denom, z);
   }

   template <class T>
   static T lanczos_sum_expG_scaled(const T& z)
   {
      static const T num[6] = {
         static_cast<T>(32.81244541029783471623665933780748627823L),
         static_cast<T>(32.12388941444332003446077108933558534361L),
         static_cast<T>(12.58034729455216106950851080138931470954L),
         static_cast<T>(2.463444478353241423633780693218408889251L),
         static_cast<T>(0.2412010548258800231126240760264822486599L),
         static_cast<T>(0.009446967704539249494420221613134244048319L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint16_t) denom[6] = {
         static_cast<boost::uint16_t>(0u),
         static_cast<boost::uint16_t>(24u),
         static_cast<boost::uint16_t>(50u),
         static_cast<boost::uint16_t>(35u),
         static_cast<boost::uint16_t>(10u),
         static_cast<boost::uint16_t>(1u)
      };
      return boost::math::tools::evaluate_rational(num, denom, z);
   }


   template<class T>
   static T lanczos_sum_near_1(const T& dz)
   {
      static const T d[5] = {
         static_cast<T>(2.044879010930422922760429926121241330235L),
         static_cast<T>(-2.751366405578505366591317846728753993668L),
         static_cast<T>(1.02282965224225004296750609604264824677L),
         static_cast<T>(-0.09786124911582813985028889636665335893627L),
         static_cast<T>(0.0009829742267506615183144364420540766510112L),
      };
      T result = 0;
      for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
      {
         result += (-d[k-1]*dz)/(k*dz + k*k);
      }
      return result;
   }

   template<class T>
   static T lanczos_sum_near_2(const T& dz)
   {
      static const T d[5] = {
         static_cast<T>(5.748142489536043490764289256167080091892L),
         static_cast<T>(-7.734074268282457156081021756682138251825L),
         static_cast<T>(2.875167944990511006997713242805893543947L),
         static_cast<T>(-0.2750873773533504542306766137703788781776L),
         static_cast<T>(0.002763134585812698552178368447708846850353L),
      };
      T result = 0;
      T z = dz + 2;
      for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
      {
         result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
      }
      return result;
   }

   static double g(){ return 5.581000000000000405009359383257105946541; }
};

//
// Lanczos Coefficients for N=11 G=10.900511
// Max experimental error (with arbitary precision arithmetic) 2.16676e-19
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos11 : public mpl::int_<60>
{
   //
   // Produces slightly better than double precision when evaluated at
   // extended-double precision:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[11] = {
         static_cast<T>(38474670393.31776828316099004518914832218L),
         static_cast<T>(36857665043.51950660081971227404959150474L),
         static_cast<T>(15889202453.72942008945006665994637853242L),
         static_cast<T>(4059208354.298834770194507810788393801607L),
         static_cast<T>(680547661.1834733286087695557084801366446L),
         static_cast<T>(78239755.00312005289816041245285376206263L),
         static_cast<T>(6246580.776401795264013335510453568106366L),
         static_cast<T>(341986.3488721347032223777872763188768288L),
         static_cast<T>(12287.19451182455120096222044424100527629L),
         static_cast<T>(261.6140441641668190791708576058805625502L),
         static_cast<T>(2.506628274631000502415573855452633787834L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint32_t) denom[11] = {
         static_cast<boost::uint32_t>(0u),
         static_cast<boost::uint32_t>(362880u),
         static_cast<boost::uint32_t>(1026576u),
         static_cast<boost::uint32_t>(1172700u),
         static_cast<boost::uint32_t>(723680u),
         static_cast<boost::uint32_t>(269325u),
         static_cast<boost::uint32_t>(63273u),
         static_cast<boost::uint32_t>(9450u),
         static_cast<boost::uint32_t>(870u),
         static_cast<boost::uint32_t>(45u),
         static_cast<boost::uint32_t>(1u)
      };
      return boost::math::tools::evaluate_rational(num, denom, z);
   }

   template <class T>
   static T lanczos_sum_expG_scaled(const T& z)
   {
      static const T num[11] = {
         static_cast<T>(709811.662581657956893540610814842699825L),
         static_cast<T>(679979.847415722640161734319823103390728L),
         static_cast<T>(293136.785721159725251629480984140341656L),
         static_cast<T>(74887.5403291467179935942448101441897121L),
         static_cast<T>(12555.29058241386295096255111537516768137L),
         static_cast<T>(1443.42992444170669746078056942194198252L),
         static_cast<T>(115.2419459613734722083208906727972935065L),
         static_cast<T>(6.30923920573262762719523981992008976989L),
         static_cast<T>(0.2266840463022436475495508977579735223818L),
         static_cast<T>(0.004826466289237661857584712046231435101741L),
         static_cast<T>(0.4624429436045378766270459638520555557321e-4L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint32_t) denom[11] = {
         static_cast<boost::uint32_t>(0u),
         static_cast<boost::uint32_t>(362880u),
         static_cast<boost::uint32_t>(1026576u),
         static_cast<boost::uint32_t>(1172700u),
         static_cast<boost::uint32_t>(723680u),
         static_cast<boost::uint32_t>(269325u),
         static_cast<boost::uint32_t>(63273u),
         static_cast<boost::uint32_t>(9450u),
         static_cast<boost::uint32_t>(870u),
         static_cast<boost::uint32_t>(45u),
         static_cast<boost::uint32_t>(1u)
      };
      return boost::math::tools::evaluate_rational(num, denom, z);
   }


   template<class T>
   static T lanczos_sum_near_1(const T& dz)
   {
      static const T d[10] = {
         static_cast<T>(4.005853070677940377969080796551266387954L),
         static_cast<T>(-13.17044315127646469834125159673527183164L),
         static_cast<T>(17.19146865350790353683895137079288129318L),
         static_cast<T>(-11.36446409067666626185701599196274701126L),
         static_cast<T>(4.024801119349323770107694133829772634737L),
         static_cast<T>(-0.7445703262078094128346501724255463005006L),
         static_cast<T>(0.06513861351917497265045550019547857713172L),
         static_cast<T>(-0.00217899958561830354633560009312512312758L),
         static_cast<T>(0.17655204574495137651670832229571934738e-4L),
         static_cast<T>(-0.1036282091079938047775645941885460820853e-7L),
      };
      T result = 0;
      for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
      {
         result += (-d[k-1]*dz)/(k*dz + k*k);
      }
      return result;
   }

   template<class T>
   static T lanczos_sum_near_2(const T& dz)
   {
      static const T d[10] = {
         static_cast<T>(19.05889633808148715159575716844556056056L),
         static_cast<T>(-62.66183664701721716960978577959655644762L),
         static_cast<T>(81.7929198065004751699057192860287512027L),
         static_cast<T>(-54.06941772964234828416072865069196553015L),
         static_cast<T>(19.14904664790693019642068229478769661515L),
         static_cast<T>(-3.542488556926667589704590409095331790317L),
         static_cast<T>(0.3099140334815639910894627700232804503017L),
         static_cast<T>(-0.01036716187296241640634252431913030440825L),
         static_cast<T>(0.8399926504443119927673843789048514017761e-4L),
         static_cast<T>(-0.493038376656195010308610694048822561263e-7L),
      };
      T result = 0;
      T z = dz + 2;
      for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
      {
         result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
      }
      return result;
   }

   static double g(){ return 10.90051099999999983936049829935654997826; }
};

//
// Lanczos Coefficients for N=13 G=13.144565
// Max experimental error (with arbitary precision arithmetic) 9.2213e-23
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos13 : public mpl::int_<72>
{
   //
   // Produces slightly better than extended-double precision when evaluated at
   // higher precision:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[13] = {
         static_cast<T>(44012138428004.60895436261759919070125699L),
         static_cast<T>(41590453358593.20051581730723108131357995L),
         static_cast<T>(18013842787117.99677796276038389462742949L),
         static_cast<T>(4728736263475.388896889723995205703970787L),
         static_cast<T>(837910083628.4046470415724300225777912264L),
         static_cast<T>(105583707273.4299344907359855510105321192L),
         static_cast<T>(9701363618.494999493386608345339104922694L),
         static_cast<T>(654914397.5482052641016767125048538245644L),
         static_cast<T>(32238322.94213356530668889463945849409184L),
         static_cast<T>(1128514.219497091438040721811544858643121L),
         static_cast<T>(26665.79378459858944762533958798805525125L),
         static_cast<T>(381.8801248632926870394389468349331394196L),
         static_cast<T>(2.506628274631000502415763426076722427007L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint32_t) denom[13] = {
         static_cast<boost::uint32_t>(0u),
         static_cast<boost::uint32_t>(39916800u),
         static_cast<boost::uint32_t>(120543840u),
         static_cast<boost::uint32_t>(150917976u),
         static_cast<boost::uint32_t>(105258076u),
         static_cast<boost::uint32_t>(45995730u),
         static_cast<boost::uint32_t>(13339535u),
         static_cast<boost::uint32_t>(2637558u),
         static_cast<boost::uint32_t>(357423u),
         static_cast<boost::uint32_t>(32670u),
         static_cast<boost::uint32_t>(1925u),
         static_cast<boost::uint32_t>(66u),
         static_cast<boost::uint32_t>(1u)
      };
      return boost::math::tools::evaluate_rational(num, denom, z);
   }

   template <class T>
   static T lanczos_sum_expG_scaled(const T& z)
   {
      static const T num[13] = {
         static_cast<T>(86091529.53418537217994842267760536134841L),
         static_cast<T>(81354505.17858011242874285785316135398567L),
         static_cast<T>(35236626.38815461910817650960734605416521L),
         static_cast<T>(9249814.988024471294683815872977672237195L),