The Riemann zeta function, its zeros and L-functions

The Riemann zeta function is the most studied special function in computational number theory, and acb_dirichlet.h contains the tools to compute with it rigorously: values and derivatives of $\zeta(s)$ and of Dirichlet $L$-functions $L(s, \chi)$ anywhere in the complex plane, and the zeros on the critical line, which are isolated with a proof that none is missed. The algorithms are described in Zeta functions, polylogarithms and L-functions.

Zeros on the critical line#

The nontrivial zeros of $\zeta(s)$ lie in the critical strip $0 < \operatorname{Re}(s) < 1$, and the Riemann hypothesis asserts that they lie on the line $\operatorname{Re}(s) = 1/2$. On this line,

is real for real $t$, where $\theta(t)$ is the Riemann–Siegel theta function, so the zeros $\tfrac12 + i t$ on the line correspond to the real zeros of the Hardy $Z$-function, which can be found from sign changes. Sign changes only show that there are at least so many zeros; to know that all zeros in the strip up to a given height are on the line and have been found, one compares with the number of zeros $N(T)$ with $0 < \operatorname{Im}(s) \le T$ in the whole strip, which is known through the argument principle. In practice this is done with Gram points (where $\theta(t)$ is a multiple of $\pi$, so that sign changes are expected between consecutive ones) and Turing's method [Tur1953], [Leh1970], [Tru2014], which certifies $N(T)$ from the signs of $Z$ at a few Gram points beyond $T$, without evaluating $\zeta$ off the line.

acb_dirichlet_hardy_z_zeros(res, n, len, prec) uses this to compute the zeros number $n, n + 1, \ldots, n + \text{len} - 1$ (counted from 1 by increasing height; the zeros are simple, which is part of what the computation verifies), as balls arb_t of about prec bits of accuracy containing the imaginary parts. It isolates them as just described, evaluating $Z(t)$ with the Riemann–Siegel formula, and refines each one by sign-checked secant and bisection steps. At large heights ($10^{11}$ and beyond) it may instead use Platt's method [Pla2017], which isolates thousands of consecutive zeros at once with a fast multipoint evaluation of $Z$, when a cost model predicts that this is faster.

Values, derivatives and L-functions#

The program below computes three zeros near the bottom of the critical strip and three near the millionth zero, evaluates $Z$ and $Z'$ at the first zero, computes the first Taylor coefficients of $\zeta(s) - 1/(s - 1)$ at $s = 1$, and evaluates the Dirichlet $L$-function of the nontrivial character modulo 4 at $s = 2$, which is Catalan's constant. The two last results are checked against the constants of arb.h.

#include <flint/fmpz.h>
#include <flint/arb.h>
#include <flint/acb.h>
#include <flint/acb_dirichlet.h>

int main(void)
{
    arb_ptr zeros;
    acb_ptr z;
    acb_t s;
    arb_t t;
    fmpz_t n;
    dirichlet_group_t G;
    dirichlet_char_t chi;
    slong i, prec = 128;

    zeros = _arb_vec_init(3);
    z = _acb_vec_init(3);
    acb_init(s);
    arb_init(t);
    fmpz_init(n);

    /* zeros number 1, 2, 3 and 10^6, 10^6 + 1, 10^6 + 2 */
    fmpz_one(n);
    acb_dirichlet_hardy_z_zeros(zeros, n, 3, prec);
    for (i = 0; i < 3; i++)
    {
        flint_printf("rho_%wd = 1/2 + i ", i + 1);
        arb_printn(zeros + i, 30, 0);
        flint_printf("\n");
    }
    fmpz_set_ui(n, 1000000);
    acb_dirichlet_hardy_z_zeros(zeros, n, 3, prec);
    for (i = 0; i < 3; i++)
    {
        flint_printf("rho_%wd = 1/2 + i ", 1000000 + i);
        arb_printn(zeros + i, 30, 0);
        flint_printf("\n");
    }

    /* Z(t) and Z'(t) at the first zero */
    fmpz_one(n);
    acb_dirichlet_hardy_z_zeros(zeros, n, 1, prec);
    acb_set_arb(s, zeros + 0);
    acb_dirichlet_hardy_z(z, s, NULL, NULL, 2, prec);
    flint_printf("Z(t_1)  = "); acb_printn(z + 0, 10, 0); flint_printf("\n");
    flint_printf("Z'(t_1) = "); acb_printn(z + 1, 30, 0); flint_printf("\n");

    /* Taylor coefficients of zeta(s) - 1/(s-1) at s = 1 */
    acb_one(s);
    acb_dirichlet_l_jet(z, s, NULL, NULL, 1, 3, prec);
    for (i = 0; i < 3; i++)
    {
        flint_printf("[(s-1)^%wd] = ", i);
        acb_printn(z + i, 30, 0);
        flint_printf("\n");
    }
    arb_const_euler(t, prec);
    flint_printf("equals Euler's constant: %d\n", arb_overlaps(acb_realref(z), t));

    /* L(2, chi) for the nontrivial character mod 4 is Catalan's constant */
    dirichlet_group_init(G, 4);
    dirichlet_char_init(chi, G);
    dirichlet_char_log(chi, G, 3);
    acb_set_ui(s, 2);
    acb_dirichlet_l(z, s, G, chi, prec);
    flint_printf("L(2, chi_4) = "); acb_printn(z, 30, 0); flint_printf("\n");
    arb_const_catalan(t, prec);
    flint_printf("equals Catalan's constant: %d\n", arb_overlaps(acb_realref(z), t));
    dirichlet_char_clear(chi);
    dirichlet_group_clear(G);

    _arb_vec_clear(zeros, 3);
    _acb_vec_clear(z, 3);
    acb_clear(s);
    arb_clear(t);
    fmpz_clear(n);
    flint_cleanup();
    return 0;
}

