74 -.5 - 0.03532739323390276872, 0.3442898999246284869,
75 0.03597993651536150163, 0.00126461541144692592,
76 0.00002286212103119451, 0.00000025347910790261,
77 0.00000000190451637722, 0.00000000001034969525,
78 0.00000000000004259816, 0.00000000000000013744,
79 0.00000000000000000035};
91 2.5 - 0.07643947903327941, -0.02235652605699819, 0.00077341811546938,
92 -0.00004281006688886, 0.00000308170017386, -0.00000026393672220,
93 0.00000002563713036, -0.00000000274270554, 0.00000000031694296,
94 -0.00000000003902353, 0.00000000000506804, -0.00000000000068895,
95 0.00000000000009744, -0.00000000000001427, 0.00000000000000215,
96 -0.00000000000000033, 0.00000000000000005};
108 2.5 - 0.01201869826307592, -0.00917485269102569, 0.00014445509317750,
109 -0.00000401361417543, 0.00000015678318108, -0.00000000777011043,
110 0.00000000046111825, -0.00000000003158592, 0.00000000000243501,
111 -0.00000000000020743, 0.00000000000001925, -0.00000000000000192,
112 0.00000000000000020, -0.00000000000000002};
127 5.5 - .07660547252839144951, 1.92733795399380827000, .22826445869203013390,
128 .01304891466707290428, .00043442709008164874, .00000942265768600193,
129 .00000014340062895106, .00000000161384906966, .00000000001396650044,
130 .00000000000009579451, .00000000000000053339, .00000000000000000245};
145 1.5 + 0.0253002273389477705, -0.3531559607765448760, -0.1226111808226571480,
146 -0.0069757238596398643, -0.0001730288957513052, -0.0000024334061415659,
147 -0.0000000221338763073, -0.0000000001411488392, -0.0000000000006666901,
148 -0.0000000000000024274, -0.0000000000000000070};
160 2.5 + 0.27443134069738830, 0.07571989953199368, -0.00144105155647540,
161 0.00006650116955125, -0.00000436998470952, 0.00000035402774997,
162 -0.00000003311163779, 0.00000000344597758, -0.00000000038989323,
163 0.00000000004720819, -0.00000000000604783, 0.00000000000081284,
164 -0.00000000000011386, 0.00000000000001654, -0.00000000000000248,
165 0.00000000000000038, -0.00000000000000006};
177 2.5 + 0.06379308343739001, 0.02832887813049721, -0.00024753706739052,
178 0.00000577197245160, -0.00000020689392195, 0.00000000973998344,
179 -0.00000000055853361, 0.00000000003732996, -0.00000000000282505,
180 0.00000000000023720, -0.00000000000002176, 0.00000000000000215,
181 -0.00000000000000022, 0.00000000000000002};
196 1.75 - 0.001971713261099859, 0.407348876675464810, 0.034838994299959456,
197 0.001545394556300123, 0.000041888521098377, 0.000000764902676483,
198 0.000000010042493924, 0.000000000099322077, 0.000000000000766380,
199 0.000000000000004741, 0.000000000000000024};
206 1.00000000000000000000000000000, 0.083333333333333333333333333333,
207 -0.00138888888888888888888888888889, 0.000033068783068783068783068783069,
208 -8.2671957671957671957671957672e-07, 2.0876756987868098979210090321e-08,
209 -5.2841901386874931848476822022e-10, 1.3382536530684678832826980975e-11,
210 -3.3896802963225828668301953912e-13, 8.5860620562778445641359054504e-15,
211 -2.1748686985580618730415164239e-16, 5.5090028283602295152026526089e-18,
212 -1.3954464685812523340707686264e-19, 3.5347070396294674716932299778e-21,
213 -8.9535174270375468504026113181e-23};
216 constexpr auto max_bits = 54.0;
217 constexpr auto jmax = 12;
218 constexpr auto kmax = 10;
220 if ((s > max_bits and q < 1.0) or (s > 0.5 * max_bits and q < 0.25))
221 return std::pow(q, -s);
222 if (s > 0.5 * max_bits and q < 1.0) {
223 auto const p1 = std::pow(q, -s);
224 auto const p2 = std::pow(q / (1.0 + q), s);
225 auto const p3 = std::pow(q / (2.0 + q), s);
226 return p1 * (1.0 + p2 + p3);
231 auto const kmax_q =
static_cast<double>(kmax) + q;
232 auto const pmax = std::pow(kmax_q, -s);
234 auto pcp = pmax / kmax_q;
235 auto ans = pmax * (kmax_q / (s - 1.0) + 0.5);
237 for (
int k = 0; k < kmax; k++)
238 ans += std::pow(
static_cast<double>(k) + q, -s);
240 for (
int j = 0; j <= jmax; j++) {
241 auto const delta =
hzeta_c[j + 1] * scp * pcp;
243 scp *= (s + 2. * j + 1.) * (s + 2. * j + 2.);
244 pcp /=
Utils::sqr(
static_cast<double>(kmax) + q);
259 return std::exp(-x) * c / std::sqrt(x);
271 return std::exp(-x) * c / std::sqrt(x);
284 11, 11, 10, 10, 9, 9,
286 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 2, 2};
290 auto const tmp = .5 * std::exp(-x) / std::sqrt(x);
294 auto const tmp = std::exp(-x) / std::sqrt(x);
295 auto const xx = (16. / 3.) / x - 5. / 3.;
304 x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.;
307 x2 = (2. * 16.) / x - 2.;
310 auto d0 = x2 * dd0 + s0[j - 1];
311 for (j -= 2; j >= 1; j--) {
312 auto const tmp0 = d0;
313 d0 = x2 * d0 - dd0 + s0[j];
316 auto const tmp = std::exp(-x) / std::sqrt(x);
317 return tmp * (0.5 * (s0[0] + x2 * d0) - dd0);
323 auto x2 = (2. / 4.5) * x * x - 2.;
325 auto d0 = x2 * dd0 +
bi0_cs[j - 1];
326 for (j -= 2; j >= 1; j--) {
327 auto const tmp0 = d0;
328 d0 = x2 * d0 - dd0 +
bi0_cs[j];
332 auto const ret = -tmp * (0.5 * (
bi0_cs[0] + x2 * d0) - dd0);
338 d0 = x2 * dd0 +
bk0_cs[j - 1];
339 for (j -= 2; j >= 1; j--) {
340 auto const tmp0 = d0;
341 d0 = x2 * d0 - dd0 +
bk0_cs[j];
344 return ret + (0.5 * (x2 * d0 +
bk0_cs[0]) - dd0);
350 auto const tmp = .5 * std::exp(-x) / std::sqrt(x);
354 auto const tmp = std::exp(-x) / std::sqrt(x);
355 auto const xx = (16. / 3.) / x - 5. / 3.;
364 x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.;
367 x2 = (2. * 16.) / x - 2.;
370 auto d1 = x2 * dd1 + s1[j - 1];
371 for (j -= 2; j >= 1; j--) {
372 auto const tmp1 = d1;
373 d1 = x2 * d1 - dd1 + s1[j];
376 auto const tmp = std::exp(-x) / std::sqrt(x);
377 return tmp * (0.5 * (s1[0] + x2 * d1) - dd1);
383 auto x2 = (2. / 4.5) * x * x - 2.;
385 auto d1 = x2 * dd1 +
bi1_cs[j - 1];
386 for (j -= 2; j >= 1; j--) {
387 auto const tmp1 = d1;
388 d1 = x2 * d1 - dd1 +
bi1_cs[j];
392 auto const ret = x * tmp * (0.5 * (
bi1_cs[0] + x2 * d1) - dd1);
398 d1 = x2 * dd1 +
bk1_cs[j - 1];
399 for (j -= 2; j >= 1; j--) {
400 auto const tmp1 = d1;
401 d1 = x2 * d1 - dd1 +
bk1_cs[j];
404 return ret + (0.5 * (x2 * d1 +
bk1_cs[0]) - dd1) / x;
408std::pair<double, double>
LPK01(
double x) {
410 auto const tmp = .5 * std::exp(-x) / std::sqrt(x);
411 auto const k0 = tmp *
ak0_cs[0];
412 auto const k1 = tmp *
ak1_cs[0];
416 auto const tmp = std::exp(-x) / std::sqrt(x);
417 auto const xx = (16. / 3.) / x - 5. / 3.;
418 auto const k0 = tmp * (xx *
ak0_cs[1] + 0.5 *
ak0_cs[0]);
419 auto const k1 = tmp * (xx *
ak1_cs[1] + 0.5 *
ak1_cs[0]);
429 x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.;
433 x2 = (2. * 16.) / x - 2.;
437 auto d0 = x2 * dd0 + s0[j - 1];
438 auto d1 = x2 * dd1 + s1[j - 1];
439 for (j -= 2; j >= 1; j--) {
440 auto const tmp0 = d0, tmp1 = d1;
441 d0 = x2 * d0 - dd0 + s0[j];
442 d1 = x2 * d1 - dd1 + s1[j];
446 auto const tmp = std::exp(-x) / std::sqrt(x);
447 auto const k0 = tmp * (0.5 * (s0[0] + x2 * d0) - dd0);
448 auto const k1 = tmp * (0.5 * (s1[0] + x2 * d1) - dd1);
455 auto x2 = (2. / 4.5) * x * x - 2.;
458 auto d0 = x2 * dd0 +
bi0_cs[j - 1];
459 auto d1 = x2 * dd1 +
bi1_cs[j - 1];
460 for (j -= 2; j >= 1; j--) {
461 auto const tmp0 = d0, tmp1 = d1;
462 d0 = x2 * d0 - dd0 +
bi0_cs[j];
463 d1 = x2 * d1 - dd1 +
bi1_cs[j];
468 auto k0 = -tmp * (0.5 * (
bi0_cs[0] + x2 * d0) - dd0);
469 auto k1 = x * tmp * (0.5 * (
bi1_cs[0] + x2 * d1) - dd1);
476 d0 = x2 * dd0 +
bk0_cs[j - 1];
477 d1 = x2 * dd1 +
bk1_cs[j - 1];
478 for (j -= 2; j >= 1; j--) {
479 auto const tmp0 = d0, tmp1 = d1;
480 d0 = x2 * d0 - dd0 +
bk0_cs[j];
481 d1 = x2 * d1 - dd1 +
bk1_cs[j];
485 k0 += (0.5 * (x2 * d0 +
bk0_cs[0]) - dd0);
486 k1 += (0.5 * (x2 * d1 +
bk1_cs[0]) - dd1) / x;
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
DEVICE_QUALIFIER constexpr T ln_2()
Natural logarithm of 2.
double hzeta(double s, double q)
Hurwitz zeta function.
static int ak01_orders[]
necessary orders for K0/1 from 2 up to 22 for 10^-14 precision.
static double ak02_cs[14]
Series for ak02.
double LPK1(double x)
Modified Bessel function of second kind, order 1, low precision.
static double bi0_cs[12]
Series for bi0.
static double bk1_cs[11]
Series for bk1.
std::pair< double, double > LPK01(double x)
Modified Bessel functions of second kind, order 0 and 1, low precision.
double K1(double x)
Modified Bessel function of second kind, order 1.
static double ak12_cs[14]
Series for ak12.
static double ak0_cs[17]
Series for ak0.
static double ak1_cs[17]
Series for ak1.
static double bk0_cs[11]
Series for bk0.
static double bi1_cs[11]
Series for bi1.
double K0(double x)
Modified Bessel function of second kind, order 0.
static double const hzeta_c[15]
Coefficients for Maclaurin summation in hzeta().
double LPK0(double x)
Modified Bessel function of second kind, order 0, low precision.
__device__ float evaluateAsChebychevSeriesAt(float const *c, int n, float x)
This file contains implementations for some special functions which are needed by the MMM family of a...