zeta.hpp

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

   else if(s < 56)
   {
      result = 1 + pow(T(2), -s);
   }
   else
   {
      result = 1;
   }
   return result;
}

template <class T, class Policy>
T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<64>&)
{
   BOOST_MATH_STD_USING
   T result;
   if(s < 1)
   {
      // Rational Approximation
      // Maximum Deviation Found:                     3.099e-20
      // Expected Error Term:                         3.099e-20
      // Max error found at long double precision:    5.890498e-20
      static const T P[6] = {    
         0.243392944335937499969L,
         -0.496837806864865688082L,
         0.0680008039723709987107L,
         -0.00511620413006619942112L,
         0.000455369899250053003335L,
         -0.279496685273033761927e-4L,
        };
      static const T Q[7] = {    
         1L,
         -0.30425480068225790522L,
         0.050052748580371598736L,
         -0.00519355671064700627862L,
         0.000360623385771198350257L,
         -0.159600883054550987633e-4L,
         0.339770279812410586032e-6L,
      };
      result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
      result -= 1.2433929443359375F;
      result += (sc);
      result /= (sc);
   }
   else if(s <= 2)
   {
      // Maximum Deviation Found:                     1.059e-21
      // Expected Error Term:                         1.059e-21
      // Max error found at long double precision:    1.626303e-19

      static const T P[6] = {    
         0.577215664901532860605L,
         0.222537368917162139445L,
         0.0356286324033215682729L,
         0.00304465292366350081446L,
         0.000178102511649069421904L,
         0.700867470265983665042e-5L,
      };
      static const T Q[7] = {    
         1L,
         0.259385759149531030085L,
         0.0373974962106091316854L,
         0.00332735159183332820617L,
         0.000188690420706998606469L,
         0.635994377921861930071e-5L,
         0.226583954978371199405e-7L,
      };
      result = tools::evaluate_polynomial(P, -sc) / tools::evaluate_polynomial(Q, -sc);
      result += 1 / (-sc);
   }
   else if(s <= 4)
   {
      // Maximum Deviation Found:          5.946e-22
      // Expected Error Term:              -5.946e-22
      static const float Y = 0.6986598968505859375;
      static const T P[7] = {    
         -0.053725830002359501027L,
         0.0470551187571475844778L,
         0.0101339410415759517471L,
         0.00100240326666092854528L,
         0.685027119098122814867e-4L,
         0.390972820219765942117e-5L,
         0.540319769113543934483e-7L,
      };
      static const T Q[8] = {    
         1,
         0.286577739726542730421L,
         0.0447355811517733225843L,
         0.00430125107610252363302L,
         0.000284956969089786662045L,
         0.116188101609848411329e-4L,
         0.278090318191657278204e-6L,
         -0.19683620233222028478e-8L,
      };
      result = tools::evaluate_polynomial(P, s - 2) / tools::evaluate_polynomial(Q, s - 2);
      result += Y + 1 / (-sc);
   }
   else if(s <= 7)
   {
      // Max error found at long double precision: 8.132216e-19
      static const T P[8] = {    
         -2.49710190602259407065L,
         -3.36664913245960625334L,
         -1.77180020623777595452L,
         -0.464717885249654313933L,
         -0.0643694921293579472583L,
         -0.00464265386202805715487L,
         -0.000165556579779704340166L,
         -0.252884970740994069582e-5L,
      };
      static const T Q[9] = {    
         1,
         1.01300131390690459085L,
         0.387898115758643503827L,
         0.0695071490045701135188L,
         0.00586908595251442839291L,
         0.000217752974064612188616L,
         0.397626583349419011731e-5L,
         -0.927884739284359700764e-8L,
         0.119810501805618894381e-9L,
      };
      result = tools::evaluate_polynomial(P, s - 4) / tools::evaluate_polynomial(Q, s - 4);
      result = 1 + exp(result);
   }
   else if(s < 15)
   {
      // Max error in interpolated form:              1.133e-18
      // Max error found at long double precision:    2.183198e-18
      static const T P[9] = {    
         -4.78558028495135548083L,
         -3.23873322238609358947L,
         -0.892338582881021799922L,
         -0.131326296217965913809L,
         -0.0115651591773783712996L,
         -0.000657728968362695775205L,
         -0.252051328129449973047e-4L,
         -0.626503445372641798925e-6L,
         -0.815696314790853893484e-8L,
        };
      static const T Q[9] = {    
         1,
         0.525765665400123515036L,
         0.10852641753657122787L,
         0.0115669945375362045249L,
         0.000732896513858274091966L,
         0.30683952282420248448e-4L,
         0.819649214609633126119e-6L,
         0.117957556472335968146e-7L,
         -0.193432300973017671137e-12L,
        };
      result = tools::evaluate_polynomial(P, s - 7) / tools::evaluate_polynomial(Q, s - 7);
      result = 1 + exp(result);
   }
   else if(s < 42)
   {
      // Max error in interpolated form:             1.668e-17
      // Max error found at long double precision:   1.669714e-17
      static const T P[9] = {    
         -10.3948950573308861781L,
         -2.82646012777913950108L,
         -0.342144362739570333665L,
         -0.0249285145498722647472L,
         -0.00122493108848097114118L,
         -0.423055371192592850196e-4L,
         -0.1025215577185967488e-5L,
         -0.165096762663509467061e-7L,
         -0.145392555873022044329e-9L,
      };
      static const T Q[10] = {    
         1,
         0.205135978585281988052L,
         0.0192359357875879453602L,
         0.00111496452029715514119L,
         0.434928449016693986857e-4L,
         0.116911068726610725891e-5L,
         0.206704342290235237475e-7L,
         0.209772836100827647474e-9L,
         -0.939798249922234703384e-16L,
         0.264584017421245080294e-18L,
      };
      result = tools::evaluate_polynomial(P, s - 15) / tools::evaluate_polynomial(Q, s - 15);
      result = 1 + exp(result);
   }
   else if(s < 63)
   {
      result = 1 + pow(T(2), -s);
   }
   else
   {
      result = 1;
   }
   return result;
}