rho_1 = 1/2 + i [14.1347251417346937904572519836 +/- 3.76e-29]
rho_2 = 1/2 + i [21.0220396387715549926284795939 +/- 3.10e-30]
rho_3 = 1/2 + i [25.0108575801456887632137909926 +/- 3.72e-29]
rho_1000000 = 1/2 + i [600269.677012444955521233914270 +/- 4.91e-25]
rho_1000001 = 1/2 + i [600270.301090711698660098258485 +/- 1.76e-25]
rho_1000002 = 1/2 + i [600270.747870594366133557973222 +/- 1.38e-25]
Z(t_1)  = [+/- 8.06e-37]
Z'(t_1) = [0.793160433356506116013897565274 +/- 3.53e-31]
[(s-1)^0] = [0.577215664901532860606512090082 +/- 4.03e-31]
[(s-1)^1] = [0.0728158454836767248605863758749 +/- 1.32e-33]
[(s-1)^2] = [-0.00484518159643615924226519301761 +/- 3.74e-33]
equals Euler's constant: 1
L(2, chi_4) = [0.915965594177219015054603514932 +/- 3.85e-31]
equals Catalan's constant: 1

Some points worth noting:

• The zero numbered $10^6$ is computed directly, without finding the zeros below it: the isolation only needs Gram points near the target and the count $N(T)$ from Turing's method.
• $Z(t_1)$ evaluated on the ball for the first zero is a ball around zero, as it must be; its tiny radius reflects the accuracy of the ball for $t_1$.
• acb_dirichlet_hardy_z and acb_dirichlet_l_jet compute len Taylor coefficients, $f(s), f'(s), f''(s)/2!, \ldots$ With the deflate flag, acb_dirichlet_l_jet removes the pole at $s = 1$, so the coefficients above are those of the Laurent expansion $\zeta(s) = \frac{1}{s - 1} + \sum_{k \ge 0} \frac{(-1)^k \gamma_k}{k!} (s - 1)^k$, with the Stieltjes constants $\gamma_0 = \gamma$, $\gamma_1 = -0.0728\ldots$, $\gamma_2 = -0.0096\ldots$ (see also acb_dirichlet_stieltjes, which computes $\gamma_k$ directly, also for large $k$ [JB2018]).
• Passing NULL for the group and the character selects the Riemann zeta function. A Dirichlet character is given by a group dirichlet_group_t for the modulus $q$ and a character dirichlet_char_t; dirichlet_char_log(chi, G, m) selects the character with Conrey label $(q, m)$, here $\chi_4(3) = -1$. See dirichlet.h for characters and acb_dirichlet.h for the $L$-functions.

The example programs#

examples/zeta_zeros.c prints consecutive zeros of $Z(t)$, with their index:

zeta_zeros [-n n] [-count n] [-prec n] [-digits n] [-threads n] [-platt] [-noplatt] [-v] [-verbose] [-h] [-help]

It starts at the zero number -n (default 1) and prints -count zeros (without -count it runs until interrupted). The precision is given in bits with -prec or in digits with -digits; by default it is $64 + \lceil \log_2 n \rceil$ bits, enough to separate zeros at that height comfortably. By default the program lets acb_dirichlet_hardy_z_zeros choose the method; -noplatt forces the Riemann–Siegel method and -platt the large height method (valid from the 10000th zero, with the precision then used as a working precision), which is best for many consecutive zeros at a large height. -threads uses several threads, and -v prints the parameters chosen.

> build/examples/zeta_zeros -n 1 -count 5 -digits 30
1   14.1347251417346937904572519836
2   21.0220396387715549926284795939
3   25.0108575801456887632137909926
4   30.4248761258595132103118975306
5   32.9350615877391896906623689641
cpu/wall(s): 0.008 0.008
virt/peak/rss/peak: 24.8 MB / 24.8 MB / 9.30 MB / 9.30 MB
> build/examples/zeta_zeros -n 1048449114 -count 2 -digits 25
1048449114  388858886.0022851217767971
1048449115  388858886.0023936897027167
cpu/wall(s): 0.114 0.133
virt/peak/rss/peak: 26.0 MB / 26.9 MB / 11.4 MB / 11.9 MB

