hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

implement the Riemann-Siegel formula

Browse source
This commit is contained in:
Fredrik Johansson 2016-10-24 15:36:39 +02:00
parent 1b92a14213
commit 15b9d10818
10 changed files with 774 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
acb_dirichlet.h
View file

@ -33,6 +33,12 @@ void acb_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev,
void acb_dirichlet_powsum_sieved(acb_ptr z, const acb_t s, ulong n, slong len, slong prec);
void acb_dirichlet_powsum_smooth(acb_ptr z, const acb_t s, ulong n, slong len, slong prec);
void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr c, const arb_t p, slong N, slong prec);
void acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong prec);
void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K);
void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec);
void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec);
typedef struct
{
acb_struct s;

84
acb_dirichlet/test/t-zeta_rs.c Normal file
View file

@ -0,0 +1,84 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("zeta_rs....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t z1, z2, s1, s2;
slong prec1, prec2, K1;
acb_init(z1);
acb_init(z2);
acb_init(s1);
acb_init(s2);
if (n_randint(state, 2))
arb_set_d(acb_realref(s1), 0.5);
else
arb_randtest(acb_realref(s1), state, 2 + n_randint(state, 100), 2);
arb_randtest(acb_imagref(s1), state, 2 + n_randint(state, 100), 2);
arb_add_ui(acb_imagref(s1), acb_imagref(s1), n_randtest(state) % 1000, 100);
acb_randtest(z1, state, 2 + n_randint(state, 100), 3);
acb_randtest(z2, state, 2 + n_randint(state, 100), 3);
prec1 = 2 + n_randint(state, 150);
prec2 = 2 + n_randint(state, 150);
K1 = n_randint(state, 10);
if (n_randint(state, 2))
{
mag_zero(arb_radref(acb_realref(s1)));
mag_zero(arb_radref(acb_imagref(s1)));
}
acb_set(s2, s1);
acb_dirichlet_zeta_rs(z1, s1, K1, prec1);
acb_zeta(z2, s2, prec2); /* should be Euler-Maclaurin */
if (!acb_overlaps(z1, z2))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd\n", iter);
flint_printf("K1 = %wd\n\n", K1);
flint_printf("s1 = "); acb_printn(s1, 50, 0); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 50, 0); flint_printf("\n\n");
flint_printf("z1 = "); acb_printn(z1, 50, 0); flint_printf("\n\n");
flint_printf("z2 = "); acb_printn(z2, 50, 0); flint_printf("\n\n");
abort();
}
acb_clear(s1);
acb_clear(s2);
acb_clear(z1);
acb_clear(z2);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

94
acb_dirichlet/test/t-zeta_rs_r.c Normal file
View file

@ -0,0 +1,94 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("zeta_rs_r....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
{
acb_t z1, z2, s1, s2;
slong prec1, prec2, K1, K2;
acb_init(z1);
acb_init(z2);
acb_init(s1);
acb_init(s2);
if (n_randint(state, 2))
arb_set_d(acb_realref(s1), 0.5);
else
arb_randtest(acb_realref(s1), state, 2 + n_randint(state, 100), 2);
arb_randtest(acb_imagref(s1), state, 2 + n_randint(state, 100), 2);
arb_abs(acb_imagref(s1), acb_imagref(s1));
arb_add_ui(acb_imagref(s1), acb_imagref(s1), n_randtest(state) % 1000000, 100);
acb_randtest(z1, state, 2 + n_randint(state, 100), 3);
acb_randtest(z2, state, 2 + n_randint(state, 100), 3);
prec1 = 2 + n_randint(state, 150);
prec2 = 2 + n_randint(state, 150);
K1 = n_randint(state, 10);
K2 = n_randint(state, 10);
if (n_randint(state, 2))
{
mag_zero(arb_radref(acb_realref(s1)));
mag_zero(arb_radref(acb_imagref(s1)));
}
if (n_randint(state, 2))
{
acb_set(s2, s1);
}
else
{
acb_add(s2, s1, z1, prec2);
acb_sub(s2, s2, z1, prec2);
}
acb_dirichlet_zeta_rs_r(z1, s1, K1, prec1);
acb_dirichlet_zeta_rs_r(z2, s2, K2, prec2);
if (!acb_overlaps(z1, z2))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd\n", iter);
flint_printf("K1 = %wd, K2 = %wd\n\n", K1, K2);
flint_printf("s1 = "); acb_printn(s1, 50, 0); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 50, 0); flint_printf("\n\n");
flint_printf("z1 = "); acb_printn(z1, 50, 0); flint_printf("\n\n");
flint_printf("z2 = "); acb_printn(z2, 50, 0); flint_printf("\n\n");
abort();
}
acb_clear(s1);
acb_clear(s2);
acb_clear(z1);
acb_clear(z2);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

