lanczos.hpp

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

         static_cast<T>(14.0261432874996476619570577285003839357L),
         static_cast<T>(43.74732405540314316089531289293124360129L),
         static_cast<T>(50.59547402616588964511581430025589038612L),
         static_cast<T>(26.90456680562548195593733429204228910299L),
         static_cast<T>(6.595765571169314946316366571954421695196L),
         static_cast<T>(0.6007854010515290065101128585795542383721L)
      };
      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>(0.4922488055204602807654354732674868442106L),
         static_cast<T>(0.004954497451132152436631238060933905650346L),
         static_cast<T>(-0.003374784572167105840686977985330859371848L),
         static_cast<T>(0.001924276018962061937026396537786414831385L),
         static_cast<T>(-0.00056533046336427583708166383712907694434L),
      };
      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>(0.6534966888520080645505805298901130485464L),
         static_cast<T>(0.006577461728560758362509168026049182707101L),
         static_cast<T>(-0.004480276069269967207178373559014835978161L),
         static_cast<T>(0.00255461870648818292376982818026706528842L),
         static_cast<T>(-0.000750517993690428370380996157470900204524L),
      };
      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 1.428456135094165802001953125; }
};

//
// Lanczos Coefficients for N=13 G=6.024680040776729583740234375
// Max experimental error (with arbitary precision arithmetic) 1.196214e-17
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos13m53 : public mpl::int_<53>
{
   //
   // Use for double precision, when evaluated as a double:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[13] = {
         static_cast<T>(23531376880.41075968857200767445163675473L),
         static_cast<T>(42919803642.64909876895789904700198885093L),
         static_cast<T>(35711959237.35566804944018545154716670596L),
         static_cast<T>(17921034426.03720969991975575445893111267L),
         static_cast<T>(6039542586.35202800506429164430729792107L),
         static_cast<T>(1439720407.311721673663223072794912393972L),
         static_cast<T>(248874557.8620541565114603864132294232163L),
         static_cast<T>(31426415.58540019438061423162831820536287L),
         static_cast<T>(2876370.628935372441225409051620849613599L),
         static_cast<T>(186056.2653952234950402949897160456992822L),
         static_cast<T>(8071.672002365816210638002902272250613822L),
         static_cast<T>(210.8242777515793458725097339207133627117L),
         static_cast<T>(2.506628274631000270164908177133837338626L)
      };
      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>(56906521.91347156388090791033559122686859L),
         static_cast<T>(103794043.1163445451906271053616070238554L),
         static_cast<T>(86363131.28813859145546927288977868422342L),
         static_cast<T>(43338889.32467613834773723740590533316085L),
         static_cast<T>(14605578.08768506808414169982791359218571L),
         static_cast<T>(3481712.15498064590882071018964774556468L),
         static_cast<T>(601859.6171681098786670226533699352302507L),
         static_cast<T>(75999.29304014542649875303443598909137092L),
         static_cast<T>(6955.999602515376140356310115515198987526L),
         static_cast<T>(449.9445569063168119446858607650988409623L),
         static_cast<T>(19.51992788247617482847860966235652136208L),
         static_cast<T>(0.5098416655656676188125178644804694509993L),
         static_cast<T>(0.006061842346248906525783753964555936883222L)
      };
      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_near_1(const T& dz)
   {
      static const T d[12] = {
         static_cast<T>(2.208709979316623790862569924861841433016L),
         static_cast<T>(-3.327150580651624233553677113928873034916L),
         static_cast<T>(1.483082862367253753040442933770164111678L),
         static_cast<T>(-0.1993758927614728757314233026257810172008L),
         static_cast<T>(0.004785200610085071473880915854204301886437L),
         static_cast<T>(-0.1515973019871092388943437623825208095123e-5L),
         static_cast<T>(-0.2752907702903126466004207345038327818713e-7L),
         static_cast<T>(0.3075580174791348492737947340039992829546e-7L),
         static_cast<T>(-0.1933117898880828348692541394841204288047e-7L),
         static_cast<T>(0.8690926181038057039526127422002498960172e-8L),
         static_cast<T>(-0.2499505151487868335680273909354071938387e-8L),
         static_cast<T>(0.3394643171893132535170101292240837927725e-9L),
      };
      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[12] = {
         static_cast<T>(6.565936202082889535528455955485877361223L),
         static_cast<T>(-9.8907772644920670589288081640128194231L),
         static_cast<T>(4.408830289125943377923077727900630927902L),
         static_cast<T>(-0.5926941084905061794445733628891024027949L),
         static_cast<T>(0.01422519127192419234315002746252160965831L),
         static_cast<T>(-0.4506604409707170077136555010018549819192e-5L),
         static_cast<T>(-0.8183698410724358930823737982119474130069e-7L),
         static_cast<T>(0.9142922068165324132060550591210267992072e-7L),
         static_cast<T>(-0.5746670642147041587497159649318454348117e-7L),
         static_cast<T>(0.2583592566524439230844378948704262291927e-7L),
         static_cast<T>(-0.7430396708998719707642735577238449585822e-8L),
         static_cast<T>(0.1009141566987569892221439918230042368112e-8L),
      };
      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 6.024680040776729583740234375; }
};

