ntl_rr_lanczos.hpp

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

//  (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_LARGE_LANCZOS_HPP
#define BOOST_LARGE_LANCZOS_HPP

#include <boost/math/bindings/rr.hpp>
#include <boost/math/special_functions/lanczos.hpp>
#include <boost/lexical_cast.hpp>

struct lanczos13UDT
{
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[13] = {
         boost::lexical_cast<T>("44012138428004.60895436261759919070125699"),
         boost::lexical_cast<T>("41590453358593.20051581730723108131357995"),
         boost::lexical_cast<T>("18013842787117.99677796276038389462742949"),
         boost::lexical_cast<T>("4728736263475.388896889723995205703970787"),
         boost::lexical_cast<T>("837910083628.4046470415724300225777912264"),
         boost::lexical_cast<T>("105583707273.4299344907359855510105321192"),
         boost::lexical_cast<T>("9701363618.494999493386608345339104922694"),
         boost::lexical_cast<T>("654914397.5482052641016767125048538245644"),
         boost::lexical_cast<T>("32238322.94213356530668889463945849409184"),
         boost::lexical_cast<T>("1128514.219497091438040721811544858643121"),
         boost::lexical_cast<T>("26665.79378459858944762533958798805525125"),
         boost::lexical_cast<T>("381.8801248632926870394389468349331394196"),
         boost::lexical_cast<T>("2.506628274631000502415763426076722427007"),
      };
      static const 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, 13);
   }

   template <class T>
   static T lanczos_sum_expG_scaled(const T& z)
   {
      static const T num[13] = {
         boost::lexical_cast<T>("86091529.53418537217994842267760536134841"),
         boost::lexical_cast<T>("81354505.17858011242874285785316135398567"),
         boost::lexical_cast<T>("35236626.38815461910817650960734605416521"),
         boost::lexical_cast<T>("9249814.988024471294683815872977672237195"),
         boost::lexical_cast<T>("1639024.216687146960253839656643518985826"),
         boost::lexical_cast<T>("206530.8157641225032631778026076868855623"),
         boost::lexical_cast<T>("18976.70193530288915698282139308582105936"),
         boost::lexical_cast<T>("1281.068909912559479885759622791374106059"),
         boost::lexical_cast<T>("63.06093343420234536146194868906771599354"),
         boost::lexical_cast<T>("2.207470909792527638222674678171050209691"),
         boost::lexical_cast<T>("0.05216058694613505427476207805814960742102"),
         boost::lexical_cast<T>("0.0007469903808915448316510079585999893674101"),
         boost::lexical_cast<T>("0.4903180573459871862552197089738373164184e-5"),
      };
      static const 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, 13);
   }


   template<class T>
   static T lanczos_sum_near_1(const T& dz)
   {
      static const T d[12] = {
         boost::lexical_cast<T>("4.832115561461656947793029596285626840312"),
         boost::lexical_cast<T>("-19.86441536140337740383120735104359034688"),
         boost::lexical_cast<T>("33.9927422807443239927197864963170585331"),
         boost::lexical_cast<T>("-31.41520692249765980987427413991250886138"),
         boost::lexical_cast<T>("17.0270866009599345679868972409543597821"),
         boost::lexical_cast<T>("-5.5077216950865501362506920516723682167"),
         boost::lexical_cast<T>("1.037811741948214855286817963800439373362"),
         boost::lexical_cast<T>("-0.106640468537356182313660880481398642811"),
         boost::lexical_cast<T>("0.005276450526660653288757565778182586742831"),
         boost::lexical_cast<T>("-0.0001000935625597121545867453746252064770029"),
         boost::lexical_cast<T>("0.462590910138598083940803704521211569234e-6"),
         boost::lexical_cast<T>("-0.1735307814426389420248044907765671743012e-9"),
      };
      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] = {
         boost::lexical_cast<T>("26.96979819614830698367887026728396466395"),
         boost::lexical_cast<T>("-110.8705424709385114023884328797900204863"),
         boost::lexical_cast<T>("189.7258846119231466417015694690434770085"),
         boost::lexical_cast<T>("-175.3397202971107486383321670769397356553"),
         boost::lexical_cast<T>("95.03437648691551457087250340903980824948"),
         boost::lexical_cast<T>("-30.7406022781665264273675797983497141978"),
         boost::lexical_cast<T>("5.792405601630517993355102578874590410552"),
         boost::lexical_cast<T>("-0.5951993240669148697377539518639997795831"),
         boost::lexical_cast<T>("0.02944979359164017509944724739946255067671"),
         boost::lexical_cast<T>("-0.0005586586555377030921194246330399163602684"),
         boost::lexical_cast<T>("0.2581888478270733025288922038673392636029e-5"),
         boost::lexical_cast<T>("-0.9685385411006641478305219367315965391289e-9"),
      };
      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 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 lanczos22UDT
{
   template <class T>
   static T lanczos_sum(const T& z)
   {
      static const T num[22] = {
         boost::lexical_cast<T>("46198410803245094237463011094.12173081986"),
         boost::lexical_cast<T>("43735859291852324413622037436.321513777"),
         boost::lexical_cast<T>("19716607234435171720534556386.97481377748"),
         boost::lexical_cast<T>("5629401471315018442177955161.245623932129"),
         boost::lexical_cast<T>("1142024910634417138386281569.245580222392"),
         boost::lexical_cast<T>("175048529315951173131586747.695329230778"),
         boost::lexical_cast<T>("21044290245653709191654675.41581372963167"),
         boost::lexical_cast<T>("2033001410561031998451380.335553678782601"),
         boost::lexical_cast<T>("160394318862140953773928.8736211601848891"),
         boost::lexical_cast<T>("10444944438396359705707.48957290388740896"),
         boost::lexical_cast<T>("565075825801617290121.1466393747967538948"),
         boost::lexical_cast<T>("25475874292116227538.99448534450411942597"),
         boost::lexical_cast<T>("957135055846602154.6720835535232270205725"),
         boost::lexical_cast<T>("29874506304047462.23662392445173880821515"),
         boost::lexical_cast<T>("769651310384737.2749087590725764959689181"),
         boost::lexical_cast<T>("16193289100889.15989633624378404096011797"),
         boost::lexical_cast<T>("273781151680.6807433264462376754578933261"),
         boost::lexical_cast<T>("3630485900.32917021712188739762161583295"),
         boost::lexical_cast<T>("36374352.05577334277856865691538582936484"),
         boost::lexical_cast<T>("258945.7742115532455441786924971194951043"),
         boost::lexical_cast<T>("1167.501919472435718934219997431551246996"),
         boost::lexical_cast<T>("2.50662827463100050241576528481104525333"),
      };
      static const boost::uint64_t denom[22] = {
         (0uLL),
         (2432902008176640000uLL),
         (8752948036761600000uLL),
         (13803759753640704000uLL),