84
acb_dirichlet/zeta_rs.c Normal file
View file

@ -0,0 +1,84 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
/* todo: decompose s into midpoint+radius */
void
acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)
{
acb_t R1, R2, X, t;
slong wp;
if (arf_sgn(arb_midref(acb_imagref(s))) < 0)
{
acb_init(t);
acb_conj(t, s);
acb_dirichlet_zeta_rs(res, t, K, prec);
acb_conj(res, res);
acb_clear(t);
return;
}
acb_init(R1);
acb_init(R2);
acb_init(X);
acb_init(t);
/* rs_r increases the precision internally */
wp = prec;
acb_dirichlet_zeta_rs_r(R1, s, K, wp);
if (arb_is_exact(acb_realref(s)) &&
(arf_cmp_2exp_si(arb_midref(acb_realref(s)), -1) == 0))
{
acb_conj(R2, R1);
}
else
{
/* conj(R(conj(1-s))) */
arb_sub_ui(acb_realref(t), acb_realref(s), 1, 10 * wp);
arb_neg(acb_realref(t), acb_realref(t));
arb_set(acb_imagref(t), acb_imagref(s));
acb_dirichlet_zeta_rs_r(R2, t, K, wp);
acb_conj(R2, R2);
}
if (acb_is_finite(R1) && acb_is_finite(R2))
{
wp += 10 + arf_abs_bound_lt_2exp_si(arb_midref(acb_imagref(s)));
wp = FLINT_MAX(wp, 10);
/* X = pi^(s-1/2) gamma((1-s)/2) rgamma(s/2)
= (2 pi)^s rgamma(s) / (2 cos(pi s / 2)) */
acb_rgamma(X, s, wp);
acb_const_pi(t, wp);
acb_mul_2exp_si(t, t, 1);
acb_pow(t, t, s, wp);
acb_mul(X, X, t, wp);
acb_mul_2exp_si(t, s, -1);
acb_cos_pi(t, t, wp);
acb_mul_2exp_si(t, t, 1);
acb_div(X, X, t, wp);
acb_mul(R2, R2, X, wp);
}
/* R1 + X * R2 */
acb_add(res, R1, R2, prec);
acb_clear(R1);
acb_clear(R2);
acb_clear(X);
acb_clear(t);
}

90
acb_dirichlet/zeta_rs_bound.c Normal file
View file

