1082 lines · plain
1// Copyright John Maddock 2007, 2014.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_ZETA_HPP7#define BOOST_MATH_ZETA_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/special_functions/math_fwd.hpp>14#include <boost/math/tools/precision.hpp>15#include <boost/math/tools/series.hpp>16#include <boost/math/tools/big_constant.hpp>17#include <boost/math/policies/error_handling.hpp>18#include <boost/math/special_functions/gamma.hpp>19#include <boost/math/special_functions/factorials.hpp>20#include <boost/math/special_functions/sin_pi.hpp>21 22#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)23//24// This is the only way we can avoid25// warning: non-standard suffix on floating constant [-Wpedantic]26// when building with -Wall -pedantic. Neither __extension__27// nor #pragma diagnostic ignored work :(28//29#pragma GCC system_header30#endif31 32namespace boost{ namespace math{ namespace detail{33 34#if 035//36// This code is commented out because we have a better more rapidly converging series37// now. Retained for future reference and in case the new code causes any issues down the line....38//39 40template <class T, class Policy>41struct zeta_series_cache_size42{43 //44 // Work how large to make our cache size when evaluating the series45 // evaluation: normally this is just large enough for the series46 // to have converged, but for arbitrary precision types we need a47 // really large cache to achieve reasonable precision in a reasonable48 // time. This is important when constructing rational approximations49 // to zeta for example.50 //51 typedef typename boost::math::policies::precision<T,Policy>::type precision_type;52 typedef typename mpl::if_<53 mpl::less_equal<precision_type, std::integral_constant<int, 0> >,54 std::integral_constant<int, 5000>,55 typename mpl::if_<56 mpl::less_equal<precision_type, std::integral_constant<int, 64> >,57 std::integral_constant<int, 70>,58 typename mpl::if_<59 mpl::less_equal<precision_type, std::integral_constant<int, 113> >,60 std::integral_constant<int, 100>,61 std::integral_constant<int, 5000>62 >::type63 >::type64 >::type type;65};66 67template <class T, class Policy>68T zeta_series_imp(T s, T sc, const Policy&)69{70 //71 // Series evaluation from:72 // Havil, J. Gamma: Exploring Euler's Constant.73 // Princeton, NJ: Princeton University Press, 2003.74 //75 // See also http://mathworld.wolfram.com/RiemannZetaFunction.html76 //77 BOOST_MATH_STD_USING78 T sum = 0;79 T mult = 0.5;80 T change;81 typedef typename zeta_series_cache_size<T,Policy>::type cache_size;82 T powers[cache_size::value] = { 0, };83 unsigned n = 0;84 do{85 T binom = -static_cast<T>(n);86 T nested_sum = 1;87 if(n < sizeof(powers) / sizeof(powers[0]))88 powers[n] = pow(static_cast<T>(n + 1), -s);89 for(unsigned k = 1; k <= n; ++k)90 {91 T p;92 if(k < sizeof(powers) / sizeof(powers[0]))93 {94 p = powers[k];95 //p = pow(k + 1, -s);96 }97 else98 p = pow(static_cast<T>(k + 1), -s);99 nested_sum += binom * p;100 binom *= (k - static_cast<T>(n)) / (k + 1);101 }102 change = mult * nested_sum;103 sum += change;104 mult /= 2;105 ++n;106 }while(fabs(change / sum) > tools::epsilon<T>());107 108 return sum * 1 / -boost::math::powm1(T(2), sc);109}110 111//112// Classical p-series:113//114template <class T>115struct zeta_series2116{117 typedef T result_type;118 zeta_series2(T _s) : s(-_s), k(1){}119 T operator()()120 {121 BOOST_MATH_STD_USING122 return pow(static_cast<T>(k++), s);123 }124private:125 T s;126 unsigned k;127};128 129template <class T, class Policy>130inline T zeta_series2_imp(T s, const Policy& pol)131{132 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();;133 zeta_series2<T> f(s);134 T result = tools::sum_series(135 f,136 policies::get_epsilon<T, Policy>(),137 max_iter);138 policies::check_series_iterations<T>("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);139 return result;140}141#endif142 143template <class T, class Policy>144T zeta_polynomial_series(T s, T sc, Policy const &)145{146 //147 // This is algorithm 3 from:148 //149 // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein,150 // Canadian Mathematical Society, Conference Proceedings.151 // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf152 //153 BOOST_MATH_STD_USING154 int n = itrunc(T(log(boost::math::tools::epsilon<T>()) / -2));155 T sum = 0; // LCOV_EXCL_LINE spurious miss as surrounding lines hit.156 T two_n = ldexp(T(1), n);157 int ej_sign = 1; // LCOV_EXCL_LINE spurious miss as surrounding lines hit.158 for(int j = 0; j < n; ++j)159 {160 sum += ej_sign * -two_n / pow(T(j + 1), s);161 ej_sign = -ej_sign;162 }163 T ej_sum = 1; // LCOV_EXCL_LINE spurious miss as surrounding lines hit.164 T ej_term = 1; // LCOV_EXCL_LINE spurious miss as surrounding lines hit.165 for(int j = n; j <= 2 * n - 1; ++j)166 {167 sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s);168 ej_sign = -ej_sign;169 ej_term *= 2 * n - j;170 ej_term /= j - n + 1;171 ej_sum += ej_term; // LCOV_EXCL_LINE spurious miss as surrounding lines hit.172 }173 return -sum / (two_n * (-powm1(T(2), sc)));174}175//176// MP only, verified as covered by the full tests:177// LCOV_EXCL_START178template <class T, class Policy>179T zeta_imp_prec(T s, T sc, const Policy& pol, const std::integral_constant<int, 0>&)180{181 BOOST_MATH_STD_USING182 T result;183 if(s >= policies::digits<T, Policy>())184 return 1;185 result = zeta_polynomial_series(s, sc, pol);186#if 0187 // Old code archived for future reference:188 189 //190 // Only use power series if it will converge in 100191 // iterations or less: the more iterations it consumes192 // the slower convergence becomes so we have to be very193 // careful in it's usage.194 //195 if (s > -log(tools::epsilon<T>()) / 4.5)196 result = detail::zeta_series2_imp(s, pol);197 else198 result = detail::zeta_series_imp(s, sc, pol);199#endif200 return result;201}202// LCOV_EXCL_STOP203 204template <class T, class Policy>205inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 53>&)206{207 BOOST_MATH_STD_USING208 T result;209 if(s < 1)210 {211 // Rational Approximation212 // Maximum Deviation Found: 2.020e-18213 // Expected Error Term: -2.020e-18214 // Max error found at double precision: 3.994987e-17215 // LCOV_EXCL_START216 static const T P[6] = {217 static_cast<T>(0.24339294433593750202L),218 static_cast<T>(-0.49092470516353571651L),219 static_cast<T>(0.0557616214776046784287L),220 static_cast<T>(-0.00320912498879085894856L),221 static_cast<T>(0.000451534528645796438704L),222 static_cast<T>(-0.933241270357061460782e-5L),223 };224 static const T Q[6] = {225 static_cast<T>(1L),226 static_cast<T>(-0.279960334310344432495L),227 static_cast<T>(0.0419676223309986037706L),228 static_cast<T>(-0.00413421406552171059003L),229 static_cast<T>(0.00024978985622317935355L),230 static_cast<T>(-0.101855788418564031874e-4L),231 };232 // LCOV_EXCL_STOP233 result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);234 result -= 1.2433929443359375F;235 result += (sc);236 result /= (sc);237 }238 else if(s <= 2)239 {240 // Maximum Deviation Found: 9.007e-20241 // Expected Error Term: 9.007e-20242 // LCOV_EXCL_START243 static const T P[6] = {244 static_cast<T>(0.577215664901532860516L),245 static_cast<T>(0.243210646940107164097L),246 static_cast<T>(0.0417364673988216497593L),247 static_cast<T>(0.00390252087072843288378L),248 static_cast<T>(0.000249606367151877175456L),249 static_cast<T>(0.110108440976732897969e-4L),250 };251 static const T Q[6] = {252 static_cast<T>(1.0),253 static_cast<T>(0.295201277126631761737L),254 static_cast<T>(0.043460910607305495864L),255 static_cast<T>(0.00434930582085826330659L),256 static_cast<T>(0.000255784226140488490982L),257 static_cast<T>(0.10991819782396112081e-4L),258 };259 // LCOV_EXCL_STOP260 result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));261 result += 1 / (-sc);262 }263 else if(s <= 4)264 {265 // Maximum Deviation Found: 5.946e-22266 // Expected Error Term: -5.946e-22267 // LCOV_EXCL_START268 static const float Y = 0.6986598968505859375;269 static const T P[6] = {270 static_cast<T>(-0.0537258300023595030676L),271 static_cast<T>(0.0445163473292365591906L),272 static_cast<T>(0.0128677673534519952905L),273 static_cast<T>(0.00097541770457391752726L),274 static_cast<T>(0.769875101573654070925e-4L),275 static_cast<T>(0.328032510000383084155e-5L),276 };277 static const T Q[7] = {278 1.0f,279 static_cast<T>(0.33383194553034051422L),280 static_cast<T>(0.0487798431291407621462L),281 static_cast<T>(0.00479039708573558490716L),282 static_cast<T>(0.000270776703956336357707L),283 static_cast<T>(0.106951867532057341359e-4L),284 static_cast<T>(0.236276623974978646399e-7L),285 };286 // LCOV_EXCL_STOP287 result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));288 result += Y + 1 / (-sc);289 }290 else if(s <= 7)291 {292 // Maximum Deviation Found: 2.955e-17293 // Expected Error Term: 2.955e-17294 // Max error found at double precision: 2.009135e-16295 // LCOV_EXCL_START296 static const T P[6] = {297 static_cast<T>(-2.49710190602259410021L),298 static_cast<T>(-2.60013301809475665334L),299 static_cast<T>(-0.939260435377109939261L),300 static_cast<T>(-0.138448617995741530935L),301 static_cast<T>(-0.00701721240549802377623L),302 static_cast<T>(-0.229257310594893932383e-4L),303 };304 static const T Q[9] = {305 1.0f,306 static_cast<T>(0.706039025937745133628L),307 static_cast<T>(0.15739599649558626358L),308 static_cast<T>(0.0106117950976845084417L),309 static_cast<T>(-0.36910273311764618902e-4L),310 static_cast<T>(0.493409563927590008943e-5L),311 static_cast<T>(-0.234055487025287216506e-6L),312 static_cast<T>(0.718833729365459760664e-8L),313 static_cast<T>(-0.1129200113474947419e-9L),314 };315 // LCOV_EXCL_STOP316 result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));317 result = 1 + exp(result);318 }319 else if(s < 15)320 {321 // Maximum Deviation Found: 7.117e-16322 // Expected Error Term: 7.117e-16323 // Max error found at double precision: 9.387771e-16324 // LCOV_EXCL_START325 static const T P[7] = {326 static_cast<T>(-4.78558028495135619286L),327 static_cast<T>(-1.89197364881972536382L),328 static_cast<T>(-0.211407134874412820099L),329 static_cast<T>(-0.000189204758260076688518L),330 static_cast<T>(0.00115140923889178742086L),331 static_cast<T>(0.639949204213164496988e-4L),332 static_cast<T>(0.139348932445324888343e-5L),333 };334 static const T Q[9] = {335 1.0f,336 static_cast<T>(0.244345337378188557777L),337 static_cast<T>(0.00873370754492288653669L),338 static_cast<T>(-0.00117592765334434471562L),339 static_cast<T>(-0.743743682899933180415e-4L),340 static_cast<T>(-0.21750464515767984778e-5L),341 static_cast<T>(0.471001264003076486547e-8L),342 static_cast<T>(-0.833378440625385520576e-10L),343 static_cast<T>(0.699841545204845636531e-12L),344 };345 // LCOV_EXCL_STOP346 result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));347 result = 1 + exp(result);348 }349 else if(s < 36)350 {351 // Max error in interpolated form: 1.668e-17352 // Max error found at long double precision: 1.669714e-17353 // LCOV_EXCL_START354 static const T P[8] = {355 static_cast<T>(-10.3948950573308896825L),356 static_cast<T>(-2.85827219671106697179L),357 static_cast<T>(-0.347728266539245787271L),358 static_cast<T>(-0.0251156064655346341766L),359 static_cast<T>(-0.00119459173416968685689L),360 static_cast<T>(-0.382529323507967522614e-4L),361 static_cast<T>(-0.785523633796723466968e-6L),362 static_cast<T>(-0.821465709095465524192e-8L),363 };364 static const T Q[10] = {365 1.0f,366 static_cast<T>(0.208196333572671890965L),367 static_cast<T>(0.0195687657317205033485L),368 static_cast<T>(0.00111079638102485921877L),369 static_cast<T>(0.408507746266039256231e-4L),370 static_cast<T>(0.955561123065693483991e-6L),371 static_cast<T>(0.118507153474022900583e-7L),372 static_cast<T>(0.222609483627352615142e-14L),373 };374 // LCOV_EXCL_STOP375 result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));376 result = 1 + exp(result);377 }378 else379 {380 result = 1 + pow(T(2), -s);381 }382 return result;383}384 385template <class T, class Policy>386T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 64>&)387{388 BOOST_MATH_STD_USING389 T result;390 if(s < 1)391 {392 // Rational Approximation393 // Maximum Deviation Found: 3.099e-20394 // Expected Error Term: 3.099e-20395 // Max error found at long double precision: 5.890498e-20396 // LCOV_EXCL_START397 static const T P[6] = {398 BOOST_MATH_BIG_CONSTANT(T, 64, 0.243392944335937499969),399 BOOST_MATH_BIG_CONSTANT(T, 64, -0.496837806864865688082),400 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0680008039723709987107),401 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00511620413006619942112),402 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000455369899250053003335),403 BOOST_MATH_BIG_CONSTANT(T, 64, -0.279496685273033761927e-4),404 };405 static const T Q[7] = {406 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),407 BOOST_MATH_BIG_CONSTANT(T, 64, -0.30425480068225790522),408 BOOST_MATH_BIG_CONSTANT(T, 64, 0.050052748580371598736),409 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00519355671064700627862),410 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000360623385771198350257),411 BOOST_MATH_BIG_CONSTANT(T, 64, -0.159600883054550987633e-4),412 BOOST_MATH_BIG_CONSTANT(T, 64, 0.339770279812410586032e-6),413 };414 // LCOV_EXCL_STOP415 result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);416 result -= 1.2433929443359375F;417 result += (sc);418 result /= (sc);419 }420 else if(s <= 2)421 {422 // Maximum Deviation Found: 1.059e-21423 // Expected Error Term: 1.059e-21424 // Max error found at long double precision: 1.626303e-19425 // LCOV_EXCL_START426 static const T P[6] = {427 BOOST_MATH_BIG_CONSTANT(T, 64, 0.577215664901532860605),428 BOOST_MATH_BIG_CONSTANT(T, 64, 0.222537368917162139445),429 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0356286324033215682729),430 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00304465292366350081446),431 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000178102511649069421904),432 BOOST_MATH_BIG_CONSTANT(T, 64, 0.700867470265983665042e-5),433 };434 static const T Q[7] = {435 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),436 BOOST_MATH_BIG_CONSTANT(T, 64, 0.259385759149531030085),437 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0373974962106091316854),438 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00332735159183332820617),439 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000188690420706998606469),440 BOOST_MATH_BIG_CONSTANT(T, 64, 0.635994377921861930071e-5),441 BOOST_MATH_BIG_CONSTANT(T, 64, 0.226583954978371199405e-7),442 };443 // LCOV_EXCL_STOP444 result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));445 result += 1 / (-sc);446 }447 else if(s <= 4)448 {449 // Maximum Deviation Found: 5.946e-22450 // Expected Error Term: -5.946e-22451 // LCOV_EXCL_START452 static const float Y = 0.6986598968505859375;453 static const T P[7] = {454 BOOST_MATH_BIG_CONSTANT(T, 64, -0.053725830002359501027),455 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0470551187571475844778),456 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0101339410415759517471),457 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00100240326666092854528),458 BOOST_MATH_BIG_CONSTANT(T, 64, 0.685027119098122814867e-4),459 BOOST_MATH_BIG_CONSTANT(T, 64, 0.390972820219765942117e-5),460 BOOST_MATH_BIG_CONSTANT(T, 64, 0.540319769113543934483e-7),461 };462 static const T Q[8] = {463 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),464 BOOST_MATH_BIG_CONSTANT(T, 64, 0.286577739726542730421),465 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0447355811517733225843),466 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00430125107610252363302),467 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000284956969089786662045),468 BOOST_MATH_BIG_CONSTANT(T, 64, 0.116188101609848411329e-4),469 BOOST_MATH_BIG_CONSTANT(T, 64, 0.278090318191657278204e-6),470 BOOST_MATH_BIG_CONSTANT(T, 64, -0.19683620233222028478e-8),471 };472 // LCOV_EXCL_STOP473 result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));474 result += Y + 1 / (-sc);475 }476 else if(s <= 7)477 {478 // Max error found at long double precision: 8.132216e-19479 // LCOV_EXCL_START480 static const T P[8] = {481 BOOST_MATH_BIG_CONSTANT(T, 64, -2.49710190602259407065),482 BOOST_MATH_BIG_CONSTANT(T, 64, -3.36664913245960625334),483 BOOST_MATH_BIG_CONSTANT(T, 64, -1.77180020623777595452),484 BOOST_MATH_BIG_CONSTANT(T, 64, -0.464717885249654313933),485 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0643694921293579472583),486 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00464265386202805715487),487 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000165556579779704340166),488 BOOST_MATH_BIG_CONSTANT(T, 64, -0.252884970740994069582e-5),489 };490 static const T Q[9] = {491 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),492 BOOST_MATH_BIG_CONSTANT(T, 64, 1.01300131390690459085),493 BOOST_MATH_BIG_CONSTANT(T, 64, 0.387898115758643503827),494 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0695071490045701135188),495 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00586908595251442839291),496 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000217752974064612188616),497 BOOST_MATH_BIG_CONSTANT(T, 64, 0.397626583349419011731e-5),498 BOOST_MATH_BIG_CONSTANT(T, 64, -0.927884739284359700764e-8),499 BOOST_MATH_BIG_CONSTANT(T, 64, 0.119810501805618894381e-9),500 };501 // LCOV_EXCL_STOP502 result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));503 result = 1 + exp(result);504 }505 else if(s < 15)506 {507 // Max error in interpolated form: 1.133e-18508 // Max error found at long double precision: 2.183198e-18509 // LCOV_EXCL_START510 static const T P[9] = {511 BOOST_MATH_BIG_CONSTANT(T, 64, -4.78558028495135548083),512 BOOST_MATH_BIG_CONSTANT(T, 64, -3.23873322238609358947),513 BOOST_MATH_BIG_CONSTANT(T, 64, -0.892338582881021799922),514 BOOST_MATH_BIG_CONSTANT(T, 64, -0.131326296217965913809),515 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0115651591773783712996),516 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000657728968362695775205),517 BOOST_MATH_BIG_CONSTANT(T, 64, -0.252051328129449973047e-4),518 BOOST_MATH_BIG_CONSTANT(T, 64, -0.626503445372641798925e-6),519 BOOST_MATH_BIG_CONSTANT(T, 64, -0.815696314790853893484e-8),520 };521 static const T Q[9] = {522 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),523 BOOST_MATH_BIG_CONSTANT(T, 64, 0.525765665400123515036),524 BOOST_MATH_BIG_CONSTANT(T, 64, 0.10852641753657122787),525 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0115669945375362045249),526 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000732896513858274091966),527 BOOST_MATH_BIG_CONSTANT(T, 64, 0.30683952282420248448e-4),528 BOOST_MATH_BIG_CONSTANT(T, 64, 0.819649214609633126119e-6),529 BOOST_MATH_BIG_CONSTANT(T, 64, 0.117957556472335968146e-7),530 BOOST_MATH_BIG_CONSTANT(T, 64, -0.193432300973017671137e-12),531 };532 // LCOV_EXCL_STOP533 result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7));534 result = 1 + exp(result);535 }536 else if(s < 42)537 {538 // Max error in interpolated form: 1.668e-17539 // Max error found at long double precision: 1.669714e-17540 // LCOV_EXCL_START541 static const T P[9] = {542 BOOST_MATH_BIG_CONSTANT(T, 64, -10.3948950573308861781),543 BOOST_MATH_BIG_CONSTANT(T, 64, -2.82646012777913950108),544 BOOST_MATH_BIG_CONSTANT(T, 64, -0.342144362739570333665),545 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0249285145498722647472),546 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00122493108848097114118),547 BOOST_MATH_BIG_CONSTANT(T, 64, -0.423055371192592850196e-4),548 BOOST_MATH_BIG_CONSTANT(T, 64, -0.1025215577185967488e-5),549 BOOST_MATH_BIG_CONSTANT(T, 64, -0.165096762663509467061e-7),550 BOOST_MATH_BIG_CONSTANT(T, 64, -0.145392555873022044329e-9),551 };552 static const T Q[10] = {553 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),554 BOOST_MATH_BIG_CONSTANT(T, 64, 0.205135978585281988052),555 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0192359357875879453602),556 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00111496452029715514119),557 BOOST_MATH_BIG_CONSTANT(T, 64, 0.434928449016693986857e-4),558 BOOST_MATH_BIG_CONSTANT(T, 64, 0.116911068726610725891e-5),559 BOOST_MATH_BIG_CONSTANT(T, 64, 0.206704342290235237475e-7),560 BOOST_MATH_BIG_CONSTANT(T, 64, 0.209772836100827647474e-9),561 BOOST_MATH_BIG_CONSTANT(T, 64, -0.939798249922234703384e-16),562 BOOST_MATH_BIG_CONSTANT(T, 64, 0.264584017421245080294e-18),563 };564 // LCOV_EXCL_STOP565 result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15));566 result = 1 + exp(result);567 }568 else569 {570 result = 1 + pow(T(2), -s);571 }572 return result;573}574 575template <class T, class Policy>576T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant<int, 113>&)577{578 BOOST_MATH_STD_USING579 T result;580 if(s < 1)581 {582 // Rational Approximation583 // Maximum Deviation Found: 9.493e-37584 // Expected Error Term: 9.492e-37585 // Max error found at long double precision: 7.281332e-31586 // LCOV_EXCL_START587 static const T P[10] = {588 BOOST_MATH_BIG_CONSTANT(T, 113, -1.0),589 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0353008629988648122808504280990313668),590 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0107795651204927743049369868548706909),591 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000523961870530500751114866884685172975),592 BOOST_MATH_BIG_CONSTANT(T, 113, -0.661805838304910731947595897966487515e-4),593 BOOST_MATH_BIG_CONSTANT(T, 113, -0.658932670403818558510656304189164638e-5),594 BOOST_MATH_BIG_CONSTANT(T, 113, -0.103437265642266106533814021041010453e-6),595 BOOST_MATH_BIG_CONSTANT(T, 113, 0.116818787212666457105375746642927737e-7),596 BOOST_MATH_BIG_CONSTANT(T, 113, 0.660690993901506912123512551294239036e-9),597 BOOST_MATH_BIG_CONSTANT(T, 113, 0.113103113698388531428914333768142527e-10),598 };599 static const T Q[11] = {600 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),601 BOOST_MATH_BIG_CONSTANT(T, 113, -0.387483472099602327112637481818565459),602 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0802265315091063135271497708694776875),603 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0110727276164171919280036408995078164),604 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00112552716946286252000434849173787243),605 BOOST_MATH_BIG_CONSTANT(T, 113, -0.874554160748626916455655180296834352e-4),606 BOOST_MATH_BIG_CONSTANT(T, 113, 0.530097847491828379568636739662278322e-5),607 BOOST_MATH_BIG_CONSTANT(T, 113, -0.248461553590496154705565904497247452e-6),608 BOOST_MATH_BIG_CONSTANT(T, 113, 0.881834921354014787309644951507523899e-8),609 BOOST_MATH_BIG_CONSTANT(T, 113, -0.217062446168217797598596496310953025e-9),610 BOOST_MATH_BIG_CONSTANT(T, 113, 0.315823200002384492377987848307151168e-11),611 };612 // LCOV_EXCL_STOP613 result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);614 result += (sc);615 result /= (sc);616 }617 else if(s <= 2)618 {619 // Maximum Deviation Found: 1.616e-37620 // Expected Error Term: -1.615e-37621 // LCOV_EXCL_START622 static const T P[10] = {623 BOOST_MATH_BIG_CONSTANT(T, 113, 0.577215664901532860606512090082402431),624 BOOST_MATH_BIG_CONSTANT(T, 113, 0.255597968739771510415479842335906308),625 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0494056503552807274142218876983542205),626 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00551372778611700965268920983472292325),627 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00043667616723970574871427830895192731),628 BOOST_MATH_BIG_CONSTANT(T, 113, 0.268562259154821957743669387915239528e-4),629 BOOST_MATH_BIG_CONSTANT(T, 113, 0.109249633923016310141743084480436612e-5),630 BOOST_MATH_BIG_CONSTANT(T, 113, 0.273895554345300227466534378753023924e-7),631 BOOST_MATH_BIG_CONSTANT(T, 113, 0.583103205551702720149237384027795038e-9),632 BOOST_MATH_BIG_CONSTANT(T, 113, -0.835774625259919268768735944711219256e-11),633 };634 static const T Q[11] = {635 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),636 BOOST_MATH_BIG_CONSTANT(T, 113, 0.316661751179735502065583176348292881),637 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0540401806533507064453851182728635272),638 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00598621274107420237785899476374043797),639 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000474907812321704156213038740142079615),640 BOOST_MATH_BIG_CONSTANT(T, 113, 0.272125421722314389581695715835862418e-4),641 BOOST_MATH_BIG_CONSTANT(T, 113, 0.112649552156479800925522445229212933e-5),642 BOOST_MATH_BIG_CONSTANT(T, 113, 0.301838975502992622733000078063330461e-7),643 BOOST_MATH_BIG_CONSTANT(T, 113, 0.422960728687211282539769943184270106e-9),644 BOOST_MATH_BIG_CONSTANT(T, 113, -0.377105263588822468076813329270698909e-11),645 BOOST_MATH_BIG_CONSTANT(T, 113, -0.581926559304525152432462127383600681e-13),646 };647 // LCOV_EXCL_STOP648 result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc));649 result += 1 / (-sc);650 }651 else if(s <= 4)652 {653 // Maximum Deviation Found: 1.891e-36654 // Expected Error Term: -1.891e-36655 // Max error found: 2.171527e-35656 // LCOV_EXCL_START657 static const float Y = 0.6986598968505859375;658 static const T P[11] = {659 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0537258300023595010275848333539748089),660 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0429086930802630159457448174466342553),661 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0136148228754303412510213395034056857),662 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00190231601036042925183751238033763915),663 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000186880390916311438818302549192456581),664 BOOST_MATH_BIG_CONSTANT(T, 113, 0.145347370745893262394287982691323657e-4),665 BOOST_MATH_BIG_CONSTANT(T, 113, 0.805843276446813106414036600485884885e-6),666 BOOST_MATH_BIG_CONSTANT(T, 113, 0.340818159286739137503297172091882574e-7),667 BOOST_MATH_BIG_CONSTANT(T, 113, 0.115762357488748996526167305116837246e-8),668 BOOST_MATH_BIG_CONSTANT(T, 113, 0.231904754577648077579913403645767214e-10),669 BOOST_MATH_BIG_CONSTANT(T, 113, 0.340169592866058506675897646629036044e-12),670 };671 static const T Q[12] = {672 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),673 BOOST_MATH_BIG_CONSTANT(T, 113, 0.363755247765087100018556983050520554),674 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0696581979014242539385695131258321598),675 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00882208914484611029571547753782014817),676 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000815405623261946661762236085660996718),677 BOOST_MATH_BIG_CONSTANT(T, 113, 0.571366167062457197282642344940445452e-4),678 BOOST_MATH_BIG_CONSTANT(T, 113, 0.309278269271853502353954062051797838e-5),679 BOOST_MATH_BIG_CONSTANT(T, 113, 0.12822982083479010834070516053794262e-6),680 BOOST_MATH_BIG_CONSTANT(T, 113, 0.397876357325018976733953479182110033e-8),681 BOOST_MATH_BIG_CONSTANT(T, 113, 0.8484432107648683277598472295289279e-10),682 BOOST_MATH_BIG_CONSTANT(T, 113, 0.105677416606909614301995218444080615e-11),683 BOOST_MATH_BIG_CONSTANT(T, 113, 0.547223964564003701979951154093005354e-15),684 };685 // LCOV_EXCL_STOP686 result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2));687 result += Y + 1 / (-sc);688 }689 else if(s <= 6)690 {691 // Max error in interpolated form: 1.510e-37692 // Max error found at long double precision: 2.769266e-34693 // LCOV_EXCL_START694 static const T Y = 3.28348541259765625F;695 696 static const T P[13] = {697 BOOST_MATH_BIG_CONSTANT(T, 113, 0.786383506575062179339611614117697622),698 BOOST_MATH_BIG_CONSTANT(T, 113, 0.495766593395271370974685959652073976),699 BOOST_MATH_BIG_CONSTANT(T, 113, -0.409116737851754766422360889037532228),700 BOOST_MATH_BIG_CONSTANT(T, 113, -0.57340744006238263817895456842655987),701 BOOST_MATH_BIG_CONSTANT(T, 113, -0.280479899797421910694892949057963111),702 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0753148409447590257157585696212649869),703 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0122934003684672788499099362823748632),704 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126148398446193639247961370266962927),705 BOOST_MATH_BIG_CONSTANT(T, 113, -0.828465038179772939844657040917364896e-4),706 BOOST_MATH_BIG_CONSTANT(T, 113, -0.361008916706050977143208468690645684e-5),707 BOOST_MATH_BIG_CONSTANT(T, 113, -0.109879825497910544424797771195928112e-6),708 BOOST_MATH_BIG_CONSTANT(T, 113, -0.214539416789686920918063075528797059e-8),709 BOOST_MATH_BIG_CONSTANT(T, 113, -0.15090220092460596872172844424267351e-10),710 };711 static const T Q[14] = {712 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),713 BOOST_MATH_BIG_CONSTANT(T, 113, 1.69490865837142338462982225731926485),714 BOOST_MATH_BIG_CONSTANT(T, 113, 1.22697696630994080733321401255942464),715 BOOST_MATH_BIG_CONSTANT(T, 113, 0.495409420862526540074366618006341533),716 BOOST_MATH_BIG_CONSTANT(T, 113, 0.122368084916843823462872905024259633),717 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0191412993625268971656513890888208623),718 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00191401538628980617753082598351559642),719 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000123318142456272424148930280876444459),720 BOOST_MATH_BIG_CONSTANT(T, 113, 0.531945488232526067889835342277595709e-5),721 BOOST_MATH_BIG_CONSTANT(T, 113, 0.161843184071894368337068779669116236e-6),722 BOOST_MATH_BIG_CONSTANT(T, 113, 0.305796079600152506743828859577462778e-8),723 BOOST_MATH_BIG_CONSTANT(T, 113, 0.233582592298450202680170811044408894e-10),724 BOOST_MATH_BIG_CONSTANT(T, 113, -0.275363878344548055574209713637734269e-13),725 BOOST_MATH_BIG_CONSTANT(T, 113, 0.221564186807357535475441900517843892e-15),726 };727 // LCOV_EXCL_STOP728 result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4));729 result -= Y;730 result = 1 + exp(result);731 }732 else if(s < 10)733 {734 // Max error in interpolated form: 1.999e-34735 // Max error found at long double precision: 2.156186e-33736 // LCOV_EXCL_START737 static const T P[13] = {738 BOOST_MATH_BIG_CONSTANT(T, 113, -4.0545627381873738086704293881227365),739 BOOST_MATH_BIG_CONSTANT(T, 113, -4.70088348734699134347906176097717782),740 BOOST_MATH_BIG_CONSTANT(T, 113, -2.36921550900925512951976617607678789),741 BOOST_MATH_BIG_CONSTANT(T, 113, -0.684322583796369508367726293719322866),742 BOOST_MATH_BIG_CONSTANT(T, 113, -0.126026534540165129870721937592996324),743 BOOST_MATH_BIG_CONSTANT(T, 113, -0.015636903921778316147260572008619549),744 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00135442294754728549644376325814460807),745 BOOST_MATH_BIG_CONSTANT(T, 113, -0.842793965853572134365031384646117061e-4),746 BOOST_MATH_BIG_CONSTANT(T, 113, -0.385602133791111663372015460784978351e-5),747 BOOST_MATH_BIG_CONSTANT(T, 113, -0.130458500394692067189883214401478539e-6),748 BOOST_MATH_BIG_CONSTANT(T, 113, -0.315861074947230418778143153383660035e-8),749 BOOST_MATH_BIG_CONSTANT(T, 113, -0.500334720512030826996373077844707164e-10),750 BOOST_MATH_BIG_CONSTANT(T, 113, -0.420204769185233365849253969097184005e-12),751 };752 static const T Q[14] = {753 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),754 BOOST_MATH_BIG_CONSTANT(T, 113, 0.97663511666410096104783358493318814),755 BOOST_MATH_BIG_CONSTANT(T, 113, 0.40878780231201806504987368939673249),756 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0963890666609396058945084107597727252),757 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0142207619090854604824116070866614505),758 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139010220902667918476773423995750877),759 BOOST_MATH_BIG_CONSTANT(T, 113, 0.940669540194694997889636696089994734e-4),760 BOOST_MATH_BIG_CONSTANT(T, 113, 0.458220848507517004399292480807026602e-5),761 BOOST_MATH_BIG_CONSTANT(T, 113, 0.16345521617741789012782420625435495e-6),762 BOOST_MATH_BIG_CONSTANT(T, 113, 0.414007452533083304371566316901024114e-8),763 BOOST_MATH_BIG_CONSTANT(T, 113, 0.68701473543366328016953742622661377e-10),764 BOOST_MATH_BIG_CONSTANT(T, 113, 0.603461891080716585087883971886075863e-12),765 BOOST_MATH_BIG_CONSTANT(T, 113, 0.294670713571839023181857795866134957e-16),766 BOOST_MATH_BIG_CONSTANT(T, 113, -0.147003914536437243143096875069813451e-18),767 };768 // LCOV_EXCL_STOP769 result = tools::evaluate_polynomial(P, T(s - 6)) / tools::evaluate_polynomial(Q, T(s - 6));770 result = 1 + exp(result);771 }772 else if(s < 17)773 {774 // Max error in interpolated form: 1.641e-32775 // Max error found at long double precision: 1.696121e-32776 // LCOV_EXCL_START777 static const T P[13] = {778 BOOST_MATH_BIG_CONSTANT(T, 113, -6.91319491921722925920883787894829678),779 BOOST_MATH_BIG_CONSTANT(T, 113, -3.65491257639481960248690596951049048),780 BOOST_MATH_BIG_CONSTANT(T, 113, -0.813557553449954526442644544105257881),781 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0994317301685870959473658713841138083),782 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00726896610245676520248617014211734906),783 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000317253318715075854811266230916762929),784 BOOST_MATH_BIG_CONSTANT(T, 113, -0.66851422826636750855184211580127133e-5),785 BOOST_MATH_BIG_CONSTANT(T, 113, 0.879464154730985406003332577806849971e-7),786 BOOST_MATH_BIG_CONSTANT(T, 113, 0.113838903158254250631678791998294628e-7),787 BOOST_MATH_BIG_CONSTANT(T, 113, 0.379184410304927316385211327537817583e-9),788 BOOST_MATH_BIG_CONSTANT(T, 113, 0.612992858643904887150527613446403867e-11),789 BOOST_MATH_BIG_CONSTANT(T, 113, 0.347873737198164757035457841688594788e-13),790 BOOST_MATH_BIG_CONSTANT(T, 113, -0.289187187441625868404494665572279364e-15),791 };792 static const T Q[14] = {793 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),794 BOOST_MATH_BIG_CONSTANT(T, 113, 0.427310044448071818775721584949868806),795 BOOST_MATH_BIG_CONSTANT(T, 113, 0.074602514873055756201435421385243062),796 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00688651562174480772901425121653945942),797 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000360174847635115036351323894321880445),798 BOOST_MATH_BIG_CONSTANT(T, 113, 0.973556847713307543918865405758248777e-5),799 BOOST_MATH_BIG_CONSTANT(T, 113, -0.853455848314516117964634714780874197e-8),800 BOOST_MATH_BIG_CONSTANT(T, 113, -0.118203513654855112421673192194622826e-7),801 BOOST_MATH_BIG_CONSTANT(T, 113, -0.462521662511754117095006543363328159e-9),802 BOOST_MATH_BIG_CONSTANT(T, 113, -0.834212591919475633107355719369463143e-11),803 BOOST_MATH_BIG_CONSTANT(T, 113, -0.5354594751002702935740220218582929e-13),804 BOOST_MATH_BIG_CONSTANT(T, 113, 0.406451690742991192964889603000756203e-15),805 BOOST_MATH_BIG_CONSTANT(T, 113, 0.887948682401000153828241615760146728e-19),806 BOOST_MATH_BIG_CONSTANT(T, 113, -0.34980761098820347103967203948619072e-21),807 };808 // LCOV_EXCL_STOP809 result = tools::evaluate_polynomial(P, T(s - 10)) / tools::evaluate_polynomial(Q, T(s - 10));810 result = 1 + exp(result);811 }812 else if(s < 30)813 {814 // Max error in interpolated form: 1.563e-31815 // Max error found at long double precision: 1.562725e-31816 // LCOV_EXCL_START817 static const T P[13] = {818 BOOST_MATH_BIG_CONSTANT(T, 113, -11.7824798233959252791987402769438322),819 BOOST_MATH_BIG_CONSTANT(T, 113, -4.36131215284987731928174218354118102),820 BOOST_MATH_BIG_CONSTANT(T, 113, -0.732260980060982349410898496846972204),821 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0744985185694913074484248803015717388),822 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00517228281320594683022294996292250527),823 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000260897206152101522569969046299309939),824 BOOST_MATH_BIG_CONSTANT(T, 113, -0.989553462123121764865178453128769948e-5),825 BOOST_MATH_BIG_CONSTANT(T, 113, -0.286916799741891410827712096608826167e-6),826 BOOST_MATH_BIG_CONSTANT(T, 113, -0.637262477796046963617949532211619729e-8),827 BOOST_MATH_BIG_CONSTANT(T, 113, -0.106796831465628373325491288787760494e-9),828 BOOST_MATH_BIG_CONSTANT(T, 113, -0.129343095511091870860498356205376823e-11),829 BOOST_MATH_BIG_CONSTANT(T, 113, -0.102397936697965977221267881716672084e-13),830 BOOST_MATH_BIG_CONSTANT(T, 113, -0.402663128248642002351627980255756363e-16),831 };832 static const T Q[14] = {833 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),834 BOOST_MATH_BIG_CONSTANT(T, 113, 0.311288325355705609096155335186466508),835 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0438318468940415543546769437752132748),836 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00374396349183199548610264222242269536),837 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000218707451200585197339671707189281302),838 BOOST_MATH_BIG_CONSTANT(T, 113, 0.927578767487930747532953583797351219e-5),839 BOOST_MATH_BIG_CONSTANT(T, 113, 0.294145760625753561951137473484889639e-6),840 BOOST_MATH_BIG_CONSTANT(T, 113, 0.704618586690874460082739479535985395e-8),841 BOOST_MATH_BIG_CONSTANT(T, 113, 0.126333332872897336219649130062221257e-9),842 BOOST_MATH_BIG_CONSTANT(T, 113, 0.16317315713773503718315435769352765e-11),843 BOOST_MATH_BIG_CONSTANT(T, 113, 0.137846712823719515148344938160275695e-13),844 BOOST_MATH_BIG_CONSTANT(T, 113, 0.580975420554224366450994232723910583e-16),845 BOOST_MATH_BIG_CONSTANT(T, 113, -0.291354445847552426900293580511392459e-22),846 BOOST_MATH_BIG_CONSTANT(T, 113, 0.73614324724785855925025452085443636e-25),847 };848 // LCOV_EXCL_STOP849 result = tools::evaluate_polynomial(P, T(s - 17)) / tools::evaluate_polynomial(Q, T(s - 17));850 result = 1 + exp(result);851 }852 else if(s < 74)853 {854 // Max error in interpolated form: 2.311e-27855 // Max error found at long double precision: 2.297544e-27856 // LCOV_EXCL_START857 static const T P[14] = {858 BOOST_MATH_BIG_CONSTANT(T, 113, -20.7944102007844314586649688802236072),859 BOOST_MATH_BIG_CONSTANT(T, 113, -4.95759941987499442499908748130192187),860 BOOST_MATH_BIG_CONSTANT(T, 113, -0.563290752832461751889194629200298688),861 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0406197001137935911912457120706122877),862 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0020846534789473022216888863613422293),863 BOOST_MATH_BIG_CONSTANT(T, 113, -0.808095978462109173749395599401375667e-4),864 BOOST_MATH_BIG_CONSTANT(T, 113, -0.244706022206249301640890603610060959e-5),865 BOOST_MATH_BIG_CONSTANT(T, 113, -0.589477682919645930544382616501666572e-7),866 BOOST_MATH_BIG_CONSTANT(T, 113, -0.113699573675553496343617442433027672e-8),867 BOOST_MATH_BIG_CONSTANT(T, 113, -0.174767860183598149649901223128011828e-10),868 BOOST_MATH_BIG_CONSTANT(T, 113, -0.210051620306761367764549971980026474e-12),869 BOOST_MATH_BIG_CONSTANT(T, 113, -0.189187969537370950337212675466400599e-14),870 BOOST_MATH_BIG_CONSTANT(T, 113, -0.116313253429564048145641663778121898e-16),871 BOOST_MATH_BIG_CONSTANT(T, 113, -0.376708747782400769427057630528578187e-19),872 };873 static const T Q[16] = {874 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),875 BOOST_MATH_BIG_CONSTANT(T, 113, 0.205076752981410805177554569784219717),876 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0202526722696670378999575738524540269),877 BOOST_MATH_BIG_CONSTANT(T, 113, 0.001278305290005994980069466658219057),878 BOOST_MATH_BIG_CONSTANT(T, 113, 0.576404779858501791742255670403304787e-4),879 BOOST_MATH_BIG_CONSTANT(T, 113, 0.196477049872253010859712483984252067e-5),880 BOOST_MATH_BIG_CONSTANT(T, 113, 0.521863830500876189501054079974475762e-7),881 BOOST_MATH_BIG_CONSTANT(T, 113, 0.109524209196868135198775445228552059e-8),882 BOOST_MATH_BIG_CONSTANT(T, 113, 0.181698713448644481083966260949267825e-10),883 BOOST_MATH_BIG_CONSTANT(T, 113, 0.234793316975091282090312036524695562e-12),884 BOOST_MATH_BIG_CONSTANT(T, 113, 0.227490441461460571047545264251399048e-14),885 BOOST_MATH_BIG_CONSTANT(T, 113, 0.151500292036937400913870642638520668e-16),886 BOOST_MATH_BIG_CONSTANT(T, 113, 0.543475775154780935815530649335936121e-19),887 BOOST_MATH_BIG_CONSTANT(T, 113, 0.241647013434111434636554455083309352e-28),888 BOOST_MATH_BIG_CONSTANT(T, 113, -0.557103423021951053707162364713587374e-31),889 BOOST_MATH_BIG_CONSTANT(T, 113, 0.618708773442584843384712258199645166e-34),890 };891 // LCOV_EXCL_STOP892 result = tools::evaluate_polynomial(P, T(s - 30)) / tools::evaluate_polynomial(Q, T(s - 30));893 result = 1 + exp(result);894 }895 else896 {897 result = 1 + pow(T(2), -s);898 }899 return result;900}901 902template <class T, class Policy>903T zeta_imp_odd_integer(int s, const T&, const Policy&, const std::true_type&)904{905 // LCOV_EXCL_START906 static const T results[] = {907 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),908 };909 // LCOV_EXCL_STOP910 return s > 113 ? 1 : results[(s - 3) / 2];911}912 913template <class T, class Policy>914T zeta_imp_odd_integer(int s, const T& sc, const Policy& pol, const std::false_type&)915{916#ifdef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES917 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/");918#endif919 // LCOV_EXCL_START920 static BOOST_MATH_THREAD_LOCAL bool is_init = false;921 static BOOST_MATH_THREAD_LOCAL T results[50] = {};922 static BOOST_MATH_THREAD_LOCAL int digits = tools::digits<T>();923 // LCOV_EXCL_STOP924 int current_digits = tools::digits<T>(); // LCOV_EXCL_LINE spurious miss as surrounding lines hit.925 if(digits != current_digits)926 {927 // Oh my precision has changed...928 is_init = false; // LCOV_EXCL_LINE variable precision MP case only, not included in coverage tests.929 }930 if(!is_init)931 {932 is_init = true;933 digits = current_digits;934 for(unsigned k = 0; k < sizeof(results) / sizeof(results[0]); ++k)935 {936 T arg = k * 2 + 3;937 T c_arg = 1 - arg; // LCOV_EXCL_LINE spurious miss as surrounding lines hit.938 results[k] = zeta_polynomial_series(arg, c_arg, pol);939 }940 }941 const unsigned index = static_cast<unsigned>((s - 3) / 2);942 return index >= sizeof(results) / sizeof(results[0]) ? zeta_polynomial_series(T(s), sc, pol): results[index];943}944 945template <class T, class Policy, class Tag>946T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)947{948 BOOST_MATH_STD_USING949 static const char* function = "boost::math::zeta<%1%>";950 if(sc == 0)951 return policies::raise_pole_error<T>(function, "Evaluation of zeta function at pole %1%", s, pol);952 T result; // LCOV_EXCL_LINE953 //954 // Trivial case:955 //956 if(s > policies::digits<T, Policy>())957 return 1;958 //959 // Start by seeing if we have a simple closed form:960 //961 if(floor(s) == s)962 {963#ifndef BOOST_MATH_NO_EXCEPTIONS964 // Without exceptions we expect itrunc to return INT_MAX on overflow965 // and we fall through anyway.966 try967 {968#endif969 int v = itrunc(s);970 if(v == s)971 {972 if(v < 0)973 {974 if(((-v) & 1) == 0)975 return 0;976 int n = (-v + 1) / 2;977 if(n <= (int)boost::math::max_bernoulli_b2n<T>::value)978 return T((-v & 1) ? -1 : 1) * boost::math::unchecked_bernoulli_b2n<T>(n) / (1 - v);979 }980 else if((v & 1) == 0)981 {982 if(((v / 2) <= (int)boost::math::max_bernoulli_b2n<T>::value) && (v <= (int)boost::math::max_factorial<T>::value))983 return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * static_cast<T>(pow(constants::pi<T, Policy>(), T(v))) *984 boost::math::unchecked_bernoulli_b2n<T>(v / 2) / boost::math::unchecked_factorial<T>(v);985 return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * static_cast<T>(pow(constants::pi<T, Policy>(), T(v))) *986 boost::math::bernoulli_b2n<T>(v / 2) / boost::math::factorial<T>(v, pol);987 }988 else989 return zeta_imp_odd_integer(v, sc, pol, std::integral_constant<bool, (Tag::value <= 113) && Tag::value>());990 }991#ifndef BOOST_MATH_NO_EXCEPTIONS992 }993 catch(const boost::math::rounding_error&){} // Just fall through, s is too large to round994 catch(const std::overflow_error&){} // LCOV_EXCL_LINE We can only get here for "strange" MP types with small exponents and very large digit counts.995#endif996 }997 998 if(fabs(s) < tools::root_epsilon<T>())999 {1000 result = -0.5f - constants::log_root_two_pi<T, Policy>() * s;1001 }1002 else if(s < 0)1003 {1004 std::swap(s, sc);1005 if(floor(sc/2) == sc/2)1006 result = 0;1007 else1008 {1009 if(s > max_factorial<T>::value)1010 {1011 T mult = boost::math::sin_pi(0.5f * sc, pol) * 2 * zeta_imp(s, sc, pol, tag);1012 result = boost::math::lgamma(s, pol);1013 result -= s * log(2 * constants::pi<T>());1014 if(result > tools::log_max_value<T>())1015 return sign(mult) * policies::raise_overflow_error<T>(function, nullptr, pol);1016 result = exp(result);1017 //1018 // Whether this if branch can be triggered is very type dependent, we need1019 // result to be just on the verge of overflow when /s/ is very close to a1020 // half integer.1021 //1022 if(tools::max_value<T>() / fabs(mult) < result)1023 return boost::math::sign(mult) * policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE1024 result *= mult;1025 }1026 else1027 {1028 result = boost::math::sin_pi(0.5f * sc, pol)1029 * 2 * pow(2 * constants::pi<T>(), -s)1030 * boost::math::tgamma(s, pol)1031 * zeta_imp(s, sc, pol, tag);1032 }1033 }1034 }1035 else1036 {1037 result = zeta_imp_prec(s, sc, pol, tag);1038 }1039 return result;1040}1041 1042} // detail1043 1044template <class T, class Policy>1045inline typename tools::promote_args<T>::type zeta(T s, const Policy&)1046{1047 typedef typename tools::promote_args<T>::type result_type;1048 typedef typename policies::evaluation<result_type, Policy>::type value_type;1049 typedef typename policies::precision<result_type, Policy>::type precision_type;1050 typedef typename policies::normalise<1051 Policy,1052 policies::promote_float<false>,1053 policies::promote_double<false>,1054 policies::discrete_quantile<>,1055 policies::assert_undefined<> >::type forwarding_policy;1056 typedef std::integral_constant<int,1057 precision_type::value <= 0 ? 0 :1058 precision_type::value <= 53 ? 53 :1059 precision_type::value <= 64 ? 64 :1060 precision_type::value <= 113 ? 113 : 01061 > tag_type;1062 1063 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(1064 static_cast<value_type>(s),1065 static_cast<value_type>(1 - static_cast<value_type>(s)),1066 forwarding_policy(),1067 tag_type()), "boost::math::zeta<%1%>(%1%)");1068}1069 1070template <class T>1071inline typename tools::promote_args<T>::type zeta(T s)1072{1073 return zeta(s, policies::policy<>());1074}1075 1076}} // namespaces1077 1078#endif // BOOST_MATH_ZETA_HPP1079 1080 1081 1082