math/special_functions/zeta.hpp

0001 //  Copyright John Maddock 2007, 2014.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_ZETA_HPP
0007 #define BOOST_MATH_ZETA_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/special_functions/math_fwd.hpp>
0014 #include <boost/math/tools/precision.hpp>
0015 #include <boost/math/tools/series.hpp>
0016 #include <boost/math/tools/big_constant.hpp>
0017 #include <boost/math/policies/error_handling.hpp>
0018 #include <boost/math/special_functions/gamma.hpp>
0019 #include <boost/math/special_functions/factorials.hpp>
0020 #include <boost/math/special_functions/sin_pi.hpp>
0021
0022 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0023 //
0024 // This is the only way we can avoid
0025 // warning: non-standard suffix on floating constant [-Wpedantic]
0026 // when building with -Wall -pedantic.  Neither __extension__
0027 // nor #pragma diagnostic ignored work :(
0028 //
0029 #pragma GCC system_header
0030 #endif
0031
0032 namespace boost{ namespace math{ namespace detail{
0033
0034 #if 0
0035 //
0036 // This code is commented out because we have a better more rapidly converging series
0037 // now.  Retained for future reference and in case the new code causes any issues down the line....
0038 //
0039
0040 template <class T, class Policy>
0041 struct zeta_series_cache_size
0042 {
0043    //
0044    // Work how large to make our cache size when evaluating the series
0045    // evaluation:  normally this is just large enough for the series
0046    // to have converged, but for arbitrary precision types we need a
0047    // really large cache to achieve reasonable precision in a reasonable
0048    // time.  This is important when constructing rational approximations
0049    // to zeta for example.
0050    //
0051    typedef typename boost::math::policies::precision<T,Policy>::type precision_type;
0052    typedef typename mpl::if_<
0053       mpl::less_equal<precision_type, std::integral_constant<int, 0> >,
0054       std::integral_constant<int, 5000>,
0055       typename mpl::if_<
0056          mpl::less_equal<precision_type, std::integral_constant<int, 64> >,
0057          std::integral_constant<int, 70>,
0058          typename mpl::if_<
0059             mpl::less_equal<precision_type, std::integral_constant<int, 113> >,
0060             std::integral_constant<int, 100>,
0061             std::integral_constant<int, 5000>
0062          >::type
0063       >::type
0064    >::type type;
0065 };
0066
0067 template <class T, class Policy>
0068 T zeta_series_imp(T s, T sc, const Policy&)
0069 {
0070    //
0071    // Series evaluation from:
0072    // Havil, J. Gamma: Exploring Euler's Constant.
0073    // Princeton, NJ: Princeton University Press, 2003.
0074    //
0075    // See also http://mathworld.wolfram.com/RiemannZetaFunction.html
0076    //
0077    BOOST_MATH_STD_USING
0078    T sum = 0;
0079    T mult = 0.5;
0080    T change;
0081    typedef typename zeta_series_cache_size<T,Policy>::type cache_size;
0082    T powers[cache_size::value] = { 0, };
0083    unsigned n = 0;
0084    do{
0085       T binom = -static_cast<T>(n);
0086       T nested_sum = 1;
0087       if(n < sizeof(powers) / sizeof(powers[0]))
0088          powers[n] = pow(static_cast<T>(n + 1), -s);
0089       for(unsigned k = 1; k <= n; ++k)
0090       {
0091          T p;
0092          if(k < sizeof(powers) / sizeof(powers[0]))
0093          {
0094             p = powers[k];
0095             //p = pow(k + 1, -s);
0096          }
0097          else
0098             p = pow(static_cast<T>(k + 1), -s);
0099          nested_sum += binom * p;
0100         binom *= (k - static_cast<T>(n)) / (k + 1);
0101       }
0102       change = mult * nested_sum;
0103       sum += change;
0104       mult /= 2;
0105       ++n;
0106    }while(fabs(change / sum) > tools::epsilon<T>());
0107
0108    return sum * 1 / -boost::math::powm1(T(2), sc);
0109 }
0110
0111 //
0112 // Classical p-series:
0113 //
0114 template <class T>
0115 struct zeta_series2
0116 {
0117    typedef T result_type;
0118    zeta_series2(T _s) : s(-_s), k(1){}
0119    T operator()()
0120    {
0121       BOOST_MATH_STD_USING
0122       return pow(static_cast<T>(k++), s);
0123    }
0124 private:
0125    T s;
0126    unsigned k;
0127 };
0128
0129 template <class T, class Policy>
0130 inline T zeta_series2_imp(T s, const Policy& pol)
0131 {
0132    std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();;
0133    zeta_series2<T> f(s);
0134    T result = tools::sum_series(
0135       f,
0136       policies::get_epsilon<T, Policy>(),
0137       max_iter);
0138    policies::check_series_iterations<T>("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
0139    return result;
0140 }
0141 #endif
0142
0143 template <class T, class Policy>
0144 T zeta_polynomial_series(T s, T sc, Policy const &)
0145 {
0146    //
0147    // This is algorithm 3 from:
0148    //
0149    // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein,
0150    // Canadian Mathematical Society, Conference Proceedings.
0151    // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
0152    //
0153    BOOST_MATH_STD_USING
0154    int n = itrunc(T(log(boost::math::tools::epsilon<T>()) / -2));
0155    T sum = 0;
0156    T two_n = ldexp(T(1), n);
0157    int ej_sign = 1;
0158    for(int j = 0; j < n; ++j)
0159    {
0160       sum += ej_sign * -two_n / pow(T(j + 1), s);
0161       ej_sign = -ej_sign;
0162    }
0163    T ej_sum = 1;
0164    T ej_term = 1;
0165    for(int j = n; j <= 2 * n - 1; ++j)
0166    {
0167       sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s);
0168       ej_sign = -ej_sign;
0169       ej_term *= 2 * n - j;
0170       ej_term /= j - n + 1;
0171       ej_sum += ej_term;
0172    }
0173    return -sum / (two_n * (-powm1(T(2), sc)));
0174 }
0175
0176 template <class T, class Policy>
0177 T zeta_imp_prec(T s, T sc, const Policy& pol, const std::integral_constant<int, 0>&)
0178 {
0179    BOOST_MATH_STD_USING
0180    T result;
0181    if(s >= policies::digits<T, Policy>())
0182       return 1;
0183    result = zeta_polynomial_series(s, sc, pol);
0184 #if 0
0185    // Old code archived for future reference:
0186
0187    //
0188    // Only use power series if it will converge in 100
0189    // iterations or less: the more iterations it consumes
0190    // the slower convergence becomes so we have to be very
0191    // careful in it's usage.
0192    //
0193    if (s > -log(tools::epsilon<T>()) / 4.5)
0194       result = detail::zeta_series2_imp(s, pol);
0195    else
0196       result = detail::zeta_series_imp(s, sc, pol);
0197 #endif
0198    return result;
0199 }
0200
0201 template <class T, class Policy>
0202 inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 53>&)
0203 {
0204    BOOST_MATH_STD_USING
0205    T result;
0206    if(s < 1)
0207    {
0208       // Rational Approximation
0209       // Maximum Deviation Found:                     2.020e-18
0210       // Expected Error Term:                         -2.020e-18
0211       // Max error found at double precision:         3.994987e-17
0212       static const T P[6] = {
0213          static_cast<T>(0.24339294433593750202L),
0214          static_cast<T>(-0.49092470516353571651L),
0215          static_cast<T>(0.0557616214776046784287L),
0216          static_cast<T>(-0.00320912498879085894856L),
0217          static_cast<T>(0.000451534528645796438704L),
0218          static_cast<T>(-0.933241270357061460782e-5L),
0219         };
0220       static const T Q[6] = {
0221          static_cast<T>(1L),
0222          static_cast<T>(-0.279960334310344432495L),
0223          static_cast<T>(0.0419676223309986037706L),
0224          static_cast<T>(-0.00413421406552171059003L),
0225          static_cast<T>(0.00024978985622317935355L),
0226          static_cast<T>(-0.101855788418564031874e-4L),
0227       };
0228       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
0229       result -= 1.2433929443359375F;
0230       result += (sc);
0231       result /= (sc);
0232    }
0233    else if(s <= 2)
0234    {
0235       // Maximum Deviation Found:        9.007e-20
0236       // Expected Error Term:            9.007e-20
0237       static const T P[6] = {
0238          static_cast<T>(0.577215664901532860516L),
0239          static_cast<T>(0.243210646940107164097L),
0240          static_cast<T>(0.0417364673988216497593L),
0241          static_cast<T>(0.00390252087072843288378L),
0242          static_cast<T>(0.000249606367151877175456L),
0243          static_cast<T>(0.110108440976732897969e-4L),
0244       };
0245       static const T Q[6] = {
0246          static_cast<T>(1.0),
0247          static_cast<T>(0.295201277126631761737L),
0248          static_cast<T>(0.043460910607305495864L),
0249          static_cast<T>(0.00434930582085826330659L),
0250          static_cast<T>(0.000255784226140488490982L),
0251          static_cast<T>(0.10991819782396112081e-4L),
0252       };
0253       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
0254       result += 1 / (-sc);
0255    }
0256    else if(s <= 4)
0257    {
0258       // Maximum Deviation Found:          5.946e-22
0259       // Expected Error Term:              -5.946e-22
0260       static const float Y = 0.6986598968505859375;
0261       static const T P[6] = {
0262          static_cast<T>(-0.0537258300023595030676L),
0263          static_cast<T>(0.0445163473292365591906L),
0264          static_cast<T>(0.0128677673534519952905L),
0265          static_cast<T>(0.00097541770457391752726L),
0266          static_cast<T>(0.769875101573654070925e-4L),
0267          static_cast<T>(0.328032510000383084155e-5L),
0268       };
0269       static const T Q[7] = {
0270          1.0f,
0271          static_cast<T>(0.33383194553034051422L),
0272          static_cast<T>(0.0487798431291407621462L),
0273          static_cast<T>(0.00479039708573558490716L),
0274          static_cast<T>(0.000270776703956336357707L),
0275          static_cast<T>(0.106951867532057341359e-4L),
0276          static_cast<T>(0.236276623974978646399e-7L),
0277       };
0278       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
0279       result += Y + 1 / (-sc);
0280    }
0281    else if(s <= 7)
0282    {
0283       // Maximum Deviation Found:                     2.955e-17
0284       // Expected Error Term:                         2.955e-17
0285       // Max error found at double precision:         2.009135e-16
0286
0287       static const T P[6] = {
0288          static_cast<T>(-2.49710190602259410021L),
0289          static_cast<T>(-2.60013301809475665334L),
0290          static_cast<T>(-0.939260435377109939261L),
0291          static_cast<T>(-0.138448617995741530935L),
0292          static_cast<T>(-0.00701721240549802377623L),
0293          static_cast<T>(-0.229257310594893932383e-4L),
0294       };
0295       static const T Q[9] = {
0296          1.0f,
0297          static_cast<T>(0.706039025937745133628L),
0298          static_cast<T>(0.15739599649558626358L),
0299          static_cast<T>(0.0106117950976845084417L),
0300          static_cast<T>(-0.36910273311764618902e-4L),
0301          static_cast<T>(0.493409563927590008943e-5L),
0302          static_cast<T>(-0.234055487025287216506e-6L),
0303          static_cast<T>(0.718833729365459760664e-8L),
0304          static_cast<T>(-0.1129200113474947419e-9L),
0305       };
0306       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
0307       result = 1 + exp(result);
0308    }
0309    else if(s < 15)
0310    {
0311       // Maximum Deviation Found:                     7.117e-16
0312       // Expected Error Term:                         7.117e-16
0313       // Max error found at double precision:         9.387771e-16
0314       static const T P[7] = {
0315          static_cast<T>(-4.78558028495135619286L),
0316          static_cast<T>(-1.89197364881972536382L),
0317          static_cast<T>(-0.211407134874412820099L),
0318          static_cast<T>(-0.000189204758260076688518L),
0319          static_cast<T>(0.00115140923889178742086L),
0320          static_cast<T>(0.639949204213164496988e-4L),
0321          static_cast<T>(0.139348932445324888343e-5L),
0322         };
0323       static const T Q[9] = {
0324          1.0f,
0325          static_cast<T>(0.244345337378188557777L),
0326          static_cast<T>(0.00873370754492288653669L),
0327          static_cast<T>(-0.00117592765334434471562L),
0328          static_cast<T>(-0.743743682899933180415e-4L),
0329          static_cast<T>(-0.21750464515767984778e-5L),
0330          static_cast<T>(0.471001264003076486547e-8L),
0331          static_cast<T>(-0.833378440625385520576e-10L),
0332          static_cast<T>(0.699841545204845636531e-12L),
0333         };
0334       result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));
0335       result = 1 + exp(result);
0336    }
0337    else if(s < 36)
0338    {
0339       // Max error in interpolated form:             1.668e-17
0340       // Max error found at long double precision:   1.669714e-17
0341       static const T P[8] = {
0342          static_cast<T>(-10.3948950573308896825L),
0343          static_cast<T>(-2.85827219671106697179L),
0344          static_cast<T>(-0.347728266539245787271L),
0345          static_cast<T>(-0.0251156064655346341766L),
0346          static_cast<T>(-0.00119459173416968685689L),
0347          static_cast<T>(-0.382529323507967522614e-4L),
0348          static_cast<T>(-0.785523633796723466968e-6L),
0349          static_cast<T>(-0.821465709095465524192e-8L),
0350       };
0351       static const T Q[10] = {
0352          1.0f,
0353          static_cast<T>(0.208196333572671890965L),
0354          static_cast<T>(0.0195687657317205033485L),
0355          static_cast<T>(0.00111079638102485921877L),
0356          static_cast<T>(0.408507746266039256231e-4L),
0357          static_cast<T>(0.955561123065693483991e-6L),
0358          static_cast<T>(0.118507153474022900583e-7L),
0359          static_cast<T>(0.222609483627352615142e-14L),
0360       };
0361       result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));
0362       result = 1 + exp(result);
0363    }
0364    else if(s < 56)
0365    {
0366       result = 1 + pow(T(2), -s);
0367    }
0368    else
0369    {
0370       result = 1;
0371    }
0372    return result;
0373 }
0374
0375 template <class T, class Policy>
0376 T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 64>&)
0377 {
0378    BOOST_MATH_STD_USING
0379    T result;
0380    if(s < 1)
0381    {
0382       // Rational Approximation
0383       // Maximum Deviation Found:                     3.099e-20
0384       // Expected Error Term:                         3.099e-20
0385       // Max error found at long double precision:    5.890498e-20
0386       static const T P[6] = {
0387          BOOST_MATH_BIG_CONSTANT(T, 64, 0.243392944335937499969),
0388          BOOST_MATH_BIG_CONSTANT(T, 64, -0.496837806864865688082),
0389          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0680008039723709987107),
0390          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00511620413006619942112),
0391          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000455369899250053003335),
0392          BOOST_MATH_BIG_CONSTANT(T, 64, -0.279496685273033761927e-4),
0393         };
0394       static const T Q[7] = {
0395          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0396          BOOST_MATH_BIG_CONSTANT(T, 64, -0.30425480068225790522),
0397          BOOST_MATH_BIG_CONSTANT(T, 64, 0.050052748580371598736),
0398          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00519355671064700627862),
0399          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000360623385771198350257),
0400          BOOST_MATH_BIG_CONSTANT(T, 64, -0.159600883054550987633e-4),
0401          BOOST_MATH_BIG_CONSTANT(T, 64, 0.339770279812410586032e-6),
0402       };
0403       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
0404       result -= 1.2433929443359375F;
0405       result += (sc);
0406       result /= (sc);
0407    }
0408    else if(s <= 2)
0409    {
0410       // Maximum Deviation Found:                     1.059e-21
0411       // Expected Error Term:                         1.059e-21
0412       // Max error found at long double precision:    1.626303e-19
0413
0414       static const T P[6] = {
0415          BOOST_MATH_BIG_CONSTANT(T, 64, 0.577215664901532860605),
0416          BOOST_MATH_BIG_CONSTANT(T, 64, 0.222537368917162139445),
0417          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0356286324033215682729),
0418          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00304465292366350081446),
0419          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000178102511649069421904),
0420          BOOST_MATH_BIG_CONSTANT(T, 64, 0.700867470265983665042e-5),
0421       };
0422       static const T Q[7] = {
0423          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0424          BOOST_MATH_BIG_CONSTANT(T, 64, 0.259385759149531030085),
0425          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0373974962106091316854),
0426          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00332735159183332820617),
0427          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000188690420706998606469),
0428          BOOST_MATH_BIG_CONSTANT(T, 64, 0.635994377921861930071e-5),
0429          BOOST_MATH_BIG_CONSTANT(T, 64, 0.226583954978371199405e-7),
0430       };
0431       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
0432       result += 1 / (-sc);
0433    }
0434    else if(s <= 4)
0435    {
0436       // Maximum Deviation Found:          5.946e-22
0437       // Expected Error Term:              -5.946e-22
0438       static const float Y = 0.6986598968505859375;
0439       static const T P[7] = {
0440          BOOST_MATH_BIG_CONSTANT(T, 64, -0.053725830002359501027),
0441          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0470551187571475844778),
0442          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0101339410415759517471),
0443          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00100240326666092854528),
0444          BOOST_MATH_BIG_CONSTANT(T, 64, 0.685027119098122814867e-4),
0445          BOOST_MATH_BIG_CONSTANT(T, 64, 0.390972820219765942117e-5),
0446          BOOST_MATH_BIG_CONSTANT(T, 64, 0.540319769113543934483e-7),
0447       };
0448       static const T Q[8] = {
0449          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0450          BOOST_MATH_BIG_CONSTANT(T, 64, 0.286577739726542730421),
0451          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0447355811517733225843),
0452          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00430125107610252363302),
0453          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000284956969089786662045),
0454          BOOST_MATH_BIG_CONSTANT(T, 64, 0.116188101609848411329e-4),
0455          BOOST_MATH_BIG_CONSTANT(T, 64, 0.278090318191657278204e-6),
0456          BOOST_MATH_BIG_CONSTANT(T, 64, -0.19683620233222028478e-8),
0457       };
0458       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
0459       result += Y + 1 / (-sc);
0460    }
0461    else if(s <= 7)
0462    {
0463       // Max error found at long double precision: 8.132216e-19
0464       static const T P[8] = {
0465          BOOST_MATH_BIG_CONSTANT(T, 64, -2.49710190602259407065),
0466          BOOST_MATH_BIG_CONSTANT(T, 64, -3.36664913245960625334),
0467          BOOST_MATH_BIG_CONSTANT(T, 64, -1.77180020623777595452),
0468          BOOST_MATH_BIG_CONSTANT(T, 64, -0.464717885249654313933),
0469          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0643694921293579472583),
0470          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00464265386202805715487),
0471          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000165556579779704340166),
0472          BOOST_MATH_BIG_CONSTANT(T, 64, -0.252884970740994069582e-5),
0473       };
0474       static const T Q[9] = {
0475          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0476          BOOST_MATH_BIG_CONSTANT(T, 64, 1.01300131390690459085),
0477          BOOST_MATH_BIG_CONSTANT(T, 64, 0.387898115758643503827),
0478          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0695071490045701135188),
0479          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00586908595251442839291),
0480          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000217752974064612188616),
0481          BOOST_MATH_BIG_CONSTANT(T, 64, 0.397626583349419011731e-5),
0482          BOOST_MATH_BIG_CONSTANT(T, 64, -0.927884739284359700764e-8),
0483          BOOST_MATH_BIG_CONSTANT(T, 64, 0.119810501805618894381e-9),
0484       };
0485       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
0486       result = 1 + exp(result);
0487    }
0488    else if(s < 15)
0489    {
0490       // Max error in interpolated form:              1.133e-18
0491       // Max error found at long double precision:    2.183198e-18
0492       static const T P[9] = {
0493          BOOST_MATH_BIG_CONSTANT(T, 64, -4.78558028495135548083),
0494          BOOST_MATH_BIG_CONSTANT(T, 64, -3.23873322238609358947),
0495          BOOST_MATH_BIG_CONSTANT(T, 64, -0.892338582881021799922),
0496          BOOST_MATH_BIG_CONSTANT(T, 64, -0.131326296217965913809),
0497          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0115651591773783712996),
0498          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000657728968362695775205),
0499          BOOST_MATH_BIG_CONSTANT(T, 64, -0.252051328129449973047e-4),
0500          BOOST_MATH_BIG_CONSTANT(T, 64, -0.626503445372641798925e-6),
0501          BOOST_MATH_BIG_CONSTANT(T, 64, -0.815696314790853893484e-8),
0502         };
0503       static const T Q[9] = {
0504          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0505          BOOST_MATH_BIG_CONSTANT(T, 64, 0.525765665400123515036),
0506          BOOST_MATH_BIG_CONSTANT(T, 64, 0.10852641753657122787),
0507          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0115669945375362045249),
0508          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000732896513858274091966),
0509          BOOST_MATH_BIG_CONSTANT(T, 64, 0.30683952282420248448e-4),
0510          BOOST_MATH_BIG_CONSTANT(T, 64, 0.819649214609633126119e-6),
0511          BOOST_MATH_BIG_CONSTANT(T, 64, 0.117957556472335968146e-7),
0512          BOOST_MATH_BIG_CONSTANT(T, 64, -0.193432300973017671137e-12),
0513         };
0514       result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));
0515       result = 1 + exp(result);
0516    }
0517    else if(s < 42)
0518    {
0519       // Max error in interpolated form:             1.668e-17
0520       // Max error found at long double precision:   1.669714e-17
0521       static const T P[9] = {
0522          BOOST_MATH_BIG_CONSTANT(T, 64, -10.3948950573308861781),
0523          BOOST_MATH_BIG_CONSTANT(T, 64, -2.82646012777913950108),
0524          BOOST_MATH_BIG_CONSTANT(T, 64, -0.342144362739570333665),
0525          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0249285145498722647472),
0526          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00122493108848097114118),
0527          BOOST_MATH_BIG_CONSTANT(T, 64, -0.423055371192592850196e-4),
0528          BOOST_MATH_BIG_CONSTANT(T, 64, -0.1025215577185967488e-5),
0529          BOOST_MATH_BIG_CONSTANT(T, 64, -0.165096762663509467061e-7),
0530          BOOST_MATH_BIG_CONSTANT(T, 64, -0.145392555873022044329e-9),
0531       };
0532       static const T Q[10] = {
0533          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0534          BOOST_MATH_BIG_CONSTANT(T, 64, 0.205135978585281988052),
0535          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0192359357875879453602),
0536          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00111496452029715514119),
0537          BOOST_MATH_BIG_CONSTANT(T, 64, 0.434928449016693986857e-4),
0538          BOOST_MATH_BIG_CONSTANT(T, 64, 0.116911068726610725891e-5),
0539          BOOST_MATH_BIG_CONSTANT(T, 64, 0.206704342290235237475e-7),
0540          BOOST_MATH_BIG_CONSTANT(T, 64, 0.209772836100827647474e-9),
0541          BOOST_MATH_BIG_CONSTANT(T, 64, -0.939798249922234703384e-16),
0542          BOOST_MATH_BIG_CONSTANT(T, 64, 0.264584017421245080294e-18),
0543       };
0544       result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));
0545       result = 1 + exp(result);
0546    }
0547    else if(s < 63)
0548    {
0549       result = 1 + pow(T(2), -s);
0550    }
0551    else
0552    {
0553       result = 1;
0554    }
0555    return result;
0556 }
0557
0558 template <class T, class Policy>
0559 T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 113>&)
0560 {
0561    BOOST_MATH_STD_USING
0562    T result;
0563    if(s < 1)
0564    {
0565       // Rational Approximation
0566       // Maximum Deviation Found:                     9.493e-37
0567       // Expected Error Term:                         9.492e-37
0568       // Max error found at long double precision:    7.281332e-31
0569
0570       static const T P[10] = {
0571          BOOST_MATH_BIG_CONSTANT(T, 113, -1.0),
0572          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0353008629988648122808504280990313668),
0573          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0107795651204927743049369868548706909),
0574          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000523961870530500751114866884685172975),
0575          BOOST_MATH_BIG_CONSTANT(T, 113, -0.661805838304910731947595897966487515e-4),
0576          BOOST_MATH_BIG_CONSTANT(T, 113, -0.658932670403818558510656304189164638e-5),
0577          BOOST_MATH_BIG_CONSTANT(T, 113, -0.103437265642266106533814021041010453e-6),
0578          BOOST_MATH_BIG_CONSTANT(T, 113, 0.116818787212666457105375746642927737e-7),
0579          BOOST_MATH_BIG_CONSTANT(T, 113, 0.660690993901506912123512551294239036e-9),
0580          BOOST_MATH_BIG_CONSTANT(T, 113, 0.113103113698388531428914333768142527e-10),
0581         };
0582       static const T Q[11] = {
0583          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0584          BOOST_MATH_BIG_CONSTANT(T, 113, -0.387483472099602327112637481818565459),
0585          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0802265315091063135271497708694776875),
0586          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0110727276164171919280036408995078164),
0587          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00112552716946286252000434849173787243),
0588          BOOST_MATH_BIG_CONSTANT(T, 113, -0.874554160748626916455655180296834352e-4),
0589          BOOST_MATH_BIG_CONSTANT(T, 113, 0.530097847491828379568636739662278322e-5),
0590          BOOST_MATH_BIG_CONSTANT(T, 113, -0.248461553590496154705565904497247452e-6),
0591          BOOST_MATH_BIG_CONSTANT(T, 113, 0.881834921354014787309644951507523899e-8),
0592          BOOST_MATH_BIG_CONSTANT(T, 113, -0.217062446168217797598596496310953025e-9),
0593          BOOST_MATH_BIG_CONSTANT(T, 113, 0.315823200002384492377987848307151168e-11),
0594       };
0595       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
0596       result += (sc);
0597       result /= (sc);
0598    }
0599    else if(s <= 2)
0600    {
0601       // Maximum Deviation Found:                     1.616e-37
0602       // Expected Error Term:                         -1.615e-37
0603
0604       static const T P[10] = {
0605          BOOST_MATH_BIG_CONSTANT(T, 113, 0.577215664901532860606512090082402431),
0606          BOOST_MATH_BIG_CONSTANT(T, 113, 0.255597968739771510415479842335906308),
0607          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0494056503552807274142218876983542205),
0608          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00551372778611700965268920983472292325),
0609          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00043667616723970574871427830895192731),
0610          BOOST_MATH_BIG_CONSTANT(T, 113, 0.268562259154821957743669387915239528e-4),
0611          BOOST_MATH_BIG_CONSTANT(T, 113, 0.109249633923016310141743084480436612e-5),
0612          BOOST_MATH_BIG_CONSTANT(T, 113, 0.273895554345300227466534378753023924e-7),
0613          BOOST_MATH_BIG_CONSTANT(T, 113, 0.583103205551702720149237384027795038e-9),
0614          BOOST_MATH_BIG_CONSTANT(T, 113, -0.835774625259919268768735944711219256e-11),
0615       };
0616       static const T Q[11] = {
0617          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0618          BOOST_MATH_BIG_CONSTANT(T, 113, 0.316661751179735502065583176348292881),
0619          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0540401806533507064453851182728635272),
0620          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00598621274107420237785899476374043797),
0621          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000474907812321704156213038740142079615),
0622          BOOST_MATH_BIG_CONSTANT(T, 113, 0.272125421722314389581695715835862418e-4),
0623          BOOST_MATH_BIG_CONSTANT(T, 113, 0.112649552156479800925522445229212933e-5),
0624          BOOST_MATH_BIG_CONSTANT(T, 113, 0.301838975502992622733000078063330461e-7),
0625          BOOST_MATH_BIG_CONSTANT(T, 113, 0.422960728687211282539769943184270106e-9),
0626          BOOST_MATH_BIG_CONSTANT(T, 113, -0.377105263588822468076813329270698909e-11),
0627          BOOST_MATH_BIG_CONSTANT(T, 113, -0.581926559304525152432462127383600681e-13),
0628       };
0629       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
0630       result += 1 / (-sc);
0631    }
0632    else if(s <= 4)
0633    {
0634       // Maximum Deviation Found:                     1.891e-36
0635       // Expected Error Term:                         -1.891e-36
0636       // Max error found: 2.171527e-35
0637
0638       static const float Y = 0.6986598968505859375;
0639       static const T P[11] = {
0640          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0537258300023595010275848333539748089),
0641          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0429086930802630159457448174466342553),
0642          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0136148228754303412510213395034056857),
0643          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00190231601036042925183751238033763915),
0644          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000186880390916311438818302549192456581),
0645          BOOST_MATH_BIG_CONSTANT(T, 113, 0.145347370745893262394287982691323657e-4),
0646          BOOST_MATH_BIG_CONSTANT(T, 113, 0.805843276446813106414036600485884885e-6),
0647          BOOST_MATH_BIG_CONSTANT(T, 113, 0.340818159286739137503297172091882574e-7),
0648          BOOST_MATH_BIG_CONSTANT(T, 113, 0.115762357488748996526167305116837246e-8),
0649          BOOST_MATH_BIG_CONSTANT(T, 113, 0.231904754577648077579913403645767214e-10),
0650          BOOST_MATH_BIG_CONSTANT(T, 113, 0.340169592866058506675897646629036044e-12),
0651       };
0652       static const T Q[12] = {
0653          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0654          BOOST_MATH_BIG_CONSTANT(T, 113, 0.363755247765087100018556983050520554),
0655          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0696581979014242539385695131258321598),
0656          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00882208914484611029571547753782014817),
0657          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000815405623261946661762236085660996718),
0658          BOOST_MATH_BIG_CONSTANT(T, 113, 0.571366167062457197282642344940445452e-4),
0659          BOOST_MATH_BIG_CONSTANT(T, 113, 0.309278269271853502353954062051797838e-5),
0660          BOOST_MATH_BIG_CONSTANT(T, 113, 0.12822982083479010834070516053794262e-6),
0661          BOOST_MATH_BIG_CONSTANT(T, 113, 0.397876357325018976733953479182110033e-8),
0662          BOOST_MATH_BIG_CONSTANT(T, 113, 0.8484432107648683277598472295289279e-10),
0663          BOOST_MATH_BIG_CONSTANT(T, 113, 0.105677416606909614301995218444080615e-11),
0664          BOOST_MATH_BIG_CONSTANT(T, 113, 0.547223964564003701979951154093005354e-15),
0665       };
0666       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
0667       result += Y + 1 / (-sc);
0668    }
0669    else if(s <= 6)
0670    {
0671       // Max error in interpolated form:             1.510e-37
0672       // Max error found at long double precision:   2.769266e-34
0673
0674       static const T Y = 3.28348541259765625F;
0675
0676       static const T P[13] = {
0677          BOOST_MATH_BIG_CONSTANT(T, 113, 0.786383506575062179339611614117697622),
0678          BOOST_MATH_BIG_CONSTANT(T, 113, 0.495766593395271370974685959652073976),
0679          BOOST_MATH_BIG_CONSTANT(T, 113, -0.409116737851754766422360889037532228),
0680          BOOST_MATH_BIG_CONSTANT(T, 113, -0.57340744006238263817895456842655987),
0681          BOOST_MATH_BIG_CONSTANT(T, 113, -0.280479899797421910694892949057963111),
0682          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0753148409447590257157585696212649869),
0683          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0122934003684672788499099362823748632),
0684          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126148398446193639247961370266962927),
0685          BOOST_MATH_BIG_CONSTANT(T, 113, -0.828465038179772939844657040917364896e-4),
0686          BOOST_MATH_BIG_CONSTANT(T, 113, -0.361008916706050977143208468690645684e-5),
0687          BOOST_MATH_BIG_CONSTANT(T, 113, -0.109879825497910544424797771195928112e-6),
0688          BOOST_MATH_BIG_CONSTANT(T, 113, -0.214539416789686920918063075528797059e-8),
0689          BOOST_MATH_BIG_CONSTANT(T, 113, -0.15090220092460596872172844424267351e-10),
0690       };
0691       static const T Q[14] = {
0692          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0693          BOOST_MATH_BIG_CONSTANT(T, 113, 1.69490865837142338462982225731926485),
0694          BOOST_MATH_BIG_CONSTANT(T, 113, 1.22697696630994080733321401255942464),
0695          BOOST_MATH_BIG_CONSTANT(T, 113, 0.495409420862526540074366618006341533),
0696          BOOST_MATH_BIG_CONSTANT(T, 113, 0.122368084916843823462872905024259633),
0697          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0191412993625268971656513890888208623),
0698          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00191401538628980617753082598351559642),
0699          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000123318142456272424148930280876444459),
0700          BOOST_MATH_BIG_CONSTANT(T, 113, 0.531945488232526067889835342277595709e-5),
0701          BOOST_MATH_BIG_CONSTANT(T, 113, 0.161843184071894368337068779669116236e-6),
0702          BOOST_MATH_BIG_CONSTANT(T, 113, 0.305796079600152506743828859577462778e-8),
0703          BOOST_MATH_BIG_CONSTANT(T, 113, 0.233582592298450202680170811044408894e-10),
0704          BOOST_MATH_BIG_CONSTANT(T, 113, -0.275363878344548055574209713637734269e-13),
0705          BOOST_MATH_BIG_CONSTANT(T, 113, 0.221564186807357535475441900517843892e-15),
0706       };
0707       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
0708       result -= Y;
0709       result = 1 + exp(result);
0710    }
0711    else if(s < 10)
0712    {
0713       // Max error in interpolated form:             1.999e-34
0714       // Max error found at long double precision:   2.156186e-33
0715
0716       static const T P[13] = {
0717          BOOST_MATH_BIG_CONSTANT(T, 113, -4.0545627381873738086704293881227365),
0718          BOOST_MATH_BIG_CONSTANT(T, 113, -4.70088348734699134347906176097717782),
0719          BOOST_MATH_BIG_CONSTANT(T, 113, -2.36921550900925512951976617607678789),
0720          BOOST_MATH_BIG_CONSTANT(T, 113, -0.684322583796369508367726293719322866),
0721          BOOST_MATH_BIG_CONSTANT(T, 113, -0.126026534540165129870721937592996324),
0722          BOOST_MATH_BIG_CONSTANT(T, 113, -0.015636903921778316147260572008619549),
0723          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00135442294754728549644376325814460807),
0724          BOOST_MATH_BIG_CONSTANT(T, 113, -0.842793965853572134365031384646117061e-4),
0725          BOOST_MATH_BIG_CONSTANT(T, 113, -0.385602133791111663372015460784978351e-5),
0726          BOOST_MATH_BIG_CONSTANT(T, 113, -0.130458500394692067189883214401478539e-6),
0727          BOOST_MATH_BIG_CONSTANT(T, 113, -0.315861074947230418778143153383660035e-8),
0728          BOOST_MATH_BIG_CONSTANT(T, 113, -0.500334720512030826996373077844707164e-10),
0729          BOOST_MATH_BIG_CONSTANT(T, 113, -0.420204769185233365849253969097184005e-12),
0730         };
0731       static const T Q[14] = {
0732          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0733          BOOST_MATH_BIG_CONSTANT(T, 113, 0.97663511666410096104783358493318814),
0734          BOOST_MATH_BIG_CONSTANT(T, 113, 0.40878780231201806504987368939673249),
0735          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0963890666609396058945084107597727252),
0736          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0142207619090854604824116070866614505),
0737          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139010220902667918476773423995750877),
0738          BOOST_MATH_BIG_CONSTANT(T, 113, 0.940669540194694997889636696089994734e-4),
0739          BOOST_MATH_BIG_CONSTANT(T, 113, 0.458220848507517004399292480807026602e-5),
0740          BOOST_MATH_BIG_CONSTANT(T, 113, 0.16345521617741789012782420625435495e-6),
0741          BOOST_MATH_BIG_CONSTANT(T, 113, 0.414007452533083304371566316901024114e-8),
0742          BOOST_MATH_BIG_CONSTANT(T, 113, 0.68701473543366328016953742622661377e-10),
0743          BOOST_MATH_BIG_CONSTANT(T, 113, 0.603461891080716585087883971886075863e-12),
0744          BOOST_MATH_BIG_CONSTANT(T, 113, 0.294670713571839023181857795866134957e-16),
0745          BOOST_MATH_BIG_CONSTANT(T, 113, -0.147003914536437243143096875069813451e-18),
0746         };
0747       result = tools::evaluate_polynomial(P, T(s - 6)) / tools::evaluate_polynomial(Q, T(s - 6));
0748       result = 1 + exp(result);
0749    }
0750    else if(s < 17)
0751    {
0752       // Max error in interpolated form:             1.641e-32
0753       // Max error found at long double precision:   1.696121e-32
0754       static const T P[13] = {
0755          BOOST_MATH_BIG_CONSTANT(T, 113, -6.91319491921722925920883787894829678),
0756          BOOST_MATH_BIG_CONSTANT(T, 113, -3.65491257639481960248690596951049048),
0757          BOOST_MATH_BIG_CONSTANT(T, 113, -0.813557553449954526442644544105257881),
0758          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0994317301685870959473658713841138083),
0759          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00726896610245676520248617014211734906),
0760          BOOST_MATH_BIG_CONSTANT(T, 113, -0.000317253318715075854811266230916762929),
0761          BOOST_MATH_BIG_CONSTANT(T, 113, -0.66851422826636750855184211580127133e-5),
0762          BOOST_MATH_BIG_CONSTANT(T, 113, 0.879464154730985406003332577806849971e-7),
0763          BOOST_MATH_BIG_CONSTANT(T, 113, 0.113838903158254250631678791998294628e-7),
0764          BOOST_MATH_BIG_CONSTANT(T, 113, 0.379184410304927316385211327537817583e-9),
0765          BOOST_MATH_BIG_CONSTANT(T, 113, 0.612992858643904887150527613446403867e-11),
0766          BOOST_MATH_BIG_CONSTANT(T, 113, 0.347873737198164757035457841688594788e-13),
0767          BOOST_MATH_BIG_CONSTANT(T, 113, -0.289187187441625868404494665572279364e-15),
0768         };
0769       static const T Q[14] = {
0770          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0771          BOOST_MATH_BIG_CONSTANT(T, 113, 0.427310044448071818775721584949868806),
0772          BOOST_MATH_BIG_CONSTANT(T, 113, 0.074602514873055756201435421385243062),
0773          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00688651562174480772901425121653945942),
0774          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000360174847635115036351323894321880445),
0775          BOOST_MATH_BIG_CONSTANT(T, 113, 0.973556847713307543918865405758248777e-5),
0776          BOOST_MATH_BIG_CONSTANT(T, 113, -0.853455848314516117964634714780874197e-8),
0777          BOOST_MATH_BIG_CONSTANT(T, 113, -0.118203513654855112421673192194622826e-7),
0778          BOOST_MATH_BIG_CONSTANT(T, 113, -0.462521662511754117095006543363328159e-9),
0779          BOOST_MATH_BIG_CONSTANT(T, 113, -0.834212591919475633107355719369463143e-11),
0780          BOOST_MATH_BIG_CONSTANT(T, 113, -0.5354594751002702935740220218582929e-13),
0781          BOOST_MATH_BIG_CONSTANT(T, 113, 0.406451690742991192964889603000756203e-15),
0782          BOOST_MATH_BIG_CONSTANT(T, 113, 0.887948682401000153828241615760146728e-19),
0783          BOOST_MATH_BIG_CONSTANT(T, 113, -0.34980761098820347103967203948619072e-21),
0784         };
0785       result = tools::evaluate_polynomial(P, T(s - 10)) / tools::evaluate_polynomial(Q, T(s - 10));
0786       result = 1 + exp(result);
0787    }
0788    else if(s < 30)
0789    {
0790       // Max error in interpolated form:             1.563e-31
0791       // Max error found at long double precision:   1.562725e-31
0792
0793       static const T P[13] = {
0794          BOOST_MATH_BIG_CONSTANT(T, 113, -11.7824798233959252791987402769438322),
0795          BOOST_MATH_BIG_CONSTANT(T, 113, -4.36131215284987731928174218354118102),
0796          BOOST_MATH_BIG_CONSTANT(T, 113, -0.732260980060982349410898496846972204),
0797          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0744985185694913074484248803015717388),
0798          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00517228281320594683022294996292250527),
0799          BOOST_MATH_BIG_CONSTANT(T, 113, -0.000260897206152101522569969046299309939),
0800          BOOST_MATH_BIG_CONSTANT(T, 113, -0.989553462123121764865178453128769948e-5),
0801          BOOST_MATH_BIG_CONSTANT(T, 113, -0.286916799741891410827712096608826167e-6),
0802          BOOST_MATH_BIG_CONSTANT(T, 113, -0.637262477796046963617949532211619729e-8),
0803          BOOST_MATH_BIG_CONSTANT(T, 113, -0.106796831465628373325491288787760494e-9),
0804          BOOST_MATH_BIG_CONSTANT(T, 113, -0.129343095511091870860498356205376823e-11),
0805          BOOST_MATH_BIG_CONSTANT(T, 113, -0.102397936697965977221267881716672084e-13),
0806          BOOST_MATH_BIG_CONSTANT(T, 113, -0.402663128248642002351627980255756363e-16),
0807       };
0808       static const T Q[14] = {
0809          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0810          BOOST_MATH_BIG_CONSTANT(T, 113, 0.311288325355705609096155335186466508),
0811          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0438318468940415543546769437752132748),
0812          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00374396349183199548610264222242269536),
0813          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000218707451200585197339671707189281302),
0814          BOOST_MATH_BIG_CONSTANT(T, 113, 0.927578767487930747532953583797351219e-5),
0815          BOOST_MATH_BIG_CONSTANT(T, 113, 0.294145760625753561951137473484889639e-6),
0816          BOOST_MATH_BIG_CONSTANT(T, 113, 0.704618586690874460082739479535985395e-8),
0817          BOOST_MATH_BIG_CONSTANT(T, 113, 0.126333332872897336219649130062221257e-9),
0818          BOOST_MATH_BIG_CONSTANT(T, 113, 0.16317315713773503718315435769352765e-11),
0819          BOOST_MATH_BIG_CONSTANT(T, 113, 0.137846712823719515148344938160275695e-13),
0820          BOOST_MATH_BIG_CONSTANT(T, 113, 0.580975420554224366450994232723910583e-16),
0821          BOOST_MATH_BIG_CONSTANT(T, 113, -0.291354445847552426900293580511392459e-22),
0822          BOOST_MATH_BIG_CONSTANT(T, 113, 0.73614324724785855925025452085443636e-25),
0823       };
0824       result = tools::evaluate_polynomial(P, T(s - 17)) / tools::evaluate_polynomial(Q, T(s - 17));
0825       result = 1 + exp(result);
0826    }
0827    else if(s < 74)
0828    {
0829       // Max error in interpolated form:             2.311e-27
0830       // Max error found at long double precision:   2.297544e-27
0831       static const T P[14] = {
0832          BOOST_MATH_BIG_CONSTANT(T, 113, -20.7944102007844314586649688802236072),
0833          BOOST_MATH_BIG_CONSTANT(T, 113, -4.95759941987499442499908748130192187),
0834          BOOST_MATH_BIG_CONSTANT(T, 113, -0.563290752832461751889194629200298688),
0835          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0406197001137935911912457120706122877),
0836          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0020846534789473022216888863613422293),
0837          BOOST_MATH_BIG_CONSTANT(T, 113, -0.808095978462109173749395599401375667e-4),
0838          BOOST_MATH_BIG_CONSTANT(T, 113, -0.244706022206249301640890603610060959e-5),
0839          BOOST_MATH_BIG_CONSTANT(T, 113, -0.589477682919645930544382616501666572e-7),
0840          BOOST_MATH_BIG_CONSTANT(T, 113, -0.113699573675553496343617442433027672e-8),
0841          BOOST_MATH_BIG_CONSTANT(T, 113, -0.174767860183598149649901223128011828e-10),
0842          BOOST_MATH_BIG_CONSTANT(T, 113, -0.210051620306761367764549971980026474e-12),
0843          BOOST_MATH_BIG_CONSTANT(T, 113, -0.189187969537370950337212675466400599e-14),
0844          BOOST_MATH_BIG_CONSTANT(T, 113, -0.116313253429564048145641663778121898e-16),
0845          BOOST_MATH_BIG_CONSTANT(T, 113, -0.376708747782400769427057630528578187e-19),
0846       };
0847       static const T Q[16] = {
0848          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0849          BOOST_MATH_BIG_CONSTANT(T, 113, 0.205076752981410805177554569784219717),
0850          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0202526722696670378999575738524540269),
0851          BOOST_MATH_BIG_CONSTANT(T, 113, 0.001278305290005994980069466658219057),
0852          BOOST_MATH_BIG_CONSTANT(T, 113, 0.576404779858501791742255670403304787e-4),
0853          BOOST_MATH_BIG_CONSTANT(T, 113, 0.196477049872253010859712483984252067e-5),
0854          BOOST_MATH_BIG_CONSTANT(T, 113, 0.521863830500876189501054079974475762e-7),
0855          BOOST_MATH_BIG_CONSTANT(T, 113, 0.109524209196868135198775445228552059e-8),
0856          BOOST_MATH_BIG_CONSTANT(T, 113, 0.181698713448644481083966260949267825e-10),
0857          BOOST_MATH_BIG_CONSTANT(T, 113, 0.234793316975091282090312036524695562e-12),
0858          BOOST_MATH_BIG_CONSTANT(T, 113, 0.227490441461460571047545264251399048e-14),
0859          BOOST_MATH_BIG_CONSTANT(T, 113, 0.151500292036937400913870642638520668e-16),
0860          BOOST_MATH_BIG_CONSTANT(T, 113, 0.543475775154780935815530649335936121e-19),
0861          BOOST_MATH_BIG_CONSTANT(T, 113, 0.241647013434111434636554455083309352e-28),
0862          BOOST_MATH_BIG_CONSTANT(T, 113, -0.557103423021951053707162364713587374e-31),
0863          BOOST_MATH_BIG_CONSTANT(T, 113, 0.618708773442584843384712258199645166e-34),
0864       };
0865       result = tools::evaluate_polynomial(P, T(s - 30)) / tools::evaluate_polynomial(Q, T(s - 30));
0866       result = 1 + exp(result);
0867    }
0868    else if(s < 117)
0869    {
0870       result = 1 + pow(T(2), -s);
0871    }
0872    else
0873    {
0874       result = 1;
0875    }
0876    return result;
0877 }
0878
0879 template <class T, class Policy>
0880 T zeta_imp_odd_integer(int s, const T&, const Policy&, const std::true_type&)
0881 {
0882    static const T results[] = {
0883       BOOST_MATH_BIG_CONSTANT(T, 113, 1.2020569031595942853997381615114500), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0369277551433699263313654864570342), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0083492773819228268397975498497968), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0020083928260822144178527692324121), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0004941886041194645587022825264699), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0001227133475784891467518365263574), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000305882363070204935517285106451), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000076371976378997622736002935630), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000019082127165539389256569577951), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000004769329867878064631167196044), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000001192199259653110730677887189), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000298035035146522801860637051), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000074507117898354294919810042), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000018626597235130490064039099), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000004656629065033784072989233), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000001164155017270051977592974), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000291038504449709968692943), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000072759598350574810145209), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000018189896503070659475848), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000004547473783042154026799), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000001136868407680227849349), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000284217097688930185546), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000071054273952108527129), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000017763568435791203275), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000004440892103143813364), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000001110223025141066134), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000277555756213612417), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000069388939045441537), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000017347234760475766), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000004336808690020650), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000001084202172494241), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000271050543122347), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000067762635780452), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000016940658945098), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000004235164736273), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000001058791184068), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000264697796017), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000066174449004), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000016543612251), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000004135903063), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000001033975766), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000258493941), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000064623485), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000016155871), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000004038968), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000001009742), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000252435), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000063109), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000015777), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000003944), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000986), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000247), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000062), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000015), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000004), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000001),
0884    };
0885    return s > 113 ? 1 : results[(s - 3) / 2];
0886 }
0887
0888 template <class T, class Policy>
0889 T zeta_imp_odd_integer(int s, const T& sc, const Policy& pol, const std::false_type&)
0890 {
0891 #ifdef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES
0892    static_assert(std::is_trivially_destructible<T>::value, "Your platform does not support thread_local with non-trivial types, last checked with Mingw-x64-8.1, Jan 2021.  Please try a Mingw build with the POSIX threading model, see https://sourceforge.net/p/mingw-w64/bugs/527/");
0893 #endif
0894    static BOOST_MATH_THREAD_LOCAL bool is_init = false;
0895    static BOOST_MATH_THREAD_LOCAL T results[50] = {};
0896    static BOOST_MATH_THREAD_LOCAL int digits = tools::digits<T>();
0897    int current_digits = tools::digits<T>();
0898    if(digits != current_digits)
0899    {
0900       // Oh my precision has changed...
0901       is_init = false;
0902    }
0903    if(!is_init)
0904    {
0905       is_init = true;
0906       digits = current_digits;
0907       for(unsigned k = 0; k < sizeof(results) / sizeof(results[0]); ++k)
0908       {
0909          T arg = k * 2 + 3;
0910          T c_arg = 1 - arg;
0911          results[k] = zeta_polynomial_series(arg, c_arg, pol);
0912       }
0913    }
0914    unsigned index = (s - 3) / 2;
0915    return index >= sizeof(results) / sizeof(results[0]) ? zeta_polynomial_series(T(s), sc, pol): results[index];
0916 }
0917
0918 template <class T, class Policy, class Tag>
0919 T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
0920 {
0921    BOOST_MATH_STD_USING
0922    static const char* function = "boost::math::zeta<%1%>";
0923    if(sc == 0)
0924       return policies::raise_pole_error<T>(
0925          function,
0926          "Evaluation of zeta function at pole %1%",
0927          s, pol);
0928    T result;
0929    //
0930    // Trivial case:
0931    //
0932    if(s > policies::digits<T, Policy>())
0933       return 1;
0934    //
0935    // Start by seeing if we have a simple closed form:
0936    //
0937    if(floor(s) == s)
0938    {
0939 #ifndef BOOST_NO_EXCEPTIONS
0940       // Without exceptions we expect itrunc to return INT_MAX on overflow
0941       // and we fall through anyway.
0942       try
0943       {
0944 #endif
0945          int v = itrunc(s);
0946          if(v == s)
0947          {
0948             if(v < 0)
0949             {
0950                if(((-v) & 1) == 0)
0951                   return 0;
0952                int n = (-v + 1) / 2;
0953                if(n <= (int)boost::math::max_bernoulli_b2n<T>::value)
0954                   return T((-v & 1) ? -1 : 1) * boost::math::unchecked_bernoulli_b2n<T>(n) / (1 - v);
0955             }
0956             else if((v & 1) == 0)
0957             {
0958                if(((v / 2) <= (int)boost::math::max_bernoulli_b2n<T>::value) && (v <= (int)boost::math::max_factorial<T>::value))
0959                   return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * static_cast<T>(pow(constants::pi<T, Policy>(), T(v))) *
0960                      boost::math::unchecked_bernoulli_b2n<T>(v / 2) / boost::math::unchecked_factorial<T>(v);
0961                return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * static_cast<T>(pow(constants::pi<T, Policy>(), T(v))) *
0962                   boost::math::bernoulli_b2n<T>(v / 2) / boost::math::factorial<T>(v, pol);
0963             }
0964             else
0965                return zeta_imp_odd_integer(v, sc, pol, std::integral_constant<bool, (Tag::value <= 113) && Tag::value>());
0966          }
0967 #ifndef BOOST_NO_EXCEPTIONS
0968       }
0969       catch(const boost::math::rounding_error&){} // Just fall through, s is too large to round
0970       catch(const std::overflow_error&){}
0971 #endif
0972    }
0973
0974    if(fabs(s) < tools::root_epsilon<T>())
0975    {
0976       result = -0.5f - constants::log_root_two_pi<T, Policy>() * s;
0977    }
0978    else if(s < 0)
0979    {
0980       std::swap(s, sc);
0981       if(floor(sc/2) == sc/2)
0982          result = 0;
0983       else
0984       {
0985          if(s > max_factorial<T>::value)
0986          {
0987             T mult = boost::math::sin_pi(0.5f * sc, pol) * 2 * zeta_imp(s, sc, pol, tag);
0988             result = boost::math::lgamma(s, pol);
0989             result -= s * log(2 * constants::pi<T>());
0990             if(result > tools::log_max_value<T>())
0991                return sign(mult) * policies::raise_overflow_error<T>(function, nullptr, pol);
0992             result = exp(result);
0993             if(tools::max_value<T>() / fabs(mult) < result)
0994                return boost::math::sign(mult) * policies::raise_overflow_error<T>(function, nullptr, pol);
0995             result *= mult;
0996          }
0997          else
0998          {
0999             result = boost::math::sin_pi(0.5f * sc, pol)
1000                * 2 * pow(2 * constants::pi<T>(), -s)
1001                * boost::math::tgamma(s, pol)
1002                * zeta_imp(s, sc, pol, tag);
1003          }
1004       }
1005    }
1006    else
1007    {
1008       result = zeta_imp_prec(s, sc, pol, tag);
1009    }
1010    return result;
1011 }
1012
1013 template <class T, class Policy, class tag>
1014 struct zeta_initializer
1015 {
1016    struct init
1017    {
1018       init()
1019       {
1020          do_init(tag());
1021       }
1022       static void do_init(const std::integral_constant<int, 0>&){ boost::math::zeta(static_cast<T>(5), Policy()); }
1023       static void do_init(const std::integral_constant<int, 53>&){ boost::math::zeta(static_cast<T>(5), Policy()); }
1024       static void do_init(const std::integral_constant<int, 64>&)
1025       {
1026          boost::math::zeta(static_cast<T>(0.5), Policy());
1027          boost::math::zeta(static_cast<T>(1.5), Policy());
1028          boost::math::zeta(static_cast<T>(3.5), Policy());
1029          boost::math::zeta(static_cast<T>(6.5), Policy());
1030          boost::math::zeta(static_cast<T>(14.5), Policy());
1031          boost::math::zeta(static_cast<T>(40.5), Policy());
1032
1033          boost::math::zeta(static_cast<T>(5), Policy());
1034       }
1035       static void do_init(const std::integral_constant<int, 113>&)
1036       {
1037          boost::math::zeta(static_cast<T>(0.5), Policy());
1038          boost::math::zeta(static_cast<T>(1.5), Policy());
1039          boost::math::zeta(static_cast<T>(3.5), Policy());
1040          boost::math::zeta(static_cast<T>(5.5), Policy());
1041          boost::math::zeta(static_cast<T>(9.5), Policy());
1042          boost::math::zeta(static_cast<T>(16.5), Policy());
1043          boost::math::zeta(static_cast<T>(25.5), Policy());
1044          boost::math::zeta(static_cast<T>(70.5), Policy());
1045
1046          boost::math::zeta(static_cast<T>(5), Policy());
1047       }
1048       void force_instantiate()const{}
1049    };
1050    static const init initializer;
1051    static void force_instantiate()
1052    {
1053       initializer.force_instantiate();
1054    }
1055 };
1056
1057 template <class T, class Policy, class tag>
1058 const typename zeta_initializer<T, Policy, tag>::init zeta_initializer<T, Policy, tag>::initializer;
1059
1060 } // detail
1061
1062 template <class T, class Policy>
1063 inline typename tools::promote_args<T>::type zeta(T s, const Policy&)
1064 {
1065    typedef typename tools::promote_args<T>::type result_type;
1066    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1067    typedef typename policies::precision<result_type, Policy>::type precision_type;
1068    typedef typename policies::normalise<
1069       Policy,
1070       policies::promote_float<false>,
1071       policies::promote_double<false>,
1072       policies::discrete_quantile<>,
1073       policies::assert_undefined<> >::type forwarding_policy;
1074    typedef std::integral_constant<int,
1075       precision_type::value <= 0 ? 0 :
1076       precision_type::value <= 53 ? 53 :
1077       precision_type::value <= 64 ? 64 :
1078       precision_type::value <= 113 ? 113 : 0
1079    > tag_type;
1080
1081    detail::zeta_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
1082
1083    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(
1084       static_cast<value_type>(s),
1085       static_cast<value_type>(1 - static_cast<value_type>(s)),
1086       forwarding_policy(),
1087       tag_type()), "boost::math::zeta<%1%>(%1%)");
1088 }
1089
1090 template <class T>
1091 inline typename tools::promote_args<T>::type zeta(T s)
1092 {
1093    return zeta(s, policies::policy<>());
1094 }
1095
1096 }} // namespaces
1097
1098 #endif // BOOST_MATH_ZETA_HPP
1099
1100
1101