lanczos.hpp

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

         static_cast<T>(0.1346609158752142460943888149156716841693e-14L),
         static_cast<T>(-0.660492688923978805315914918995410340796e-15L),
         static_cast<T>(0.1472114697343266749193617793755763792681e-15L),
         static_cast<T>(-0.1410901942033374651613542904678399264447e-16L),
      };
      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[16] = {
         static_cast<T>(23.56409085052261327114594781581930373708L),
         static_cast<T>(-88.92116338946308797946237246006238652361L),
         static_cast<T>(137.3472822086847596961177383569603988797L),
         static_cast<T>(-112.0400438263562152489272966461114852861L),
         static_cast<T>(51.98768915202973863076166956576777843805L),
         static_cast<T>(-13.78552090862799358221343319574970124948L),
         static_cast<T>(1.996371068830872830250406773917646121742L),
         static_cast<T>(-0.1423525874909934506274738563671862576161L),
         static_cast<T>(0.004098338646046865122459664947239111298524L),
         static_cast<T>(-0.3203286637326511000882086573060433529094e-4L),
         static_cast<T>(0.2641536751640138646146395939004587594407e-7L),
         static_cast<T>(-0.7777876663062235617693516558976641009819e-13L),
         static_cast<T>(0.7061443477097101636871806229515157914789e-14L),
         static_cast<T>(-0.3463537849537988455590834887691613484813e-14L),
         static_cast<T>(0.7719578215795234036320348283011129450595e-15L),
         static_cast<T>(-0.7398586479708476329563577384044188912075e-16L),
      };
      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 12.2252227365970611572265625; }
};