@ -0,0 +1,90 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
/* Arias de Reyna, Theorem 4.2 */
void
acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K)
{
arb_t a;
mag_t c1, c2, c3;
slong e;
if (!arb_is_positive(acb_imagref(s)) || K < 1 || !acb_is_finite(s))
{
mag_inf(err);
return;
}
arb_init(a);
arb_add_ui(a, acb_realref(s), K, MAG_BITS);
arb_sub_ui(a, a, 2, MAG_BITS);
if (!arb_is_nonnegative(acb_realref(s)) && !arb_is_nonnegative(a))
{
mag_inf(err);
arb_clear(a);
return;
}
mag_init(c1);
mag_init(c2);
mag_init(c3);
/* c1 = 1/7 2^(3 sigma / 2) re(sigma) >= 0 */
/* c1 < 1/2 re(sigma) < 0 */
arf_set_mag(arb_midref(a), arb_radref(acb_realref(s)));
arf_add(arb_midref(a), arb_midref(a), arb_midref(acb_realref(s)), MAG_BITS, ARF_RND_CEIL);
if (arf_sgn(arb_midref(a)) <= 0)
{
mag_set_ui_2exp_si(c1, 1, -1);
}
else if (arf_cmp_2exp_si(arb_midref(a), FLINT_BITS - 4) < 0)
{
mag_one(c1);
mag_div_ui(c1, c1, 7);
e = arf_get_si(arb_midref(a), ARF_RND_CEIL);
mag_mul_2exp_si(c1, c1, (3 * e + 1) / 2);
if (mag_cmp_2exp_si(c1, -1) < 0)
mag_set_ui_2exp_si(c1, 1, -1);
}
else
{
mag_inf(c1);
}
/* c2 = 1 / ((10/11) sqrt(t/(2pi))) = (11/10) sqrt((2pi)/t) */
arb_get_mag_lower(c3, acb_imagref(s));
mag_const_pi(c2);
mag_mul_2exp_si(c2, c2, 1);
mag_div(c2, c2, c3);
mag_sqrt(c2, c2);
mag_mul_ui(c2, c2, 11);
mag_div_ui(c2, c2, 10);
/* c2 = c2^(K+1) */
mag_pow_ui(c2, c2, K + 1);
/* c3 = gamma((K+1)/2) */
mag_fac_ui(c3, K / 2);
mag_mul(err, c1, c2);
mag_mul(err, err, c3);
mag_clear(c1);
mag_clear(c2);
mag_clear(c3);
arb_clear(a);
}

71
acb_dirichlet/zeta_rs_d_coeffs.c Normal file
View file

@ -0,0 +1,71 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
void
acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong prec)
{
slong j, r, m;
arb_t u;
arb_init(u);
arb_one(u);
arb_submul_ui(u, sigma, 2, prec);
if (k == 0)
{
arb_one(d);
arb_zero(d + 1);
return;
}
for (j = (3 * k) / 2; j >= 0; j--)
{
m = 3 * k - 2 * j;
if (m != 0)
{
arb_mul_2exp_si(d + j, d + j, -1);
if (j >= 1)
arb_addmul(d + j, d + j - 1, u, prec);
arb_div_ui(d + j, d + j, 2 * m, prec);
if (j >= 2)
arb_submul_ui(d + j, d + j - 2, m + 1, prec);
}
}
if (k % 2 == 0)
{
j = (3 * k) / 2;
arb_zero(d + j);
arb_set_ui(u, 2);
for (r = j - 1; r >= 0; r--)
{
if ((j - r) % 2 == 0)
arb_submul(d + j, d + r, u, prec);
else
arb_addmul(d + j, d + r, u, prec);
arb_mul_ui(u, u, 4 * j - 4 * r + 2, prec);
}
}
arb_zero(d + (3 * k) / 2 + 1);
arb_clear(u);
}

70
acb_dirichlet/zeta_rs_f_coeffs.c Normal file
View file