template <class T, class Policy>
T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<113>&)
{
   BOOST_MATH_STD_USING
   T result;
   if(s < 1)
   {
      // Rational Approximation
      // Maximum Deviation Found:                     9.493e-37
      // Expected Error Term:                         9.492e-37
      // Max error found at long double precision:    7.281332e-31

      static const T P[10] = {    
         -1L,
         -0.0353008629988648122808504280990313668L,
         0.0107795651204927743049369868548706909L,
         0.000523961870530500751114866884685172975L,
         -0.661805838304910731947595897966487515e-4L,
         -0.658932670403818558510656304189164638e-5L,
         -0.103437265642266106533814021041010453e-6L,
         0.116818787212666457105375746642927737e-7L,
         0.660690993901506912123512551294239036e-9L,
         0.113103113698388531428914333768142527e-10L,
        };
      static const T Q[11] = {    
         1L,
         -0.387483472099602327112637481818565459L,
         0.0802265315091063135271497708694776875L,
         -0.0110727276164171919280036408995078164L,
         0.00112552716946286252000434849173787243L,
         -0.874554160748626916455655180296834352e-4L,
         0.530097847491828379568636739662278322e-5L,
         -0.248461553590496154705565904497247452e-6L,
         0.881834921354014787309644951507523899e-8L,
         -0.217062446168217797598596496310953025e-9L,
         0.315823200002384492377987848307151168e-11L,
      };
      result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
      result += (sc);
      result /= (sc);
   }
   else if(s <= 2)
   {
      // Maximum Deviation Found:                     1.616e-37
      // Expected Error Term:                         -1.615e-37

      static const T P[10] = {    
         0.577215664901532860606512090082402431L,
         0.255597968739771510415479842335906308L,
         0.0494056503552807274142218876983542205L,
         0.00551372778611700965268920983472292325L,
         0.00043667616723970574871427830895192731L,
         0.268562259154821957743669387915239528e-4L,
         0.109249633923016310141743084480436612e-5L,
         0.273895554345300227466534378753023924e-7L,
         0.583103205551702720149237384027795038e-9L,
         -0.835774625259919268768735944711219256e-11L,
      };
      static const T Q[11] = {    
         1L,
         0.316661751179735502065583176348292881L,
         0.0540401806533507064453851182728635272L,
         0.00598621274107420237785899476374043797L,
         0.000474907812321704156213038740142079615L,
         0.272125421722314389581695715835862418e-4L,
         0.112649552156479800925522445229212933e-5L,
         0.301838975502992622733000078063330461e-7L,
         0.422960728687211282539769943184270106e-9L,
         -0.377105263588822468076813329270698909e-11L,
         -0.581926559304525152432462127383600681e-13L,
      };
      result = tools::evaluate_polynomial(P, -sc) / tools::evaluate_polynomial(Q, -sc);
      result += 1 / (-sc);
   }
   else if(s <= 4)
   {
      // Maximum Deviation Found:                     1.891e-36
      // Expected Error Term:                         -1.891e-36
      // Max error found: 2.171527e-35

      static const float Y = 0.6986598968505859375;
      static const T P[11] = {    
         -0.0537258300023595010275848333539748089L,
         0.0429086930802630159457448174466342553L,
         0.0136148228754303412510213395034056857L,
         0.00190231601036042925183751238033763915L,
         0.000186880390916311438818302549192456581L,
         0.145347370745893262394287982691323657e-4L,
         0.805843276446813106414036600485884885e-6L,
         0.340818159286739137503297172091882574e-7L,
         0.115762357488748996526167305116837246e-8L,
         0.231904754577648077579913403645767214e-10L,
         0.340169592866058506675897646629036044e-12L,
      };
      static const T Q[12] = {    
         1L,
         0.363755247765087100018556983050520554L,
         0.0696581979014242539385695131258321598L,
         0.00882208914484611029571547753782014817L,
         0.000815405623261946661762236085660996718L,
         0.571366167062457197282642344940445452e-4L,
         0.309278269271853502353954062051797838e-5L,
         0.12822982083479010834070516053794262e-6L,
         0.397876357325018976733953479182110033e-8L,
         0.8484432107648683277598472295289279e-10L,
         0.105677416606909614301995218444080615e-11L,
         0.547223964564003701979951154093005354e-15L,
      };