//
// Lanczos Coefficients for N=24 G=20.3209821879863739013671875
// Max experimental error (with arbitary precision arithmetic) 1.0541e-38
// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
//
struct lanczos24m113 : public mpl::int_<113>
{
   //
   // Use for long-double precision, when evaluated as an long-double:
   //
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[24] = {
         static_cast<T>(2029889364934367661624137213253.22102954656825019111612712252027267955023987678816620961507L),
         static_cast<T>(2338599599286656537526273232565.2727349714338768161421882478417543004440597874814359063158L),
         static_cast<T>(1288527989493833400335117708406.3953711906175960449186720680201425446299360322830739180195L),
         static_cast<T>(451779745834728745064649902914.550539158066332484594436145043388809847364393288132164411521L),
         static_cast<T>(113141284461097964029239556815.291212318665536114012605167994061291631013303788706545334708L),
         static_cast<T>(21533689802794625866812941616.7509064680880468667055339259146063256555368135236149614592432L),
         static_cast<T>(3235510315314840089932120340.71494940111731241353655381919722177496659303550321056514776757L),
         static_cast<T>(393537392344185475704891959.081297108513472083749083165179784098220158201055270548272414314L),
         static_cast<T>(39418265082950435024868801.5005452240816902251477336582325944930252142622315101857742955673L),
         static_cast<T>(3290158764187118871697791.05850632319194734270969161036889516414516566453884272345518372696L),
         static_cast<T>(230677110449632078321772.618245845856640677845629174549731890660612368500786684333975350954L),
         static_cast<T>(13652233645509183190158.5916189185218250859402806777406323001463296297553612462737044693697L),
         static_cast<T>(683661466754325350495.216655026531202476397782296585200982429378069417193575896602446904762L),
         static_cast<T>(28967871782219334117.0122379171041074970463982134039409352925258212207710168851968215545064L),
         static_cast<T>(1036104088560167006.2022834098572346459442601718514554488352117620272232373622553429728555L),
         static_cast<T>(31128490785613152.8380102669349814751268126141105475287632676569913936040772990253369753962L),
         static_cast<T>(779327504127342.536207878988196814811198475410572992436243686674896894543126229424358472541L),
         static_cast<T>(16067543181294.643350688789124777020407337133926174150582333950666044399234540521336771876L),
         static_cast<T>(268161795520.300916569439413185778557212729611517883948634711190170998896514639936969855484L),
         static_cast<T>(3533216359.10528191668842486732408440112703691790824611391987708562111396961696753452085068L),
         static_cast<T>(35378979.5479656110614685178752543826919239614088343789329169535932709470588426584501652577L),
         static_cast<T>(253034.881362204346444503097491737872930637147096453940375713745904094735506180552724766444L),
         static_cast<T>(1151.61895453463992438325318456328526085882924197763140514450975619271382783957699017875304L),
         static_cast<T>(2.50662827463100050241576528481104515966515623051532908941425544355490413900497467936202516L)
      };
      static const T denom[24] = {
         static_cast<T>(0L),
         static_cast<T>(0.112400072777760768e22L),
         static_cast<T>(0.414847677933545472e22L),
         static_cast<T>(6756146673770930688000.0L),
         static_cast<T>(6548684852703068697600.0L),
         static_cast<T>(4280722865357147142912.0L),
         static_cast<T>(2021687376910682741568.0L),
         static_cast<T>(720308216440924653696.0L),
         static_cast<T>(199321978221066137360.0L),
         static_cast<T>(43714229649594412832.0L),
         static_cast<T>(7707401101297361068.0L),
         static_cast<T>(1103230881185949736.0L),
         static_cast<T>(129006659818331295.0L),
         static_cast<T>(12363045847086207.0L),
         static_cast<T>(971250460939913.0L),
         static_cast<T>(62382416421941.0L),
         static_cast<T>(3256091103430.0L),
         static_cast<T>(136717357942.0L),
         static_cast<T>(4546047198.0L),
         static_cast<T>(116896626L),
         static_cast<T>(2240315L),
         static_cast<T>(30107L),
         static_cast<T>(253L),
         static_cast<T>(1L)
      };
      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[24] = {
         static_cast<T>(3035162425359883494754.02878223286972654682199012688209026810841953293372712802258398358538L),
         static_cast<T>(3496756894406430103600.16057175075063458536101374170860226963245118484234495645518505519827L),
         static_cast<T>(1926652656689320888654.01954015145958293168365236755537645929361841917596501251362171653478L),
         static_cast<T>(675517066488272766316.083023742440619929434602223726894748181327187670231286180156444871912L),
         static_cast<T>(169172853104918752780.086262749564831660238912144573032141700464995906149421555926000038492L),
         static_cast<T>(32197935167225605785.6444116302160245528783954573163541751756353183343357329404208062043808L),
         static_cast<T>(4837849542714083249.37587447454818124327561966323276633775195138872820542242539845253171632L),
         static_cast<T>(588431038090493242.308438203986649553459461798968819276505178004064031201740043314534404158L),
         static_cast<T>(58939585141634058.6206417889192563007809470547755357240808035714047014324843817783741669733L),
         static_cast<T>(4919561837722192.82991866530802080996138070630296720420704876654726991998309206256077395868L),
         static_cast<T>(344916580244240.407442753122831512004021081677987651622305356145640394384006997569631719101L),
         static_cast<T>(20413302960687.8250598845969238472629322716685686993835561234733641729957841485003560103066L),
         static_cast<T>(1022234822943.78400752460970689311934727763870970686747383486600540378889311406851534545789L),
         static_cast<T>(43313787191.9821354846952908076307094286897439975815501673706144217246093900159173598852503L),
         static_cast<T>(1549219505.59667418528481770869280437577581951167003505825834192510436144666564648361001914L),
         static_cast<T>(46544421.1998761919380541579358096705925369145324466147390364674998568485110045455014967149L),
         static_cast<T>(1165278.06807504975090675074910052763026564833951579556132777702952882101173607903881127542L),
         static_cast<T>(24024.759267256769471083727721827405338569868270177779485912486668586611981795179894572115L),
         static_cast<T>(400.965008113421955824358063769761286758463521789765880962939528760888853281920872064838918L),
         static_cast<T>(5.28299015654478269617039029170846385138134929147421558771949982217659507918482272439717603L),
         static_cast<T>(0.0528999024412510102409256676599360516359062802002483877724963720047531347449011629466149805L),
         static_cast<T>(0.000378346710654740685454266569593414561162134092347356968516522170279688139165340746957511115L),
         static_cast<T>(0.172194142179211139195966608011235161516824700287310869949928393345257114743230967204370963e-5L),
         static_cast<T>(0.374799931707148855771381263542708435935402853962736029347951399323367765509988401336565436e-8L)
      };
      static const T denom[24] = {
         static_cast<T>(0L),
         static_cast<T>(0.112400072777760768e22L),
         static_cast<T>(0.414847677933545472e22L),
         static_cast<T>(6756146673770930688000.0L),
         static_cast<T>(6548684852703068697600.0L),
         static_cast<T>(4280722865357147142912.0L),
         static_cast<T>(2021687376910682741568.0L),
         static_cast<T>(720308216440924653696.0L),
         static_cast<T>(199321978221066137360.0L),
         static_cast<T>(43714229649594412832.0L),
         static_cast<T>(7707401101297361068.0L),
         static_cast<T>(1103230881185949736.0L),
         static_cast<T>(129006659818331295.0L),
         static_cast<T>(12363045847086207.0L),
         static_cast<T>(971250460939913.0L),
         static_cast<T>(62382416421941.0L),
         static_cast<T>(3256091103430.0L),
         static_cast<T>(136717357942.0L),
         static_cast<T>(4546047198.0L),
         static_cast<T>(116896626L),
         static_cast<T>(2240315L),
         static_cast<T>(30107L),
         static_cast<T>(253L),
         static_cast<T>(1L)
      };
      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[23] = {
         static_cast<T>(7.4734083002469026177867421609938203388868806387315406134072298925733950040583068760685908L),
         static_cast<T>(-50.4225805042247530267317342133388132970816607563062253708655085754357843064134941138154171L),
         static_cast<T>(152.288200621747008570784082624444625293884063492396162110698238568311211546361189979357019L),
         static_cast<T>(-271.894959539150384169327513139846971255640842175739337449692360299099322742181325023644769L),
         static_cast<T>(319.240102980202312307047586791116902719088581839891008532114107693294261542869734803906793L),
         static_cast<T>(-259.493144143048088289689500935518073716201741349569864988870534417890269467336454358361499L),
         static_cast<T>(149.747518319689708813209645403067832020714660918583227716408482877303972685262557460145835L),
         static_cast<T>(-61.9261301009341333289187201425188698128684426428003249782448828881580630606817104372760037L),
         static_cast<T>(18.3077524177286961563937379403377462608113523887554047531153187277072451294845795496072365L),
         static_cast<T>(-3.82011322251948043097070160584761236869363471824695092089556195047949392738162970152230254L),
         static_cast<T>(0.549382685505691522516705902336780999493262538301283190963770663549981309645795228539620711L),
         static_cast<T>(-0.0524814679715180697633723771076668718265358076235229045603747927518423453658004287459638024L),
         static_cast<T>(0.00315392664003333528534120626687784812050217700942910879712808180705014754163256855643360698L),
         static_cast<T>(-0.000110098373127648510519799564665442121339511198561008748083409549601095293123407080388658329L),
         static_cast<T>(0.19809382866681658224945717689377373458866950897791116315219376038432014207446832310901893e-5L),
         static_cast<T>(-0.152278977408600291408265615203504153130482270424202400677280558181047344681214058227949755e-7L),
         static_cast<T>(0.364344768076106268872239259083188037615571711218395765792787047015406264051536972018235217e-10L),
         static_cast<T>(-0.148897510480440424971521542520683536298361220674662555578951242811522959610991621951203526e-13L),
         static_cast<T>(0.261199241161582662426512749820666625442516059622425213340053324061794752786482115387573582e-18L),
         static_cast<T>(-0.780072664167099103420998436901014795601783313858454665485256897090476089641613851903791529e-24L),
         static_cast<T>(0.303465867587106629530056603454807425512962762653755513440561256044986695349304176849392735e-24L),
         static_cast<T>(-0.615420597971283870342083342286977366161772327800327789325710571275345878439656918541092056e-25L),
         static_cast<T>(0.499641233843540749369110053005439398774706583601830828776209650445427083113181961630763702e-26L),
      };
      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[23] = {
         static_cast<T>(61.4165001061101455341808888883960361969557848005400286332291451422461117307237198559485365L),
         static_cast<T>(-414.372973678657049667308134761613915623353625332248315105320470271523320700386200587519147L),
         static_cast<T>(1251.50505818554680171298972755376376836161706773644771875668053742215217922228357204561873L),
         static_cast<T>(-2234.43389421602399514176336175766511311493214354568097811220122848998413358085613880612158L),
         static_cast<T>(2623.51647746991904821899989145639147785427273427135380151752779100215839537090464785708684L),
         static_cast<T>(-2132.51572435428751962745870184529534443305617818870214348386131243463614597272260797772423L),
         static_cast<T>(1230.62572059218405766499842067263311220019173335523810725664442147670956427061920234820189L),
         static_cast<T>(-508.90919151163744999377586956023909888833335885805154492270846381061182696305011395981929L),
         static_cast<T>(150.453184562246579758706538566480316921938628645961177699894388251635886834047343195475395L),
         static_cast<T>(-31.3937061525822497422230490071156186113405446381476081565548185848237169870395131828731397L),
         static_cast<T>(4.51482916590287954234936829724231512565732528859217337795452389161322923867318809206313688L),
         static_cast<T>(-0.431292919341108177524462194102701868233551186625103849565527515201492276412231365776131952L),
         static_cast<T>(0.0259189820815586225636729971503340447445001375909094681698918294680345547092233915092128323L),
         static_cast<T>(-0.000904788882557558697594884691337532557729219389814315972435534723829065673966567231504429712L),
         static_cast<T>(0.162793589759218213439218473348810982422449144393340433592232065020562974405674317564164312e-4L),
         static_cast<T>(-0.125142926178202562426432039899709511761368233479483128438847484617555752948755923647214487e-6L),
         static_cast<T>(0.299418680048132583204152682950097239197934281178261879500770485862852229898797687301941982e-9L),
         static_cast<T>(-0.122364035267809278675627784883078206654408225276233049012165202996967011873995261617995421e-12L),
         static_cast<T>(0.21465364366598631597052073538883430194257709353929022544344097235100199405814005393447785e-17L),
         static_cast<T>(-0.641064035802907518396608051803921688237330857546406669209280666066685733941549058513986818e-23L),
         static_cast<T>(0.249388374622173329690271566855185869111237201309011956145463506483151054813346819490278951e-23L),
         static_cast<T>(-0.505752900177513489906064295001851463338022055787536494321532352380960774349054239257683149e-24L),
         static_cast<T>(0.410605371184590959139968810080063542546949719163227555918846829816144878123034347778284006e-25L),
      };
      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 20.3209821879863739013671875; }
};