@ -0,0 +1,70 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
#include "acb_poly.h"
void
acb_dirichlet_zeta_rs_f_coeffs(acb_ptr c, const arb_t p, slong N, slong prec)
{
arb_ptr R, I, T, X;
slong i, len;
R = _arb_vec_init(N);
I = _arb_vec_init(N);
T = _arb_vec_init(N);
X = _arb_vec_init(2);
arb_set(X, p);
arb_one(X + 1);
/* I, R = sin,cos(pi*(X^2/2 + 3/8)) */
len = FLINT_MIN(N, 3);
_arb_poly_mullow(T, X, 2, X, 2, len, prec);
_arb_vec_scalar_mul_2exp_si(T, T, len, -1);
arb_set_d(R, 0.375);
arb_add(T, T, R, prec);
_arb_poly_sin_cos_pi_series(I, R, T, len, N, prec);
/* I -= cos(pi*x/2) * sqrt(2) */
_arb_vec_scalar_mul_2exp_si(X, X, 2, -1);
_arb_poly_cos_pi_series(T, X, 2, N, prec);
arb_sqrt_ui((arb_ptr) c, 2, prec);
_arb_vec_scalar_mul(T, T, N, (arb_ptr) c, prec);
_arb_vec_sub(I, I, T, N, prec);
_arb_vec_scalar_mul_2exp_si(X, X, 2, 1);
/* T = 1 / (2 cos(pi*x)) */
_arb_poly_cos_pi_series(T, X, 2, N, prec);
_arb_vec_scalar_mul_2exp_si(T, T, N, 1);
_arb_poly_inv_series((arb_ptr) c, T, N, N, prec);
_arb_vec_swap(T, (arb_ptr) c, N);
/* R, I *= T */
_arb_poly_mullow((arb_ptr) c, R, N, T, N, N, prec);
_arb_vec_swap(R, (arb_ptr) c, N);
_arb_poly_mullow((arb_ptr) c, I, N, T, N, N, prec);
_arb_vec_swap(I, (arb_ptr) c, N);
for (i = 0; i < N; i++)
{
arb_swap(acb_realref(c + i), R + i);
arb_swap(acb_imagref(c + i), I + i);
}
_acb_poly_inv_borel_transform(c, c, N, prec);
_arb_vec_clear(R, N);
_arb_vec_clear(I, N);
_arb_vec_clear(T, N);
_arb_vec_clear(X, 2);
}

215
acb_dirichlet/zeta_rs_r.c Normal file
View file

