lanczos.hpp

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

         static_cast<T>(1639024.216687146960253839656643518985826L),
         static_cast<T>(206530.8157641225032631778026076868855623L),
         static_cast<T>(18976.70193530288915698282139308582105936L),
         static_cast<T>(1281.068909912559479885759622791374106059L),
         static_cast<T>(63.06093343420234536146194868906771599354L),
         static_cast<T>(2.207470909792527638222674678171050209691L),
         static_cast<T>(0.05216058694613505427476207805814960742102L),
         static_cast<T>(0.0007469903808915448316510079585999893674101L),
         static_cast<T>(0.4903180573459871862552197089738373164184e-5L)
      };
      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>(4.832115561461656947793029596285626840312L),
         static_cast<T>(-19.86441536140337740383120735104359034688L),
         static_cast<T>(33.9927422807443239927197864963170585331L),
         static_cast<T>(-31.41520692249765980987427413991250886138L),
         static_cast<T>(17.0270866009599345679868972409543597821L),
         static_cast<T>(-5.5077216950865501362506920516723682167L),
         static_cast<T>(1.037811741948214855286817963800439373362L),
         static_cast<T>(-0.106640468537356182313660880481398642811L),
         static_cast<T>(0.005276450526660653288757565778182586742831L),
         static_cast<T>(-0.0001000935625597121545867453746252064770029L),
         static_cast<T>(0.462590910138598083940803704521211569234e-6L),
         static_cast<T>(-0.1735307814426389420248044907765671743012e-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>(26.96979819614830698367887026728396466395L),
         static_cast<T>(-110.8705424709385114023884328797900204863L),
         static_cast<T>(189.7258846119231466417015694690434770085L),
         static_cast<T>(-175.3397202971107486383321670769397356553L),
         static_cast<T>(95.03437648691551457087250340903980824948L),
         static_cast<T>(-30.7406022781665264273675797983497141978L),
         static_cast<T>(5.792405601630517993355102578874590410552L),
         static_cast<T>(-0.5951993240669148697377539518639997795831L),
         static_cast<T>(0.02944979359164017509944724739946255067671L),
         static_cast<T>(-0.0005586586555377030921194246330399163602684L),
         static_cast<T>(0.2581888478270733025288922038673392636029e-5L),
         static_cast<T>(-0.9685385411006641478305219367315965391289e-9L),
      };
      T result = 0;
      T z = z = 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 13.1445650000000000545696821063756942749; }
};

