71 changes: 71 additions & 0 deletions acb_dirichlet/zeta_rs_d_coeffs.c
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 changes: 70 additions & 0 deletions acb_dirichlet/zeta_rs_f_coeffs.c
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 changes: 215 additions & 0 deletions acb_dirichlet/zeta_rs_r.c
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 changes: 58 additions & 0 deletions doc/source/acb_dirichlet.rst
View file
Expand Up @@ -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
-------------------------------------------------------------------------------

Expand Down
2 changes: 2 additions & 0 deletions doc/source/credits.rst
View file
Expand Up @@ -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
Expand Down