@ -0,0 +1,215 @@
/*
Copyright (C) 2016 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_dirichlet.h"
void
acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec)
{
arb_ptr dk;
acb_ptr Fp;
arb_t a, p;
acb_t U, S, u, v;
fmpz_t N;
mag_t err;
slong j, k, wp;
/* determinate K automatically */
if (K <= 0)
{
double sigma, t, log2err, best_log2err;
slong best_K;
sigma = arf_get_d(arb_midref(acb_realref(s)), ARF_RND_DOWN);
t = arf_get_d(arb_midref(acb_imagref(s)), ARF_RND_DOWN);
if (!(sigma > -1e6 && sigma < 1e6) || !(t > 1 && t < 1e40))
{
acb_indeterminate(res);
return;
}
best_K = 1;
best_log2err = 1e300;
/* todo: also break if too slow rate of decay? */
for (K = 1; K < 10 + prec * 0.25; K++)
{
if (sigma < 0 && K + sigma < 3)
continue;
/* Asymptotic approximation of the error term */
log2err = 2.7889996532222537064 - 0.12022458674074695061 / K +
0.2419040680416126037 * K + 0.7213475204444817037 * K * log(K)
+ (-0.7213475204444817037 - 0.7213475204444817037 * K) * log(t);
if (sigma >= 0.0)
log2err += -2.8073549220576041074 + 1.5 * sigma;
if (log2err < best_log2err)
{
best_log2err = log2err;
best_K = K;
}
if (log2err < -prec)
break;
}
K = best_K;
}
mag_init(err);
acb_dirichlet_zeta_rs_bound(err, s, K);
if (!mag_is_finite(err))
{
acb_indeterminate(res);
mag_clear(err);
return;
}
arb_init(a);
arb_init(p);
acb_init(U);
acb_init(S);
acb_init(u);
acb_init(v);
fmpz_init(N);
dk = _arb_vec_init((3 * K) / 2 + 2);
Fp = _acb_vec_init(3 * K + 1);
for (wp = 2 * prec; ; wp *= 2)
{
/* a = sqrt(t / (2pi)) */
arb_const_pi(a, wp);
arb_mul_2exp_si(a, a, 1);
arb_div(a, acb_imagref(s), a, wp);
arb_sqrt(a, a, wp);
/* N = floor(a) */
arb_floor(p, a, wp);
if (!arb_get_unique_fmpz(N, p))
{
if (wp > 4 * prec && wp > arb_rel_accuracy_bits(acb_imagref(s)))
{
acb_indeterminate(res);
goto cleanup;
}
continue;
}
/* p = 1 + 2(N-a) */
arb_sub_fmpz(p, a, N, wp);
arb_neg(p, p);
arb_mul_2exp_si(p, p, 1);
arb_add_ui(p, p, 1, wp);
acb_dirichlet_zeta_rs_f_coeffs(Fp, p, 3 * K + 1, wp);
if (acb_rel_accuracy_bits(Fp + 3 * K) >= prec)
break;
if (wp > 4 * prec && wp > arb_rel_accuracy_bits(acb_imagref(s)))
break;
}
if (!fmpz_fits_si(N))
{
acb_indeterminate(res);
goto cleanup;
}
wp = prec + 10 + 3 * fmpz_bits(N); /* xxx */
wp = FLINT_MAX(wp, prec + 10);
wp = wp + FLINT_BIT_COUNT(K);
acb_zero(S);
for (k = 0; k <= K; k++)
{
acb_dirichlet_zeta_rs_d_coeffs(dk, acb_realref(s), k, wp);
acb_zero(u);
for (j = 0; j <= (3 * k) / 2; j++)
{
/* todo: precompute pi powers */
acb_const_pi(v, wp);
acb_div_onei(v, v);
acb_mul_2exp_si(v, v, -1);
acb_pow_ui(v, v, j, wp);
acb_mul_arb(v, v, dk + j, wp);
acb_addmul(u, v, Fp + 3 * k - 2 * j, wp);
}
acb_const_pi(v, wp);
acb_mul(v, v, v, wp);
acb_mul_arb(v, v, a, wp);
acb_pow_ui(v, v, k, wp);
acb_div(u, u, v, wp);
acb_add(S, S, u, wp);
}
acb_add_error_mag(S, err);
/* U = exp(-i[(t/2) log(t/(2pi)) - t/2 - pi/8]) */
arb_log(acb_realref(u), a, wp);
arb_mul_2exp_si(acb_realref(u), acb_realref(u), 1);
arb_sub_ui(acb_realref(u), acb_realref(u), 1, wp);
arb_mul(acb_realref(u), acb_realref(u), acb_imagref(s), wp);
arb_mul_2exp_si(acb_realref(u), acb_realref(u), -1);
arb_const_pi(acb_imagref(u), wp);
arb_mul_2exp_si(acb_imagref(u), acb_imagref(u), -3);
arb_sub(acb_realref(u), acb_realref(u), acb_imagref(u), wp);
arb_neg(acb_realref(u), acb_realref(u));
arb_sin_cos(acb_imagref(U), acb_realref(U), acb_realref(u), wp);
/* S = (-1)^(N-1) * U * a^(-sigma) * S */
acb_mul(S, S, U, wp);
arb_neg(acb_realref(u), acb_realref(s));
arb_pow(acb_realref(u), a, acb_realref(u), wp);
acb_mul_arb(S, S, acb_realref(u), wp);
if (fmpz_is_even(N))
acb_neg(S, S);
if (_acb_vec_estimate_allocated_bytes(fmpz_get_ui(N) / 6, wp) < 4e9)
acb_dirichlet_powsum_sieved(u, s, fmpz_get_ui(N), 1, wp);
else
acb_dirichlet_powsum_smooth(u, s, fmpz_get_ui(N), 1, wp);
acb_add(S, S, u, wp);
acb_set(res, S); /* don't set_round here; the extra precision is useful */
cleanup:
_arb_vec_clear(dk, (3 * K) / 2 + 2);
_acb_vec_clear(Fp, 3 * K + 1);
arb_clear(a);
arb_clear(p);
acb_clear(U);
acb_clear(S);
acb_clear(u);
acb_clear(v);
fmpz_clear(N);
mag_clear(err);
}

