zeta.hpp

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

//  Copyright John Maddock 2007.
//  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_MATH_ZETA_HPP
#define BOOST_MATH_ZETA_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/series.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/sin_pi.hpp>

namespace boost{ namespace math{ namespace detail{

template <class T, class Policy>
struct zeta_series_cache_size
{
   //
   // Work how large to make our cache size when evaluating the series 
   // evaluation:  normally this is just large enough for the series
   // to have converged, but for arbitrary precision types we need a 
   // really large cache to achieve reasonable precision in a reasonable 
   // time.  This is important when constructing rational approximations
   // to zeta for example.
   //
   typedef typename boost::math::policies::precision<T,Policy>::type precision_type;
   typedef typename mpl::if_<
      mpl::less_equal<precision_type, mpl::int_<0> >,
      mpl::int_<5000>,
      typename mpl::if_<
         mpl::less_equal<precision_type, mpl::int_<64> >,
         mpl::int_<70>,
         typename mpl::if_<
            mpl::less_equal<precision_type, mpl::int_<113> >,
            mpl::int_<100>,
            mpl::int_<5000>
         >::type
      >::type
   >::type type;
};

template <class T, class Policy>
T zeta_series_imp(T s, T sc, const Policy&)
{
   //
   // Series evaluation from:
   // Havil, J. Gamma: Exploring Euler's Constant. 
   // Princeton, NJ: Princeton University Press, 2003.
   //
   // See also http://mathworld.wolfram.com/RiemannZetaFunction.html
   //
   BOOST_MATH_STD_USING
   T sum = 0;
   T mult = 0.5;
   T change;
   typedef typename zeta_series_cache_size<T,Policy>::type cache_size;
   T powers[cache_size::value] = { 0, };
   unsigned n = 0;
   do{
      T binom = -static_cast<T>(n);
      T nested_sum = 1;
      if(n < sizeof(powers) / sizeof(powers[0]))
         powers[n] = pow(static_cast<T>(n + 1), -s);
      for(unsigned k = 1; k <= n; ++k)
      {
         T p;
         if(k < sizeof(powers) / sizeof(powers[0]))
         {
            p = powers[k];
            //p = pow(k + 1, -s);
         }
         else
            p = pow(static_cast<T>(k + 1), -s);
         nested_sum += binom * p;
        binom *= (k - static_cast<T>(n)) / (k + 1);
      }
      change = mult * nested_sum;
      sum += change;
      mult /= 2;
      ++n;
   }while(fabs(change / sum) > tools::epsilon<T>());

   return sum * 1 / -boost::math::powm1(T(2), sc);
}
//
// Classical p-series:
//
template <class T>
struct zeta_series2
{
   typedef T result_type;
   zeta_series2(T _s) : s(-_s), k(1){}
   T operator()()
   {
      BOOST_MATH_STD_USING
      return pow(static_cast<T>(k++), s);
   }
private:
   T s;
   unsigned k;
};

template <class T, class Policy>
inline T zeta_series2_imp(T s, const Policy& pol)
{
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();;
   zeta_series2<T> f(s);
   T result = tools::sum_series(
      f, 
      tools::digits<T>(),
      max_iter);
   policies::check_series_iterations("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
   return result;
}

template <class T, class Policy>
T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<0>&)
{
   BOOST_MATH_STD_USING
   T result; 
   //
   // Only use power series if it will converge in 100 
   // iterations or less: the more iterations it consumes
   // the slower convergence becomes so we have to be very 
   // careful in it's usage.
   //
   if (s > -log(tools::epsilon<T>()) / 4.5)
      result = detail::zeta_series2_imp(s, pol);
   else
      result = detail::zeta_series_imp(s, sc, pol);
   return result;
}

template <class T, class Policy>
inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
{
   BOOST_MATH_STD_USING
   T result;
   if(s < 1)
   {
      // Rational Approximation
      // Maximum Deviation Found:                     2.020e-18
      // Expected Error Term:                         -2.020e-18
      // Max error found at double precision:         3.994987e-17
      static const T P[6] = {    
         0.24339294433593750202L,
         -0.49092470516353571651L,
         0.0557616214776046784287L,
         -0.00320912498879085894856L,
         0.000451534528645796438704L,
         -0.933241270357061460782e-5L,
        };
      static const T Q[6] = {    
         1L,
         -0.279960334310344432495L,
         0.0419676223309986037706L,
         -0.00413421406552171059003L,
         0.00024978985622317935355L,
         -0.101855788418564031874e-4L,
      };
      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:        9.007e-20
      // Expected Error Term:            9.007e-20
      static const T P[6] = {    
         0.577215664901532860516,
         0.243210646940107164097,
         0.0417364673988216497593,
         0.00390252087072843288378,
         0.000249606367151877175456,
         0.110108440976732897969e-4,
      };
      static const T Q[6] = {    
         1,
         0.295201277126631761737,
         0.043460910607305495864,
         0.00434930582085826330659,
         0.000255784226140488490982,
         0.10991819782396112081e-4,
      };
      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[6] = {    
         -0.0537258300023595030676,
         0.0445163473292365591906,
         0.0128677673534519952905,
         0.00097541770457391752726,
         0.769875101573654070925e-4,
         0.328032510000383084155e-5,
      };
      static const T Q[7] = {    
         1,
         0.33383194553034051422,
         0.0487798431291407621462,
         0.00479039708573558490716,
         0.000270776703956336357707,
         0.106951867532057341359e-4,
         0.236276623974978646399e-7,
      };
      result = tools::evaluate_polynomial(P, s - 2) / tools::evaluate_polynomial(Q, s - 2);
      result += Y + 1 / (-sc);
   }
   else if(s <= 7)
   {
      // Maximum Deviation Found:                     2.955e-17
      // Expected Error Term:                         2.955e-17
      // Max error found at double precision:         2.009135e-16

      static const T P[6] = {    
         -2.49710190602259410021,
         -2.60013301809475665334,
         -0.939260435377109939261,
         -0.138448617995741530935,
         -0.00701721240549802377623,
         -0.229257310594893932383e-4,
      };
      static const T Q[9] = {    
         1,
         0.706039025937745133628,
         0.15739599649558626358,
         0.0106117950976845084417,
         -0.36910273311764618902e-4,
         0.493409563927590008943e-5,
         -0.234055487025287216506e-6,
         0.718833729365459760664e-8,
         -0.1129200113474947419e-9,
      };
      result = tools::evaluate_polynomial(P, s - 4) / tools::evaluate_polynomial(Q, s - 4);
      result = 1 + exp(result);
   }
   else if(s < 15)
   {
      // Maximum Deviation Found:                     7.117e-16
      // Expected Error Term:                         7.117e-16
      // Max error found at double precision:         9.387771e-16
      static const T P[7] = {    
         -4.78558028495135619286,
         -1.89197364881972536382,
         -0.211407134874412820099,
         -0.000189204758260076688518,
         0.00115140923889178742086,
         0.639949204213164496988e-4,
         0.139348932445324888343e-5,
        };
      static const T Q[9] = {    
         1,
         0.244345337378188557777,
         0.00873370754492288653669,
         -0.00117592765334434471562,
         -0.743743682899933180415e-4,
         -0.21750464515767984778e-5,
         0.471001264003076486547e-8,
         -0.833378440625385520576e-10,
         0.699841545204845636531e-12,
        };
      result = tools::evaluate_polynomial(P, s - 7) / tools::evaluate_polynomial(Q, s - 7);
      result = 1 + exp(result);
   }
   else if(s < 36)
   {
      // Max error in interpolated form:             1.668e-17
      // Max error found at long double precision:   1.669714e-17
      static const T P[8] = {    
         -10.3948950573308896825,
         -2.85827219671106697179,
         -0.347728266539245787271,
         -0.0251156064655346341766,
         -0.00119459173416968685689,
         -0.382529323507967522614e-4,
         -0.785523633796723466968e-6,
         -0.821465709095465524192e-8,
      };
      static const T Q[10] = {    
         1,
         0.208196333572671890965,
         0.0195687657317205033485,
         0.00111079638102485921877,
         0.408507746266039256231e-4,
         0.955561123065693483991e-6,
         0.118507153474022900583e-7,
         0.222609483627352615142e-14,
      };
      result = tools::evaluate_polynomial(P, s - 15) / tools::evaluate_polynomial(Q, s - 15);
      result = 1 + exp(result);
   }