//
// Lanczos Coefficients for N=17 G=12.2252227365970611572265625
// Max experimental error (with arbitary precision arithmetic) 2.7699e-26
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos17m64 : public mpl::int_<64>
{
   //
   // Use for extended-double precision, when evaluated as an extended-double:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[17] = {
         static_cast<T>(553681095419291969.2230556393350368550504L),
         static_cast<T>(731918863887667017.2511276782146694632234L),
         static_cast<T>(453393234285807339.4627124634539085143364L),
         static_cast<T>(174701893724452790.3546219631779712198035L),
         static_cast<T>(46866125995234723.82897281620357050883077L),
         static_cast<T>(9281280675933215.169109622777099699054272L),
         static_cast<T>(1403600894156674.551057997617468721789536L),
         static_cast<T>(165345984157572.7305349809894046783973837L),
         static_cast<T>(15333629842677.31531822808737907246817024L),
         static_cast<T>(1123152927963.956626161137169462874517318L),
         static_cast<T>(64763127437.92329018717775593533620578237L),
         static_cast<T>(2908830362.657527782848828237106640944457L),
         static_cast<T>(99764700.56999856729959383751710026787811L),
         static_cast<T>(2525791.604886139959837791244686290089331L),
         static_cast<T>(44516.94034970167828580039370201346554872L),
         static_cast<T>(488.0063567520005730476791712814838113252L),
         static_cast<T>(2.50662827463100050241576877135758834683L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint64_t) denom[17] = {
         (0uLL),
         (1307674368000uLL),
         (4339163001600uLL),
         (6165817614720uLL),
         (5056995703824uLL),
         (2706813345600uLL),
         (1009672107080uLL),
         (272803210680uLL),
         (54631129553uLL),
         (8207628000uLL),
         (928095740uLL),
         (78558480uLL),
         (4899622uLL),
         (218400uLL),
         (6580uLL),
         (120uLL),
         (1uLL)
      };
      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[17] = {
         static_cast<T>(2715894658327.717377557655133124376674911L),
         static_cast<T>(3590179526097.912105038525528721129550434L),
         static_cast<T>(2223966599737.814969312127353235818710172L),
         static_cast<T>(856940834518.9562481809925866825485883417L),
         static_cast<T>(229885871668.749072933597446453399395469L),
         static_cast<T>(45526171687.54610815813502794395753410032L),
         static_cast<T>(6884887713.165178784550917647709216424823L),
         static_cast<T>(811048596.1407531864760282453852372777439L),
         static_cast<T>(75213915.96540822314499613623119501704812L),
         static_cast<T>(5509245.417224265151697527957954952830126L),
         static_cast<T>(317673.5368435419126714931842182369574221L),
         static_cast<T>(14268.27989845035520147014373320337523596L),
         static_cast<T>(489.3618720403263670213909083601787814792L),
         static_cast<T>(12.38941330038454449295883217865458609584L),
         static_cast<T>(0.2183627389504614963941574507281683147897L),
         static_cast<T>(0.002393749522058449186690627996063983095463L),
         static_cast<T>(0.1229541408909435212800785616808830746135e-4L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint64_t) denom[17] = {
         (0uLL),
         (1307674368000uLL),
         (4339163001600uLL),
         (6165817614720uLL),
         (5056995703824uLL),
         (2706813345600uLL),
         (1009672107080uLL),
         (272803210680uLL),
         (54631129553uLL),
         (8207628000uLL),
         (928095740uLL),
         (78558480uLL),
         (4899622uLL),
         (218400uLL),
         (6580uLL),
         (120uLL),
         (1uLL)
      };
      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[16] = {
         static_cast<T>(4.493645054286536365763334986866616581265L),
         static_cast<T>(-16.95716370392468543800733966378143997694L),
         static_cast<T>(26.19196892983737527836811770970479846644L),
         static_cast<T>(-21.3659076437988814488356323758179283908L),
         static_cast<T>(9.913992596774556590710751047594507535764L),
         static_cast<T>(-2.62888300018780199210536267080940382158L),
         static_cast<T>(0.3807056693542503606384861890663080735588L),
         static_cast<T>(-0.02714647489697685807340312061034730486958L),
         static_cast<T>(0.0007815484715461206757220527133967191796747L),
         static_cast<T>(-0.6108630817371501052576880554048972272435e-5L),
         static_cast<T>(0.5037380238864836824167713635482801545086e-8L),
         static_cast<T>(-0.1483232144262638814568926925964858237006e-13L),