The second pair is a close pair of zeros (a Lehmer pair) with a distance of about $10^{-4}$, more than a thousand times smaller than the average spacing $2\pi / \log(t / 2\pi) \approx 0.35$ at that height; telling them apart requires evaluating $Z$ near its tiny local extremum between them, which the program does without being told. The printed digits are correct (the digits are printed with ARB_STR_NO_RADIUS, so the last digit may be off by one unit).

examples/lvalue.c evaluates $L(s, \chi)$, or with -z the $Z$-function, and optionally its Taylor series, at a point $s = a + bi$:

lvalue [-character q n] [-re a] [-im b] [-prec p] [-z] [-deflate] [-len l]

The character is given by its Conrey label $(q, n)$ (default $(1, 0)$, the Riemann zeta function), the point by -re and -im (default $s = 1/2$), and the precision in bits (default 53). With -len l it prints $l$ Taylor coefficients, and with -deflate it removes the pole at $s = 1$ as above:

> build/examples/lvalue -character 4 3 -re 2 -prec 128
L(s) = [0.91596559417721901505460351493238411077 +/- 7.86e-39]
> build/examples/lvalue -re 1 -deflate -len 4
L(s) = [0.577215664901532861 +/- 7.54e-19]
L'(s) = [0.072815845483676725 +/- 5.70e-19]
[x^2] L(s+x) = [-0.004845181596436159 +/- 8.03e-19]
[x^3] L(s+x) = [-0.000342305736717224 +/- 6.97e-19]
> build/examples/lvalue -z -re 14.134725 -len 3
Z(s) = [-1.1241835e-7 +/- 3.84e-15]
Z'(s) = [0.7931604148843 +/- 5.24e-14]
[x^2] Z(s+x) = [0.0651645864916 +/- 6.45e-14]

(each followed by timing and memory lines).

examples/lcentral.c computes the central values $L(1/2, \chi)$ for all primitive characters $\chi$ modulo each $q$ in a range, lcentral [--quiet] [--check] [--prec bits] qmin qmax (default 100 bits). It is much faster than calling acb_dirichlet_l for each character: all $L(1/2, \chi)$ modulo $q$ are obtained from the Hurwitz zeta values $\zeta(1/2, k/q)$ by a discrete Fourier transform (acb_dirichlet_l_vec_hurwitz), and the Hurwitz zeta values for all moduli in the range are computed from a common precomputation (acb_dirichlet_hurwitz_precomp_init_num). With --check, the program verifies that no computed central value contains zero (it is conjectured that $L(1/2, \chi) \ne 0$ for all primitive $\chi$); the check can only fail when the precision is too low, in which case the program reports Value could be zero! and aborts.

> build/examples/lcentral 1 5
3,2: [0.48086755769682862618122006324 +/- 7.42e-30] + [+/- 2.45e-30]*I
4,3: [0.66769145718960917665869092930 +/- 1.62e-30] + [+/- 9.34e-31]*I
5,2: [0.76374788011728687822451215264 +/- 2.90e-30] + [0.21696476751886069363858659310 +/- 3.82e-30]*I
5,4: [0.23175094750401575588338366176 +/- 2.12e-30] + [+/- 1.44e-30]*I
5,3: [0.76374788011728687822451215264 +/- 2.90e-30] + [-0.21696476751886069363858659310 +/- 3.82e-30]*I

For real characters the value is real, but the computation through the complex Fourier transform gives a ball for the imaginary part that only contains zero.

examples/keiper_li.c computes the Keiper–Li coefficients $\lambda_0, \ldots, \lambda_n$: with $\xi(s) = \tfrac12 s (s - 1) \pi^{-s/2} \Gamma(s/2) \zeta(s)$, they are defined by $\lambda_n = \frac{1}{(n-1)!} \frac{d^n}{ds^n} \left[s^{n-1} \log \xi(s)\right]_{s=1}$ for $n \ge 1$, and Li's criterion states that the Riemann hypothesis holds if and only if $\lambda_n > 0$ for all $n \ge 1$ [Kei1992], [Li1997]. The program computes the power series of $\log \xi(s)$ at $s = 0$ from those of $\zeta(s)$ (_arb_poly_zeta_series), $\log \Gamma(1 + s/2)$ (_arb_poly_lgamma_series) and $\log(1 - s)$, and obtains the $\lambda_n$ with a binomial transform (_arb_poly_binomial_transform). Since there is heavy cancellation, the default working precision is about $1.1 n + 50$ bits. It was used for the record computations of [Joh2013]:

keiper_li n [-prec prec] [-threads num_threads] [-out out_file]

-out writes all coefficients, as exact midpoints and radii, to a file. For $n = 1000$ the program prints the first and last coefficients with their error bounds, from $\lambda_1 = 0.0230957\ldots$ to $\lambda_{1000} = 2.32605\ldots$.

The zeros can also be found with the general root isolation of arb_calc.h, as in real_roots 0 (see Isolating the real roots of a function), which is slower but shows the method in its simplest form, and counted with the argument principle and numerical integration, as in integrals -i 22 (see Rigorous numerical integration).