//
// placeholder for no lanczos info available:
//
struct undefined_lanczos : public mpl::int_<INT_MAX - 1> { };

#if 0
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
#define BOOST_MATH_FLT_DIGITS ::std::numeric_limits<float>::digits
#define BOOST_MATH_DBL_DIGITS ::std::numeric_limits<double>::digits
#define BOOST_MATH_LDBL_DIGITS ::std::numeric_limits<long double>::digits
#else
#define BOOST_MATH_FLT_DIGITS FLT_MANT_DIG
#define BOOST_MATH_DBL_DIGITS DBL_MANT_DIG
#define BOOST_MATH_LDBL_DIGITS LDBL_MANT_DIG
#endif
#endif

typedef mpl::list<
   lanczos6m24, 
/*   lanczos6, */
   lanczos13m53, 
/*   lanczos13, */
   lanczos17m64, 
   lanczos24m113, 
   lanczos22, 
   undefined_lanczos> lanczos_list;

template <class Real, class Policy>
struct lanczos
{
   typedef typename mpl::if_<
      typename mpl::less_equal<
         typename policies::precision<Real, Policy>::type,
         mpl::int_<0>
      >::type,
      mpl::int_<INT_MAX - 2>,
      typename policies::precision<Real, Policy>::type
   >::type target_precision;

   typedef typename mpl::deref<typename mpl::find_if<
      lanczos_list, 
      mpl::less_equal<target_precision, mpl::_1> >::type>::type type;
};

} // namespace lanczos
} // namespace math
} // namespace boost

#if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__)
#include <boost/math/special_functions/detail/lanczos_sse2.hpp>
#endif

#endif // BOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS