view release on metacpan or  search on metacpan

src/boost/math/special_functions/zeta.hpp  view on Meta::CPAN

   // 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, 
      policies::get_epsilon<T, Policy>(),
      max_iter);
   policies::check_series_iterations<T>("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
   return result;
}
#endif

template <class T, class Policy>
T zeta_polynomial_series(T s, T sc, Policy const &)
{
   //
   // This is algorithm 3 from:
   // 
   // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein, 
   // Canadian Mathematical Society, Conference Proceedings.
   // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
   //
   BOOST_MATH_STD_USING
   int n = itrunc(T(log(boost::math::tools::epsilon<T>()) / -2));
   T sum = 0;
   T two_n = ldexp(T(1), n);
   int ej_sign = 1;
   for(int j = 0; j < n; ++j)
   {
      sum += ej_sign * -two_n / pow(T(j + 1), s);
      ej_sign = -ej_sign; 
   }
   T ej_sum = 1;
   T ej_term = 1;
   for(int j = n; j <= 2 * n - 1; ++j)
   {
      sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s);
      ej_sign = -ej_sign;
      ej_term *= 2 * n - j;
      ej_term /= j - n + 1;
      ej_sum += ej_term;
   }
   return -sum / (two_n * (-powm1(T(2), sc)));
}

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;
   result = zeta_polynomial_series(s, sc, pol); 
#if 0
   // Old code archived for future reference:

   //
   // 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);
#endif
   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,