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