boost/math/special_functions/zeta.hpp

   1 //  Copyright John Maddock 2007.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_ZETA_HPP
   7 #define BOOST_MATH_ZETA_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/precision.hpp>
  14 #include <boost/math/tools/series.hpp>
  15 #include <boost/math/tools/big_constant.hpp>
  16 #include <boost/math/policies/error_handling.hpp>
  17 #include <boost/math/special_functions/gamma.hpp>
  18 #include <boost/math/special_functions/sin_pi.hpp>
  19
  20 namespace boost{ namespace math{ namespace detail{
  21
  22 #if 0
  23 //
  24 // This code is commented out because we have a better more rapidly converging series
  25 // now.  Retained for future reference and in case the new code causes any issues down the line....
  26 //
  27
  28 template <class T, class Policy>
  29 struct zeta_series_cache_size
  30 {
  31    //
  32    // Work how large to make our cache size when evaluating the series
  33    // evaluation:  normally this is just large enough for the series
  34    // to have converged, but for arbitrary precision types we need a
  35    // really large cache to achieve reasonable precision in a reasonable
  36    // time.  This is important when constructing rational approximations
  37    // to zeta for example.
  38    //
  39    typedef typename boost::math::policies::precision<T,Policy>::type precision_type;
  40    typedef typename mpl::if_<
  41       mpl::less_equal<precision_type, mpl::int_<0> >,
  42       mpl::int_<5000>,
  43       typename mpl::if_<
  44          mpl::less_equal<precision_type, mpl::int_<64> >,
  45          mpl::int_<70>,
  46          typename mpl::if_<
  47             mpl::less_equal<precision_type, mpl::int_<113> >,
  48             mpl::int_<100>,
  49             mpl::int_<5000>
  50          >::type
  51       >::type
  52    >::type type;
  53 };
  54
  55 template <class T, class Policy>
  56 T zeta_series_imp(T s, T sc, const Policy&)
  57 {
  58    //
  59    // Series evaluation from:
  60    // Havil, J. Gamma: Exploring Euler's Constant.
  61    // Princeton, NJ: Princeton University Press, 2003.
  62    //
  63    // See also http://mathworld.wolfram.com/RiemannZetaFunction.html
  64    //
  65    BOOST_MATH_STD_USING
  66    T sum = 0;
  67    T mult = 0.5;
  68    T change;
  69    typedef typename zeta_series_cache_size<T,Policy>::type cache_size;
  70    T powers[cache_size::value] = { 0, };
  71    unsigned n = 0;
  72    do{
  73       T binom = -static_cast<T>(n);
  74       T nested_sum = 1;
  75       if(n < sizeof(powers) / sizeof(powers[0]))
  76          powers[n] = pow(static_cast<T>(n + 1), -s);
  77       for(unsigned k = 1; k <= n; ++k)
  78       {
  79          T p;
  80          if(k < sizeof(powers) / sizeof(powers[0]))
  81          {
  82             p = powers[k];
  83             //p = pow(k + 1, -s);
  84          }
  85          else
  86             p = pow(static_cast<T>(k + 1), -s);
  87          nested_sum += binom * p;
  88         binom *= (k - static_cast<T>(n)) / (k + 1);
  89       }
  90       change = mult * nested_sum;
  91       sum += change;
  92       mult /= 2;
  93       ++n;
  94    }while(fabs(change / sum) > tools::epsilon<T>());
  95
  96    return sum * 1 / -boost::math::powm1(T(2), sc);
  97 }
  98
  99 //
 100 // Classical p-series:
 101 //
 102 template <class T>
 103 struct zeta_series2
 104 {
 105    typedef T result_type;
 106    zeta_series2(T _s) : s(-_s), k(1){}
 107    T operator()()
 108    {
 109       BOOST_MATH_STD_USING
 110       return pow(static_cast<T>(k++), s);
 111    }
 112 private:
 113    T s;
 114    unsigned k;
 115 };
 116
 117 template <class T, class Policy>
 118 inline T zeta_series2_imp(T s, const Policy& pol)
 119 {
 120    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();;
 121    zeta_series2<T> f(s);
 122    T result = tools::sum_series(
 123       f,
 124       policies::get_epsilon<T, Policy>(),
 125       max_iter);
 126    policies::check_series_iterations<T>("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
 127    return result;
 128 }
 129 #endif
 130
 131 template <class T, class Policy>
 132 T zeta_polynomial_series(T s, T sc, Policy const &)
 133 {
 134    //
 135    // This is algorithm 3 from:
 136    //
 137    // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein,
 138    // Canadian Mathematical Society, Conference Proceedings.
 139    // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
 140    //
 141    BOOST_MATH_STD_USING
 142    int n = itrunc(T(log(boost::math::tools::epsilon<T>()) / -2));
 143    T sum = 0;
 144    T two_n = ldexp(T(1), n);
 145    int ej_sign = 1;
 146    for(int j = 0; j < n; ++j)
 147    {
 148       sum += ej_sign * -two_n / pow(T(j + 1), s);
 149       ej_sign = -ej_sign;
 150    }
 151    T ej_sum = 1;
 152    T ej_term = 1;
 153    for(int j = n; j <= 2 * n - 1; ++j)
 154    {
 155       sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s);
 156       ej_sign = -ej_sign;
 157       ej_term *= 2 * n - j;
 158       ej_term /= j - n + 1;
 159       ej_sum += ej_term;
 160    }
 161    return -sum / (two_n * (-powm1(T(2), sc)));
 162 }
 163
 164 template <class T, class Policy>
 165 T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<0>&)
 166 {
 167    BOOST_MATH_STD_USING
 168    T result;
 169    result = zeta_polynomial_series(s, sc, pol);
 170 #if 0
 171    // Old code archived for future reference:
 172
 173    //
 174    // Only use power series if it will converge in 100
 175    // iterations or less: the more iterations it consumes
 176    // the slower convergence becomes so we have to be very
 177    // careful in it's usage.
 178    //
 179    if (s > -log(tools::epsilon<T>()) / 4.5)
 180       result = detail::zeta_series2_imp(s, pol);
 181    else
 182       result = detail::zeta_series_imp(s, sc, pol);
 183 #endif
 184    return result;
 185 }
 186
 187 template <class T, class Policy>
 188 inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
 189 {
 190    BOOST_MATH_STD_USING
 191    T result;
 192    if(s < 1)
 193    {
 194       // Rational Approximation
 195       // Maximum Deviation Found:                     2.020e-18
 196       // Expected Error Term:                         -2.020e-18
 197       // Max error found at double precision:         3.994987e-17
 198       static const T P[6] = {
 199          0.24339294433593750202L,
 200          -0.49092470516353571651L,
 201          0.0557616214776046784287L,
 202          -0.00320912498879085894856L,
 203          0.000451534528645796438704L,
 204          -0.933241270357061460782e-5L,
 205         };
 206       static const T Q[6] = {
 207          1L,
 208          -0.279960334310344432495L,
 209          0.0419676223309986037706L,
 210          -0.00413421406552171059003L,
 211          0.00024978985622317935355L,
 212          -0.101855788418564031874e-4L,
 213       };
 214       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
 215       result -= 1.2433929443359375F;
 216       result += (sc);
 217       result /= (sc);
 218    }
 219    else if(s <= 2)
 220    {
 221       // Maximum Deviation Found:        9.007e-20
 222       // Expected Error Term:            9.007e-20
 223       static const T P[6] = {
 224          0.577215664901532860516,
 225          0.243210646940107164097,
 226          0.0417364673988216497593,
 227          0.00390252087072843288378,
 228          0.000249606367151877175456,
 229          0.110108440976732897969e-4,
 230       };
 231       static const T Q[6] = {
 232          1,
 233          0.295201277126631761737,
 234          0.043460910607305495864,
 235          0.00434930582085826330659,
 236          0.000255784226140488490982,
 237          0.10991819782396112081e-4,
 238       };
 239       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
 240       result += 1 / (-sc);
 241    }
 242    else if(s <= 4)
 243    {
 244       // Maximum Deviation Found:          5.946e-22
 245       // Expected Error Term:              -5.946e-22
 246       static const float Y = 0.6986598968505859375;
 247       static const T P[6] = {
 248          -0.0537258300023595030676,
 249          0.0445163473292365591906,
 250          0.0128677673534519952905,
 251          0.00097541770457391752726,
 252          0.769875101573654070925e-4,
 253          0.328032510000383084155e-5,
 254       };
 255       static const T Q[7] = {
 256          1,
 257          0.33383194553034051422,
 258          0.0487798431291407621462,
 259          0.00479039708573558490716,
 260          0.000270776703956336357707,
 261          0.106951867532057341359e-4,
 262          0.236276623974978646399e-7,
 263       };
 264       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
 265       result += Y + 1 / (-sc);
 266    }
 267    else if(s <= 7)
 268    {
 269       // Maximum Deviation Found:                     2.955e-17
 270       // Expected Error Term:                         2.955e-17
 271       // Max error found at double precision:         2.009135e-16
 272
 273       static const T P[6] = {
 274          -2.49710190602259410021,
 275          -2.60013301809475665334,
 276          -0.939260435377109939261,
 277          -0.138448617995741530935,
 278          -0.00701721240549802377623,
 279          -0.229257310594893932383e-4,
 280       };
 281       static const T Q[9] = {
 282          1,
 283          0.706039025937745133628,
 284          0.15739599649558626358,
 285          0.0106117950976845084417,
 286          -0.36910273311764618902e-4,
 287          0.493409563927590008943e-5,
 288          -0.234055487025287216506e-6,
 289          0.718833729365459760664e-8,
 290          -0.1129200113474947419e-9,
 291       };
 292       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
 293       result = 1 + exp(result);
 294    }
 295    else if(s < 15)
 296    {
 297       // Maximum Deviation Found:                     7.117e-16
 298       // Expected Error Term:                         7.117e-16
 299       // Max error found at double precision:         9.387771e-16
 300       static const T P[7] = {
 301          -4.78558028495135619286,
 302          -1.89197364881972536382,
 303          -0.211407134874412820099,
 304          -0.000189204758260076688518,
 305          0.00115140923889178742086,
 306          0.639949204213164496988e-4,
 307          0.139348932445324888343e-5,
 308         };
 309       static const T Q[9] = {
 310          1,
 311          0.244345337378188557777,
 312          0.00873370754492288653669,
 313          -0.00117592765334434471562,
 314          -0.743743682899933180415e-4,
 315          -0.21750464515767984778e-5,
 316          0.471001264003076486547e-8,
 317          -0.833378440625385520576e-10,
 318          0.699841545204845636531e-12,
 319         };
 320       result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));
 321       result = 1 + exp(result);
 322    }
 323    else if(s < 36)
 324    {
 325       // Max error in interpolated form:             1.668e-17
 326       // Max error found at long double precision:   1.669714e-17
 327       static const T P[8] = {
 328          -10.3948950573308896825,
 329          -2.85827219671106697179,
 330          -0.347728266539245787271,
 331          -0.0251156064655346341766,
 332          -0.00119459173416968685689,
 333          -0.382529323507967522614e-4,
 334          -0.785523633796723466968e-6,
 335          -0.821465709095465524192e-8,
 336       };
 337       static const T Q[10] = {
 338          1,
 339          0.208196333572671890965,
 340          0.0195687657317205033485,
 341          0.00111079638102485921877,
 342          0.408507746266039256231e-4,
 343          0.955561123065693483991e-6,
 344          0.118507153474022900583e-7,
 345          0.222609483627352615142e-14,
 346       };
 347       result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));
 348       result = 1 + exp(result);
 349    }
 350    else if(s < 56)
 351    {
 352       result = 1 + pow(T(2), -s);
 353    }
 354    else
 355    {
 356       result = 1;
 357    }
 358    return result;
 359 }
 360
 361 template <class T, class Policy>
 362 T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<64>&)
 363 {
 364    BOOST_MATH_STD_USING
 365    T result;
 366    if(s < 1)
 367    {
 368       // Rational Approximation
 369       // Maximum Deviation Found:                     3.099e-20
 370       // Expected Error Term:                         3.099e-20
 371       // Max error found at long double precision:    5.890498e-20
 372       static const T P[6] = {
 373          BOOST_MATH_BIG_CONSTANT(T, 64, 0.243392944335937499969),
 374          BOOST_MATH_BIG_CONSTANT(T, 64, -0.496837806864865688082),
 375          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0680008039723709987107),
 376          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00511620413006619942112),
 377          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000455369899250053003335),
 378          BOOST_MATH_BIG_CONSTANT(T, 64, -0.279496685273033761927e-4),
 379         };
 380       static const T Q[7] = {
 381          BOOST_MATH_BIG_CONSTANT(T, 64, 1),
 382          BOOST_MATH_BIG_CONSTANT(T, 64, -0.30425480068225790522),
 383          BOOST_MATH_BIG_CONSTANT(T, 64, 0.050052748580371598736),
 384          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00519355671064700627862),
 385          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000360623385771198350257),
 386          BOOST_MATH_BIG_CONSTANT(T, 64, -0.159600883054550987633e-4),
 387          BOOST_MATH_BIG_CONSTANT(T, 64, 0.339770279812410586032e-6),
 388       };
 389       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
 390       result -= 1.2433929443359375F;
 391       result += (sc);
 392       result /= (sc);
 393    }
 394    else if(s <= 2)
 395    {
 396       // Maximum Deviation Found:                     1.059e-21
 397       // Expected Error Term:                         1.059e-21
 398       // Max error found at long double precision:    1.626303e-19
 399
 400       static const T P[6] = {
 401          BOOST_MATH_BIG_CONSTANT(T, 64, 0.577215664901532860605),
 402          BOOST_MATH_BIG_CONSTANT(T, 64, 0.222537368917162139445),
 403          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0356286324033215682729),
 404          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00304465292366350081446),
 405          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000178102511649069421904),
 406          BOOST_MATH_BIG_CONSTANT(T, 64, 0.700867470265983665042e-5),
 407       };
 408       static const T Q[7] = {
 409          BOOST_MATH_BIG_CONSTANT(T, 64, 1),
 410          BOOST_MATH_BIG_CONSTANT(T, 64, 0.259385759149531030085),
 411          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0373974962106091316854),
 412          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00332735159183332820617),
 413          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000188690420706998606469),
 414          BOOST_MATH_BIG_CONSTANT(T, 64, 0.635994377921861930071e-5),
 415          BOOST_MATH_BIG_CONSTANT(T, 64, 0.226583954978371199405e-7),
 416       };
 417       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
 418       result += 1 / (-sc);
 419    }
 420    else if(s <= 4)
 421    {
 422       // Maximum Deviation Found:          5.946e-22
 423       // Expected Error Term:              -5.946e-22
 424       static const float Y = 0.6986598968505859375;
 425       static const T P[7] = {
 426          BOOST_MATH_BIG_CONSTANT(T, 64, -0.053725830002359501027),
 427          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0470551187571475844778),
 428          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0101339410415759517471),
 429          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00100240326666092854528),
 430          BOOST_MATH_BIG_CONSTANT(T, 64, 0.685027119098122814867e-4),
 431          BOOST_MATH_BIG_CONSTANT(T, 64, 0.390972820219765942117e-5),
 432          BOOST_MATH_BIG_CONSTANT(T, 64, 0.540319769113543934483e-7),
 433       };
 434       static const T Q[8] = {
 435          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 436          BOOST_MATH_BIG_CONSTANT(T, 64, 0.286577739726542730421),
 437          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0447355811517733225843),
 438          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00430125107610252363302),
 439          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000284956969089786662045),
 440          BOOST_MATH_BIG_CONSTANT(T, 64, 0.116188101609848411329e-4),
 441          BOOST_MATH_BIG_CONSTANT(T, 64, 0.278090318191657278204e-6),
 442          BOOST_MATH_BIG_CONSTANT(T, 64, -0.19683620233222028478e-8),
 443       };
 444       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
 445       result += Y + 1 / (-sc);
 446    }
 447    else if(s <= 7)
 448    {
 449       // Max error found at long double precision: 8.132216e-19
 450       static const T P[8] = {
 451          BOOST_MATH_BIG_CONSTANT(T, 64, -2.49710190602259407065),
 452          BOOST_MATH_BIG_CONSTANT(T, 64, -3.36664913245960625334),
 453          BOOST_MATH_BIG_CONSTANT(T, 64, -1.77180020623777595452),
 454          BOOST_MATH_BIG_CONSTANT(T, 64, -0.464717885249654313933),
 455          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0643694921293579472583),
 456          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00464265386202805715487),
 457          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000165556579779704340166),
 458          BOOST_MATH_BIG_CONSTANT(T, 64, -0.252884970740994069582e-5),
 459       };
 460       static const T Q[9] = {
 461          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 462          BOOST_MATH_BIG_CONSTANT(T, 64, 1.01300131390690459085),
 463          BOOST_MATH_BIG_CONSTANT(T, 64, 0.387898115758643503827),
 464          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0695071490045701135188),
 465          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00586908595251442839291),
 466          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000217752974064612188616),
 467          BOOST_MATH_BIG_CONSTANT(T, 64, 0.397626583349419011731e-5),
 468          BOOST_MATH_BIG_CONSTANT(T, 64, -0.927884739284359700764e-8),
 469          BOOST_MATH_BIG_CONSTANT(T, 64, 0.119810501805618894381e-9),
 470       };
 471       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
 472       result = 1 + exp(result);
 473    }
 474    else if(s < 15)
 475    {
 476       // Max error in interpolated form:              1.133e-18
 477       // Max error found at long double precision:    2.183198e-18
 478       static const T P[9] = {
 479          BOOST_MATH_BIG_CONSTANT(T, 64, -4.78558028495135548083),
 480          BOOST_MATH_BIG_CONSTANT(T, 64, -3.23873322238609358947),
 481          BOOST_MATH_BIG_CONSTANT(T, 64, -0.892338582881021799922),
 482          BOOST_MATH_BIG_CONSTANT(T, 64, -0.131326296217965913809),
 483          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0115651591773783712996),
 484          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000657728968362695775205),
 485          BOOST_MATH_BIG_CONSTANT(T, 64, -0.252051328129449973047e-4),
 486          BOOST_MATH_BIG_CONSTANT(T, 64, -0.626503445372641798925e-6),
 487          BOOST_MATH_BIG_CONSTANT(T, 64, -0.815696314790853893484e-8),
 488         };
 489       static const T Q[9] = {
 490          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 491          BOOST_MATH_BIG_CONSTANT(T, 64, 0.525765665400123515036),
 492          BOOST_MATH_BIG_CONSTANT(T, 64, 0.10852641753657122787),
 493          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0115669945375362045249),
 494          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000732896513858274091966),
 495          BOOST_MATH_BIG_CONSTANT(T, 64, 0.30683952282420248448e-4),
 496          BOOST_MATH_BIG_CONSTANT(T, 64, 0.819649214609633126119e-6),
 497          BOOST_MATH_BIG_CONSTANT(T, 64, 0.117957556472335968146e-7),
 498          BOOST_MATH_BIG_CONSTANT(T, 64, -0.193432300973017671137e-12),
 499         };
 500       result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));
 501       result = 1 + exp(result);
 502    }
 503    else if(s < 42)
 504    {
 505       // Max error in interpolated form:             1.668e-17
 506       // Max error found at long double precision:   1.669714e-17
 507       static const T P[9] = {
 508          BOOST_MATH_BIG_CONSTANT(T, 64, -10.3948950573308861781),
 509          BOOST_MATH_BIG_CONSTANT(T, 64, -2.82646012777913950108),
 510          BOOST_MATH_BIG_CONSTANT(T, 64, -0.342144362739570333665),
 511          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0249285145498722647472),
 512          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00122493108848097114118),
 513          BOOST_MATH_BIG_CONSTANT(T, 64, -0.423055371192592850196e-4),
 514          BOOST_MATH_BIG_CONSTANT(T, 64, -0.1025215577185967488e-5),
 515          BOOST_MATH_BIG_CONSTANT(T, 64, -0.165096762663509467061e-7),
 516          BOOST_MATH_BIG_CONSTANT(T, 64, -0.145392555873022044329e-9),
 517       };
 518       static const T Q[10] = {
 519          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 520          BOOST_MATH_BIG_CONSTANT(T, 64, 0.205135978585281988052),
 521          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0192359357875879453602),
 522          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00111496452029715514119),
 523          BOOST_MATH_BIG_CONSTANT(T, 64, 0.434928449016693986857e-4),
 524          BOOST_MATH_BIG_CONSTANT(T, 64, 0.116911068726610725891e-5),
 525          BOOST_MATH_BIG_CONSTANT(T, 64, 0.206704342290235237475e-7),
 526          BOOST_MATH_BIG_CONSTANT(T, 64, 0.209772836100827647474e-9),
 527          BOOST_MATH_BIG_CONSTANT(T, 64, -0.939798249922234703384e-16),
 528          BOOST_MATH_BIG_CONSTANT(T, 64, 0.264584017421245080294e-18),
 529       };
 530       result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));
 531       result = 1 + exp(result);
 532    }
 533    else if(s < 63)
 534    {
 535       result = 1 + pow(T(2), -s);
 536    }
 537    else
 538    {
 539       result = 1;
 540    }
 541    return result;
 542 }
 543
 544 template <class T, class Policy>
 545 T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<113>&)
 546 {
 547    BOOST_MATH_STD_USING
 548    T result;
 549    if(s < 1)
 550    {
 551       // Rational Approximation
 552       // Maximum Deviation Found:                     9.493e-37
 553       // Expected Error Term:                         9.492e-37
 554       // Max error found at long double precision:    7.281332e-31
 555
 556       static const T P[10] = {
 557          BOOST_MATH_BIG_CONSTANT(T, 113, -1),
 558          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0353008629988648122808504280990313668),
 559          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0107795651204927743049369868548706909),
 560          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000523961870530500751114866884685172975),
 561          BOOST_MATH_BIG_CONSTANT(T, 113, -0.661805838304910731947595897966487515e-4),
 562          BOOST_MATH_BIG_CONSTANT(T, 113, -0.658932670403818558510656304189164638e-5),
 563          BOOST_MATH_BIG_CONSTANT(T, 113, -0.103437265642266106533814021041010453e-6),
 564          BOOST_MATH_BIG_CONSTANT(T, 113, 0.116818787212666457105375746642927737e-7),
 565          BOOST_MATH_BIG_CONSTANT(T, 113, 0.660690993901506912123512551294239036e-9),
 566          BOOST_MATH_BIG_CONSTANT(T, 113, 0.113103113698388531428914333768142527e-10),
 567         };
 568       static const T Q[11] = {
 569          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 570          BOOST_MATH_BIG_CONSTANT(T, 113, -0.387483472099602327112637481818565459),
 571          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0802265315091063135271497708694776875),
 572          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0110727276164171919280036408995078164),
 573          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00112552716946286252000434849173787243),
 574          BOOST_MATH_BIG_CONSTANT(T, 113, -0.874554160748626916455655180296834352e-4),
 575          BOOST_MATH_BIG_CONSTANT(T, 113, 0.530097847491828379568636739662278322e-5),
 576          BOOST_MATH_BIG_CONSTANT(T, 113, -0.248461553590496154705565904497247452e-6),
 577          BOOST_MATH_BIG_CONSTANT(T, 113, 0.881834921354014787309644951507523899e-8),
 578          BOOST_MATH_BIG_CONSTANT(T, 113, -0.217062446168217797598596496310953025e-9),
 579          BOOST_MATH_BIG_CONSTANT(T, 113, 0.315823200002384492377987848307151168e-11),
 580       };
 581       result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
 582       result += (sc);
 583       result /= (sc);
 584    }
 585    else if(s <= 2)
 586    {
 587       // Maximum Deviation Found:                     1.616e-37
 588       // Expected Error Term:                         -1.615e-37
 589
 590       static const T P[10] = {
 591          BOOST_MATH_BIG_CONSTANT(T, 113, 0.577215664901532860606512090082402431),
 592          BOOST_MATH_BIG_CONSTANT(T, 113, 0.255597968739771510415479842335906308),
 593          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0494056503552807274142218876983542205),
 594          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00551372778611700965268920983472292325),
 595          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00043667616723970574871427830895192731),
 596          BOOST_MATH_BIG_CONSTANT(T, 113, 0.268562259154821957743669387915239528e-4),
 597          BOOST_MATH_BIG_CONSTANT(T, 113, 0.109249633923016310141743084480436612e-5),
 598          BOOST_MATH_BIG_CONSTANT(T, 113, 0.273895554345300227466534378753023924e-7),
 599          BOOST_MATH_BIG_CONSTANT(T, 113, 0.583103205551702720149237384027795038e-9),
 600          BOOST_MATH_BIG_CONSTANT(T, 113, -0.835774625259919268768735944711219256e-11),
 601       };
 602       static const T Q[11] = {
 603          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 604          BOOST_MATH_BIG_CONSTANT(T, 113, 0.316661751179735502065583176348292881),
 605          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0540401806533507064453851182728635272),
 606          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00598621274107420237785899476374043797),
 607          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000474907812321704156213038740142079615),
 608          BOOST_MATH_BIG_CONSTANT(T, 113, 0.272125421722314389581695715835862418e-4),
 609          BOOST_MATH_BIG_CONSTANT(T, 113, 0.112649552156479800925522445229212933e-5),
 610          BOOST_MATH_BIG_CONSTANT(T, 113, 0.301838975502992622733000078063330461e-7),
 611          BOOST_MATH_BIG_CONSTANT(T, 113, 0.422960728687211282539769943184270106e-9),
 612          BOOST_MATH_BIG_CONSTANT(T, 113, -0.377105263588822468076813329270698909e-11),
 613          BOOST_MATH_BIG_CONSTANT(T, 113, -0.581926559304525152432462127383600681e-13),
 614       };
 615       result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));
 616       result += 1 / (-sc);
 617    }
 618    else if(s <= 4)
 619    {
 620       // Maximum Deviation Found:                     1.891e-36
 621       // Expected Error Term:                         -1.891e-36
 622       // Max error found: 2.171527e-35
 623
 624       static const float Y = 0.6986598968505859375;
 625       static const T P[11] = {
 626          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0537258300023595010275848333539748089),
 627          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0429086930802630159457448174466342553),
 628          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0136148228754303412510213395034056857),
 629          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00190231601036042925183751238033763915),
 630          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000186880390916311438818302549192456581),
 631          BOOST_MATH_BIG_CONSTANT(T, 113, 0.145347370745893262394287982691323657e-4),
 632          BOOST_MATH_BIG_CONSTANT(T, 113, 0.805843276446813106414036600485884885e-6),
 633          BOOST_MATH_BIG_CONSTANT(T, 113, 0.340818159286739137503297172091882574e-7),
 634          BOOST_MATH_BIG_CONSTANT(T, 113, 0.115762357488748996526167305116837246e-8),
 635          BOOST_MATH_BIG_CONSTANT(T, 113, 0.231904754577648077579913403645767214e-10),
 636          BOOST_MATH_BIG_CONSTANT(T, 113, 0.340169592866058506675897646629036044e-12),
 637       };
 638       static const T Q[12] = {
 639          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 640          BOOST_MATH_BIG_CONSTANT(T, 113, 0.363755247765087100018556983050520554),
 641          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0696581979014242539385695131258321598),
 642          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00882208914484611029571547753782014817),
 643          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000815405623261946661762236085660996718),
 644          BOOST_MATH_BIG_CONSTANT(T, 113, 0.571366167062457197282642344940445452e-4),
 645          BOOST_MATH_BIG_CONSTANT(T, 113, 0.309278269271853502353954062051797838e-5),
 646          BOOST_MATH_BIG_CONSTANT(T, 113, 0.12822982083479010834070516053794262e-6),
 647          BOOST_MATH_BIG_CONSTANT(T, 113, 0.397876357325018976733953479182110033e-8),
 648          BOOST_MATH_BIG_CONSTANT(T, 113, 0.8484432107648683277598472295289279e-10),
 649          BOOST_MATH_BIG_CONSTANT(T, 113, 0.105677416606909614301995218444080615e-11),
 650          BOOST_MATH_BIG_CONSTANT(T, 113, 0.547223964564003701979951154093005354e-15),
 651       };
 652       result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));
 653       result += Y + 1 / (-sc);
 654    }
 655    else if(s <= 6)
 656    {
 657       // Max error in interpolated form:             1.510e-37
 658       // Max error found at long double precision:   2.769266e-34
 659
 660       static const T Y = 3.28348541259765625F;
 661
 662       static const T P[13] = {
 663          BOOST_MATH_BIG_CONSTANT(T, 113, 0.786383506575062179339611614117697622),
 664          BOOST_MATH_BIG_CONSTANT(T, 113, 0.495766593395271370974685959652073976),
 665          BOOST_MATH_BIG_CONSTANT(T, 113, -0.409116737851754766422360889037532228),
 666          BOOST_MATH_BIG_CONSTANT(T, 113, -0.57340744006238263817895456842655987),
 667          BOOST_MATH_BIG_CONSTANT(T, 113, -0.280479899797421910694892949057963111),
 668          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0753148409447590257157585696212649869),
 669          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0122934003684672788499099362823748632),
 670          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126148398446193639247961370266962927),
 671          BOOST_MATH_BIG_CONSTANT(T, 113, -0.828465038179772939844657040917364896e-4),
 672          BOOST_MATH_BIG_CONSTANT(T, 113, -0.361008916706050977143208468690645684e-5),
 673          BOOST_MATH_BIG_CONSTANT(T, 113, -0.109879825497910544424797771195928112e-6),
 674          BOOST_MATH_BIG_CONSTANT(T, 113, -0.214539416789686920918063075528797059e-8),
 675          BOOST_MATH_BIG_CONSTANT(T, 113, -0.15090220092460596872172844424267351e-10),
 676       };
 677       static const T Q[14] = {
 678          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 679          BOOST_MATH_BIG_CONSTANT(T, 113, 1.69490865837142338462982225731926485),
 680          BOOST_MATH_BIG_CONSTANT(T, 113, 1.22697696630994080733321401255942464),
 681          BOOST_MATH_BIG_CONSTANT(T, 113, 0.495409420862526540074366618006341533),
 682          BOOST_MATH_BIG_CONSTANT(T, 113, 0.122368084916843823462872905024259633),
 683          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0191412993625268971656513890888208623),
 684          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00191401538628980617753082598351559642),
 685          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000123318142456272424148930280876444459),
 686          BOOST_MATH_BIG_CONSTANT(T, 113, 0.531945488232526067889835342277595709e-5),
 687          BOOST_MATH_BIG_CONSTANT(T, 113, 0.161843184071894368337068779669116236e-6),
 688          BOOST_MATH_BIG_CONSTANT(T, 113, 0.305796079600152506743828859577462778e-8),
 689          BOOST_MATH_BIG_CONSTANT(T, 113, 0.233582592298450202680170811044408894e-10),
 690          BOOST_MATH_BIG_CONSTANT(T, 113, -0.275363878344548055574209713637734269e-13),
 691          BOOST_MATH_BIG_CONSTANT(T, 113, 0.221564186807357535475441900517843892e-15),
 692       };
 693       result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));
 694       result -= Y;
 695       result = 1 + exp(result);
 696    }
 697    else if(s < 10)
 698    {
 699       // Max error in interpolated form:             1.999e-34
 700       // Max error found at long double precision:   2.156186e-33
 701
 702       static const T P[13] = {
 703          BOOST_MATH_BIG_CONSTANT(T, 113, -4.0545627381873738086704293881227365),
 704          BOOST_MATH_BIG_CONSTANT(T, 113, -4.70088348734699134347906176097717782),
 705          BOOST_MATH_BIG_CONSTANT(T, 113, -2.36921550900925512951976617607678789),
 706          BOOST_MATH_BIG_CONSTANT(T, 113, -0.684322583796369508367726293719322866),
 707          BOOST_MATH_BIG_CONSTANT(T, 113, -0.126026534540165129870721937592996324),
 708          BOOST_MATH_BIG_CONSTANT(T, 113, -0.015636903921778316147260572008619549),
 709          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00135442294754728549644376325814460807),
 710          BOOST_MATH_BIG_CONSTANT(T, 113, -0.842793965853572134365031384646117061e-4),
 711          BOOST_MATH_BIG_CONSTANT(T, 113, -0.385602133791111663372015460784978351e-5),
 712          BOOST_MATH_BIG_CONSTANT(T, 113, -0.130458500394692067189883214401478539e-6),
 713          BOOST_MATH_BIG_CONSTANT(T, 113, -0.315861074947230418778143153383660035e-8),
 714          BOOST_MATH_BIG_CONSTANT(T, 113, -0.500334720512030826996373077844707164e-10),
 715          BOOST_MATH_BIG_CONSTANT(T, 113, -0.420204769185233365849253969097184005e-12),
 716         };
 717       static const T Q[14] = {
 718          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 719          BOOST_MATH_BIG_CONSTANT(T, 113, 0.97663511666410096104783358493318814),
 720          BOOST_MATH_BIG_CONSTANT(T, 113, 0.40878780231201806504987368939673249),
 721          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0963890666609396058945084107597727252),
 722          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0142207619090854604824116070866614505),
 723          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139010220902667918476773423995750877),
 724          BOOST_MATH_BIG_CONSTANT(T, 113, 0.940669540194694997889636696089994734e-4),
 725          BOOST_MATH_BIG_CONSTANT(T, 113, 0.458220848507517004399292480807026602e-5),
 726          BOOST_MATH_BIG_CONSTANT(T, 113, 0.16345521617741789012782420625435495e-6),
 727          BOOST_MATH_BIG_CONSTANT(T, 113, 0.414007452533083304371566316901024114e-8),
 728          BOOST_MATH_BIG_CONSTANT(T, 113, 0.68701473543366328016953742622661377e-10),
 729          BOOST_MATH_BIG_CONSTANT(T, 113, 0.603461891080716585087883971886075863e-12),
 730          BOOST_MATH_BIG_CONSTANT(T, 113, 0.294670713571839023181857795866134957e-16),
 731          BOOST_MATH_BIG_CONSTANT(T, 113, -0.147003914536437243143096875069813451e-18),
 732         };
 733       result = tools::evaluate_polynomial(P, T(s - 6)) / tools::evaluate_polynomial(Q, T(s - 6));
 734       result = 1 + exp(result);
 735    }
 736    else if(s < 17)
 737    {
 738       // Max error in interpolated form:             1.641e-32
 739       // Max error found at long double precision:   1.696121e-32
 740       static const T P[13] = {
 741          BOOST_MATH_BIG_CONSTANT(T, 113, -6.91319491921722925920883787894829678),
 742          BOOST_MATH_BIG_CONSTANT(T, 113, -3.65491257639481960248690596951049048),
 743          BOOST_MATH_BIG_CONSTANT(T, 113, -0.813557553449954526442644544105257881),
 744          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0994317301685870959473658713841138083),
 745          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00726896610245676520248617014211734906),
 746          BOOST_MATH_BIG_CONSTANT(T, 113, -0.000317253318715075854811266230916762929),
 747          BOOST_MATH_BIG_CONSTANT(T, 113, -0.66851422826636750855184211580127133e-5),
 748          BOOST_MATH_BIG_CONSTANT(T, 113, 0.879464154730985406003332577806849971e-7),
 749          BOOST_MATH_BIG_CONSTANT(T, 113, 0.113838903158254250631678791998294628e-7),
 750          BOOST_MATH_BIG_CONSTANT(T, 113, 0.379184410304927316385211327537817583e-9),
 751          BOOST_MATH_BIG_CONSTANT(T, 113, 0.612992858643904887150527613446403867e-11),
 752          BOOST_MATH_BIG_CONSTANT(T, 113, 0.347873737198164757035457841688594788e-13),
 753          BOOST_MATH_BIG_CONSTANT(T, 113, -0.289187187441625868404494665572279364e-15),
 754         };
 755       static const T Q[14] = {
 756          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 757          BOOST_MATH_BIG_CONSTANT(T, 113, 0.427310044448071818775721584949868806),
 758          BOOST_MATH_BIG_CONSTANT(T, 113, 0.074602514873055756201435421385243062),
 759          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00688651562174480772901425121653945942),
 760          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000360174847635115036351323894321880445),
 761          BOOST_MATH_BIG_CONSTANT(T, 113, 0.973556847713307543918865405758248777e-5),
 762          BOOST_MATH_BIG_CONSTANT(T, 113, -0.853455848314516117964634714780874197e-8),
 763          BOOST_MATH_BIG_CONSTANT(T, 113, -0.118203513654855112421673192194622826e-7),
 764          BOOST_MATH_BIG_CONSTANT(T, 113, -0.462521662511754117095006543363328159e-9),
 765          BOOST_MATH_BIG_CONSTANT(T, 113, -0.834212591919475633107355719369463143e-11),
 766          BOOST_MATH_BIG_CONSTANT(T, 113, -0.5354594751002702935740220218582929e-13),
 767          BOOST_MATH_BIG_CONSTANT(T, 113, 0.406451690742991192964889603000756203e-15),
 768          BOOST_MATH_BIG_CONSTANT(T, 113, 0.887948682401000153828241615760146728e-19),
 769          BOOST_MATH_BIG_CONSTANT(T, 113, -0.34980761098820347103967203948619072e-21),
 770         };
 771       result = tools::evaluate_polynomial(P, T(s - 10)) / tools::evaluate_polynomial(Q, T(s - 10));
 772       result = 1 + exp(result);
 773    }
 774    else if(s < 30)
 775    {
 776       // Max error in interpolated form:             1.563e-31
 777       // Max error found at long double precision:   1.562725e-31
 778
 779       static const T P[13] = {
 780          BOOST_MATH_BIG_CONSTANT(T, 113, -11.7824798233959252791987402769438322),
 781          BOOST_MATH_BIG_CONSTANT(T, 113, -4.36131215284987731928174218354118102),
 782          BOOST_MATH_BIG_CONSTANT(T, 113, -0.732260980060982349410898496846972204),
 783          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0744985185694913074484248803015717388),
 784          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00517228281320594683022294996292250527),
 785          BOOST_MATH_BIG_CONSTANT(T, 113, -0.000260897206152101522569969046299309939),
 786          BOOST_MATH_BIG_CONSTANT(T, 113, -0.989553462123121764865178453128769948e-5),
 787          BOOST_MATH_BIG_CONSTANT(T, 113, -0.286916799741891410827712096608826167e-6),
 788          BOOST_MATH_BIG_CONSTANT(T, 113, -0.637262477796046963617949532211619729e-8),
 789          BOOST_MATH_BIG_CONSTANT(T, 113, -0.106796831465628373325491288787760494e-9),
 790          BOOST_MATH_BIG_CONSTANT(T, 113, -0.129343095511091870860498356205376823e-11),
 791          BOOST_MATH_BIG_CONSTANT(T, 113, -0.102397936697965977221267881716672084e-13),
 792          BOOST_MATH_BIG_CONSTANT(T, 113, -0.402663128248642002351627980255756363e-16),
 793       };
 794       static const T Q[14] = {
 795          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 796          BOOST_MATH_BIG_CONSTANT(T, 113, 0.311288325355705609096155335186466508),
 797          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0438318468940415543546769437752132748),
 798          BOOST_MATH_BIG_CONSTANT(T, 113, 0.00374396349183199548610264222242269536),
 799          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000218707451200585197339671707189281302),
 800          BOOST_MATH_BIG_CONSTANT(T, 113, 0.927578767487930747532953583797351219e-5),
 801          BOOST_MATH_BIG_CONSTANT(T, 113, 0.294145760625753561951137473484889639e-6),
 802          BOOST_MATH_BIG_CONSTANT(T, 113, 0.704618586690874460082739479535985395e-8),
 803          BOOST_MATH_BIG_CONSTANT(T, 113, 0.126333332872897336219649130062221257e-9),
 804          BOOST_MATH_BIG_CONSTANT(T, 113, 0.16317315713773503718315435769352765e-11),
 805          BOOST_MATH_BIG_CONSTANT(T, 113, 0.137846712823719515148344938160275695e-13),
 806          BOOST_MATH_BIG_CONSTANT(T, 113, 0.580975420554224366450994232723910583e-16),
 807          BOOST_MATH_BIG_CONSTANT(T, 113, -0.291354445847552426900293580511392459e-22),
 808          BOOST_MATH_BIG_CONSTANT(T, 113, 0.73614324724785855925025452085443636e-25),
 809       };
 810       result = tools::evaluate_polynomial(P, T(s - 17)) / tools::evaluate_polynomial(Q, T(s - 17));
 811       result = 1 + exp(result);
 812    }
 813    else if(s < 74)
 814    {
 815       // Max error in interpolated form:             2.311e-27
 816       // Max error found at long double precision:   2.297544e-27
 817       static const T P[14] = {
 818          BOOST_MATH_BIG_CONSTANT(T, 113, -20.7944102007844314586649688802236072),
 819          BOOST_MATH_BIG_CONSTANT(T, 113, -4.95759941987499442499908748130192187),
 820          BOOST_MATH_BIG_CONSTANT(T, 113, -0.563290752832461751889194629200298688),
 821          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0406197001137935911912457120706122877),
 822          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0020846534789473022216888863613422293),
 823          BOOST_MATH_BIG_CONSTANT(T, 113, -0.808095978462109173749395599401375667e-4),
 824          BOOST_MATH_BIG_CONSTANT(T, 113, -0.244706022206249301640890603610060959e-5),
 825          BOOST_MATH_BIG_CONSTANT(T, 113, -0.589477682919645930544382616501666572e-7),
 826          BOOST_MATH_BIG_CONSTANT(T, 113, -0.113699573675553496343617442433027672e-8),
 827          BOOST_MATH_BIG_CONSTANT(T, 113, -0.174767860183598149649901223128011828e-10),
 828          BOOST_MATH_BIG_CONSTANT(T, 113, -0.210051620306761367764549971980026474e-12),
 829          BOOST_MATH_BIG_CONSTANT(T, 113, -0.189187969537370950337212675466400599e-14),
 830          BOOST_MATH_BIG_CONSTANT(T, 113, -0.116313253429564048145641663778121898e-16),
 831          BOOST_MATH_BIG_CONSTANT(T, 113, -0.376708747782400769427057630528578187e-19),
 832       };
 833       static const T Q[16] = {
 834          BOOST_MATH_BIG_CONSTANT(T, 113, 1),
 835          BOOST_MATH_BIG_CONSTANT(T, 113, 0.205076752981410805177554569784219717),
 836          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0202526722696670378999575738524540269),
 837          BOOST_MATH_BIG_CONSTANT(T, 113, 0.001278305290005994980069466658219057),
 838          BOOST_MATH_BIG_CONSTANT(T, 113, 0.576404779858501791742255670403304787e-4),
 839          BOOST_MATH_BIG_CONSTANT(T, 113, 0.196477049872253010859712483984252067e-5),
 840          BOOST_MATH_BIG_CONSTANT(T, 113, 0.521863830500876189501054079974475762e-7),
 841          BOOST_MATH_BIG_CONSTANT(T, 113, 0.109524209196868135198775445228552059e-8),
 842          BOOST_MATH_BIG_CONSTANT(T, 113, 0.181698713448644481083966260949267825e-10),
 843          BOOST_MATH_BIG_CONSTANT(T, 113, 0.234793316975091282090312036524695562e-12),
 844          BOOST_MATH_BIG_CONSTANT(T, 113, 0.227490441461460571047545264251399048e-14),
 845          BOOST_MATH_BIG_CONSTANT(T, 113, 0.151500292036937400913870642638520668e-16),
 846          BOOST_MATH_BIG_CONSTANT(T, 113, 0.543475775154780935815530649335936121e-19),
 847          BOOST_MATH_BIG_CONSTANT(T, 113, 0.241647013434111434636554455083309352e-28),
 848          BOOST_MATH_BIG_CONSTANT(T, 113, -0.557103423021951053707162364713587374e-31),
 849          BOOST_MATH_BIG_CONSTANT(T, 113, 0.618708773442584843384712258199645166e-34),
 850       };
 851       result = tools::evaluate_polynomial(P, T(s - 30)) / tools::evaluate_polynomial(Q, T(s - 30));
 852       result = 1 + exp(result);
 853    }
 854    else if(s < 117)
 855    {
 856       result = 1 + pow(T(2), -s);
 857    }
 858    else
 859    {
 860       result = 1;
 861    }
 862    return result;
 863 }
 864
 865 template <class T, class Policy, class Tag>
 866 T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
 867 {
 868    BOOST_MATH_STD_USING
 869    if(sc == 0)
 870       return policies::raise_pole_error<T>(
 871          "boost::math::zeta<%1%>",
 872          "Evaluation of zeta function at pole %1%",
 873          s, pol);
 874    T result;
 875    if(fabs(s) < tools::root_epsilon<T>())
 876    {
 877       result = -0.5f - constants::log_root_two_pi<T, Policy>() * s;
 878    }
 879    else if(s < 0)
 880    {
 881       std::swap(s, sc);
 882       if(floor(sc/2) == sc/2)
 883          result = 0;
 884       else
 885       {
 886          result = boost::math::sin_pi(0.5f * sc, pol)
 887             * 2 * pow(2 * constants::pi<T>(), -s)
 888             * boost::math::tgamma(s, pol)
 889             * zeta_imp(s, sc, pol, tag);
 890       }
 891    }
 892    else
 893    {
 894       result = zeta_imp_prec(s, sc, pol, tag);
 895    }
 896    return result;
 897 }
 898
 899 template <class T, class Policy, class tag>
 900 struct zeta_initializer
 901 {
 902    struct init
 903    {
 904       init()
 905       {
 906          do_init(tag());
 907       }
 908       static void do_init(const mpl::int_<0>&){}
 909       static void do_init(const mpl::int_<53>&){}
 910       static void do_init(const mpl::int_<64>&)
 911       {
 912          boost::math::zeta(static_cast<T>(0.5), Policy());
 913          boost::math::zeta(static_cast<T>(1.5), Policy());
 914          boost::math::zeta(static_cast<T>(3.5), Policy());
 915          boost::math::zeta(static_cast<T>(6.5), Policy());
 916          boost::math::zeta(static_cast<T>(14.5), Policy());
 917          boost::math::zeta(static_cast<T>(40.5), Policy());
 918       }
 919       static void do_init(const mpl::int_<113>&)
 920       {
 921          boost::math::zeta(static_cast<T>(0.5), Policy());
 922          boost::math::zeta(static_cast<T>(1.5), Policy());
 923          boost::math::zeta(static_cast<T>(3.5), Policy());
 924          boost::math::zeta(static_cast<T>(5.5), Policy());
 925          boost::math::zeta(static_cast<T>(9.5), Policy());
 926          boost::math::zeta(static_cast<T>(16.5), Policy());
 927          boost::math::zeta(static_cast<T>(25), Policy());
 928          boost::math::zeta(static_cast<T>(70), Policy());
 929       }
 930       void force_instantiate()const{}
 931    };
 932    static const init initializer;
 933    static void force_instantiate()
 934    {
 935       initializer.force_instantiate();
 936    }
 937 };
 938
 939 template <class T, class Policy, class tag>
 940 const typename zeta_initializer<T, Policy, tag>::init zeta_initializer<T, Policy, tag>::initializer;
 941
 942 } // detail
 943
 944 template <class T, class Policy>
 945 inline typename tools::promote_args<T>::type zeta(T s, const Policy&)
 946 {
 947    typedef typename tools::promote_args<T>::type result_type;
 948    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 949    typedef typename policies::precision<result_type, Policy>::type precision_type;
 950    typedef typename policies::normalise<
 951       Policy,
 952       policies::promote_float<false>,
 953       policies::promote_double<false>,
 954       policies::discrete_quantile<>,
 955       policies::assert_undefined<> >::type forwarding_policy;
 956    typedef typename mpl::if_<
 957       mpl::less_equal<precision_type, mpl::int_<0> >,
 958       mpl::int_<0>,
 959       typename mpl::if_<
 960          mpl::less_equal<precision_type, mpl::int_<53> >,
 961          mpl::int_<53>,  // double
 962          typename mpl::if_<
 963             mpl::less_equal<precision_type, mpl::int_<64> >,
 964             mpl::int_<64>, // 80-bit long double
 965             typename mpl::if_<
 966                mpl::less_equal<precision_type, mpl::int_<113> >,
 967                mpl::int_<113>, // 128-bit long double
 968                mpl::int_<0> // too many bits, use generic version.
 969             >::type
 970          >::type
 971       >::type
 972    >::type tag_type;
 973    //typedef mpl::int_<0> tag_type;
 974
 975    detail::zeta_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
 976
 977    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(
 978       static_cast<value_type>(s),
 979       static_cast<value_type>(1 - static_cast<value_type>(s)),
 980       forwarding_policy(),
 981       tag_type()), "boost::math::zeta<%1%>(%1%)");
 982 }
 983
 984 template <class T>
 985 inline typename tools::promote_args<T>::type zeta(T s)
 986 {
 987    return zeta(s, policies::policy<>());
 988 }
 989
 990 }} // namespaces
 991
 992 #endif // BOOST_MATH_ZETA_HPP
 993
 994
 995