58
doc/source/acb_dirichlet.rst
View file

@ -63,6 +63,64 @@ Truncated L-series and power sums
A slightly bigger gain for larger *n* could be achieved by using more
small prime factors, at the expense of space.
Riemann zeta function and Riemann-Siegel formula
-------------------------------------------------------------------------------
The Riemann-Siegel (RS) formula is implemented closely following
J. Arias de Reyna [Ari2011]_.
For `s = \sigma + it` with `t > 0`, the expansion takes the form
.. math ::
\zeta(s) = \mathcal{R}(s) + X(s) \mathcal{R}(1-s), \quad X(s) = \pi^{s-1/2} \frac{\Gamma((1-s)/2)}{\Gamma(s/2)}
where
.. math ::
\mathcal{R}(s) = \sum_{k=1}^N \frac{1}{k^s} + (-1)^{N-1} U a^{-\sigma} \left[ \sum_{k=0}^K \frac{C_k(p)}{a^k} + RS_K \right]
.. math ::
U = \exp\left(-i\left[ \frac{t}{2} \log\left(\frac{t}{2\pi}\right)-\frac{t}{2}-\frac{\pi}{8} \right]\right), \quad
a = \sqrt{\frac{t}{2\pi}}, \quad N = \lfloor a \rfloor, \quad p = 1-2(a-N).
The coefficients `C_k(p)` in the asymptotic part of the expansion
are expressed in terms of certain auxiliary coefficients `d_j^{(k)}`
and `F^{(j)}(p)`.
Because of artificial discontinuities, *s* should be exact inside
the evaluation (automatic reduction to the exact case is not yet implemented).
.. function:: void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr f, const arb_t p, slong n, slong prec)
Computes the coefficients `F^{(j)}(p)` for `0 \le j < n`.
Uses power series division. This method breaks down when `p = \pm 1/2`
(which is not problem if *s* is an exact floating-point number).
.. function:: void acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong prec)
Computes the coefficients `d_j^{(k)}` for `0 \le j \le \lfloor 3k/2 \rfloor + 1`.
On input, the array *d* must contain the coefficients for `d_j^{(k-1)}`
unless `k = 0`, and these coefficients will be updated in-place.
.. function:: void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K)
Bounds the error term `RS_K` following Theorem 4.2 in Arias de Reyna.
.. function:: void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec)
Computes `\mathcal{R}(s)` in the upper half plane. Uses precisely *K*
asymptotic terms in the RS formula if this input parameter is positive;
otherwise chooses the number of terms automatically based on *s* and the
precision.
.. function:: void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)
Computes `\zeta(s)` using the Riemann-Siegel formula. Uses precisely
*K* asymptotic terms in the RS formula if this input parameter is positive;
otherwise chooses the number of terms automatically based on *s* and the
precision.
Hurwitz zeta function
-------------------------------------------------------------------------------

2
doc/source/credits.rst
View file

@ -123,6 +123,8 @@ Bibliography
(In the PDF edition, this section is empty. See the bibliography listing at the end of the document.)
.. [Ari2011] \J. Arias de Reyna, "High precision computation of Riemann’s zeta function by the Riemann-Siegel formula, I", Mathematics of Computation 80 (2011), 995-1009
.. [Arn2010] \J. Arndt, *Matters Computational*, Springer (2010), http://www.jjj.de/fxt/#fxtbook
.. [BBC1997] \D. H. Bailey, J. M. Borwein and R. E. Crandall, "On the Khintchine constant", Mathematics of Computation 66 (1997) 417-431