//
// Lanczos Coefficients for N=22 G=22.61891
// Max experimental error (with arbitary precision arithmetic) 2.9524e-38
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos22 : public mpl::int_<120>
{
   //
   // Produces slightly better than 128-bit long-double precision when 
   // evaluated at higher precision:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[22] = {
         static_cast<T>(46198410803245094237463011094.12173081986L),
         static_cast<T>(43735859291852324413622037436.321513777L),
         static_cast<T>(19716607234435171720534556386.97481377748L),
         static_cast<T>(5629401471315018442177955161.245623932129L),
         static_cast<T>(1142024910634417138386281569.245580222392L),
         static_cast<T>(175048529315951173131586747.695329230778L),
         static_cast<T>(21044290245653709191654675.41581372963167L),
         static_cast<T>(2033001410561031998451380.335553678782601L),
         static_cast<T>(160394318862140953773928.8736211601848891L),
         static_cast<T>(10444944438396359705707.48957290388740896L),
         static_cast<T>(565075825801617290121.1466393747967538948L),
         static_cast<T>(25475874292116227538.99448534450411942597L),
         static_cast<T>(957135055846602154.6720835535232270205725L),
         static_cast<T>(29874506304047462.23662392445173880821515L),
         static_cast<T>(769651310384737.2749087590725764959689181L),
         static_cast<T>(16193289100889.15989633624378404096011797L),
         static_cast<T>(273781151680.6807433264462376754578933261L),
         static_cast<T>(3630485900.32917021712188739762161583295L),
         static_cast<T>(36374352.05577334277856865691538582936484L),
         static_cast<T>(258945.7742115532455441786924971194951043L),
         static_cast<T>(1167.501919472435718934219997431551246996L),
         static_cast<T>(2.50662827463100050241576528481104525333L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint64_t) denom[22] = {
         (0uLL),
         (2432902008176640000uLL),
         (8752948036761600000uLL),
         (13803759753640704000uLL),
         (12870931245150988800uLL),
         (8037811822645051776uLL),
         (3599979517947607200uLL),
         (1206647803780373360uLL),
         (311333643161390640uLL),
         (63030812099294896uLL),
         (10142299865511450uLL),
         (1307535010540395uLL),
         (135585182899530uLL),
         (11310276995381uLL),
         (756111184500uLL),
         (40171771630uLL),
         (1672280820uLL),
         (53327946uLL),
         (1256850uLL),
         (20615uLL),
         (210uLL),
         (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[22] = {
         static_cast<T>(6939996264376682180.277485395074954356211L),
         static_cast<T>(6570067992110214451.87201438870245659384L),
         static_cast<T>(2961859037444440551.986724631496417064121L),
         static_cast<T>(845657339772791245.3541226499766163431651L),
         static_cast<T>(171556737035449095.2475716923888737881837L),
         static_cast<T>(26296059072490867.7822441885603400926007L),
         static_cast<T>(3161305619652108.433798300149816829198706L),
         static_cast<T>(305400596026022.4774396904484542582526472L),
         static_cast<T>(24094681058862.55120507202622377623528108L),
         static_cast<T>(1569055604375.919477574824168939428328839L),
         static_cast<T>(84886558909.02047889339710230696942513159L),
         static_cast<T>(3827024985.166751989686050643579753162298L),
         static_cast<T>(143782298.9273215199098728674282885500522L),
         static_cast<T>(4487794.24541641841336786238909171265944L),
         static_cast<T>(115618.2025760830513505888216285273541959L),
         static_cast<T>(2432.580773108508276957461757328744780439L),
         static_cast<T>(41.12782532742893597168530008461874360191L),
         static_cast<T>(0.5453771709477689805460179187388702295792L),
         static_cast<T>(0.005464211062612080347167337964166505282809L),
         static_cast<T>(0.388992321263586767037090706042788910953e-4L),
         static_cast<T>(0.1753839324538447655939518484052327068859e-6L),
         static_cast<T>(0.3765495513732730583386223384116545391759e-9L)
      };
      static const BOOST_MATH_INT_TABLE_TYPE(T, boost::uint64_t) denom[22] = {
         (0uLL),
         (2432902008176640000uLL),
         (8752948036761600000uLL),
         (13803759753640704000uLL),
         (12870931245150988800uLL),
         (8037811822645051776uLL),
         (3599979517947607200uLL),
         (1206647803780373360uLL),
         (311333643161390640uLL),
         (63030812099294896uLL),
         (10142299865511450uLL),
         (1307535010540395uLL),
         (135585182899530uLL),
         (11310276995381uLL),
         (756111184500uLL),
         (40171771630uLL),
         (1672280820uLL),
         (53327946uLL),
         (1256850uLL),
         (20615uLL),
         (210uLL),
         (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[21] = {
         static_cast<T>(8.318998691953337183034781139546384476554L),
         static_cast<T>(-63.15415991415959158214140353299240638675L),
         static_cast<T>(217.3108224383632868591462242669081540163L),
         static_cast<T>(-448.5134281386108366899784093610397354889L),
         static_cast<T>(619.2903759363285456927248474593012711346L),
         static_cast<T>(-604.1630177420625418522025080080444177046L),
         static_cast<T>(428.8166750424646119935047118287362193314L),
         static_cast<T>(-224.6988753721310913866347429589434550302L),
         static_cast<T>(87.32181627555510833499451817622786940961L),
         static_cast<T>(-25.07866854821128965662498003029199058098L),
         static_cast<T>(5.264398125689025351448861011657789005392L),
         static_cast<T>(-0.792518936256495243383586076579921559914L),
         static_cast<T>(0.08317448364744713773350272460937904691566L),
         static_cast<T>(-0.005845345166274053157781068150827567998882L),
         static_cast<T>(0.0002599412126352082483326238522490030412391L),
         static_cast<T>(-0.6748102079670763884917431338234783496303e-5L),
         static_cast<T>(0.908824383434109002762325095643458603605e-7L),
         static_cast<T>(-0.5299325929309389890892469299969669579725e-9L),
         static_cast<T>(0.994306085859549890267983602248532869362e-12L),
         static_cast<T>(-0.3499893692975262747371544905820891835298e-15L),
         static_cast<T>(0.7260746353663365145454867069182884694961e-20L),
      };
      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[21] = {
         static_cast<T>(75.39272007105208086018421070699575462226L),
         static_cast<T>(-572.3481967049935412452681346759966390319L),
         static_cast<T>(1969.426202741555335078065370698955484358L),
         static_cast<T>(-4064.74968778032030891520063865996757519L),
         static_cast<T>(5612.452614138013929794736248384309574814L),
         static_cast<T>(-5475.357667500026172903620177988213902339L),
         static_cast<T>(3886.243614216111328329547926490398103492L),
         static_cast<T>(-2036.382026072125407192448069428134470564L),
         static_cast<T>(791.3727954936062108045551843636692287652L),
         static_cast<T>(-227.2808432388436552794021219198885223122L),
         static_cast<T>(47.70974355562144229897637024320739257284L),
         static_cast<T>(-7.182373807798293545187073539819697141572L),
         static_cast<T>(0.7537866989631514559601547530490976100468L),
         static_cast<T>(-0.05297470142240154822658739758236594717787L),
         static_cast<T>(0.00235577330936380542539812701472320434133L),
         static_cast<T>(-0.6115613067659273118098229498679502138802e-4L),
         static_cast<T>(0.8236417010170941915758315020695551724181e-6L),
         static_cast<T>(-0.4802628430993048190311242611330072198089e-8L),
         static_cast<T>(0.9011113376981524418952720279739624707342e-11L),
         static_cast<T>(-0.3171854152689711198382455703658589996796e-14L),
         static_cast<T>(0.6580207998808093935798753964580596673177e-19L),
      };
      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 22.61890999999999962710717227309942245483; }
};

//
// Lanczos Coefficients for N=6 G=1.428456135094165802001953125
// Max experimental error (with arbitary precision arithmetic) 8.111667e-8
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos6m24 : public mpl::int_<24>
{
   //
   // Use for float precision, when evaluated as a float:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[6] = {
         static_cast<T>(58.52061591769095910314047740215847630266L),
         static_cast<T>(182.5248962595894264831189414768236280862L),
         static_cast<T>(211.0971093028510041839168287718170827259L),
         static_cast<T>(112.2526547883668146736465390902227161763L),
         static_cast<T>(27.5192015197455403062503721613097825345L),
         static_cast<T>(2.50662858515256974113978724717473206342L)
      };
      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] = {