INFORMATICA Informatica 1822-8844 0868-4952 0868-4952 Vilnius University INFOR434 10.15388/20-INFOR434 Research Article Series with Binomial-Like Coefficients for Evaluation and 3D Visualization of Zeta Functions BelovasIgorisigoris.belovas@mif.vu.lt∗ SabaliauskasMartynas Vilnius University, Institute of Data Science and Digital Technologies, Akademijos 4, LT-04812 Vilnius, Lithuania Corresponding author. 2020 9112020 314659680 62019 102020 © 2020 Vilnius University2020 Open access article under the CC BY license.

In this paper, we continue the study of efficient algorithms for the computation of zeta functions on the complex plane, extending works of Coffey, Šleževičienė and Vepštas. We prove a central limit theorem for the coefficients of the series with binomial-like coefficients used for evaluation of the Riemann zeta function and establish the rate of convergence to the limiting distribution. An asymptotic expression is derived for the coefficients of the series. We discuss the computational complexity and numerical aspects of the implementation of the algorithm. In the last part of the paper we present our results on 3D visualizations of zeta functions, based on series with binomial-like coefficients. 3D visualizations illustrate underlying structures of surfaces and 3D curves associated with zeta functions.

 zeta functions asymptotic normality limit theorem efficient algorithms 3D visualization Introduction

Zeta functions are very important, playing a pivotal role in the analytic number theory. Properties of zeta functions are essential in the distribution of primes. Prime numbers, in turn, have a broad spectrum of applications, ranging from quantum mechanics to information security (e.g. RSA and the Diffie–Hellman key exchange algorithms or elliptic-curve cryptography). Investigating zeta functions numerically, one needs an algorithm enabling to effectively compute numerous values over the complex plane, not only in the critical strip or on the critical line.

In this paper, we continue the study of efficient algorithms for the computation of zeta functions, extending works of Borwein (2000), Šleževičienė (2004), Vepštas (2008) and Coffey (2009). These algorithms allow the computation of zeta functions over the complex plane. The idea of the method comes from the alternating series convergence and properties of Chebyshev polynomials. The algorithm is nearly optimal in the sense that there is no sequence of n-term exponential polynomials that can converge to a zeta function very much faster than those of the algorithm (see, e.g. Theorem 3.1 in Borwein, 2000, Theorem 5 in Šleževičienė, 2004). The algorithm steps aside for the Riemann–Siegel formula based algorithms (Arias de Reyna, 2011) for computations concerning zeroes on the critical line, where multiple low precision evaluations are required. However, it is more efficient than Euler–Maclaurin based algorithms for arbitrary precision computations (note that the Riemann–Siegel formula permits to approximate ζ(σ+it) in the critical strip for large values of t with a number of terms proportional to t , while Euler–Maclaurin based algorithms require a number of terms proportional to t, cf. Borwein et al., 2000; Borwein, 2000). Belovas obtained central (Belovas, 2019a) and local (Belovas, 2019b) limit theorems, which allowed to introduce an asymptotic approximation for coefficients of the algorithm, providing a considerable speedup in calculations (Belovas and Sakalauskas, 2018).

In the present study, we develop a modification of a globally convergent series for the Riemann zeta function and establish a limit theorem for coefficients of the modified series. We specify the rate of convergence to the limiting distribution, identify the error term and discuss computational complexity. The algorithm is applied to produce 3D visualizations, illustrating underlying structures of surfaces and 3D curves associated with zeta functions.

The paper is organized as follows. The first part is the introduction. In the second part, we introduce a modification of Lerch’s series. In the auxiliary third part, we establish double ordinary and double semi-exponential generating functions for coefficients of the modified series. The proof of the validity of the modification is presented in the fourth section. In the fifth part of the research, we establish analytic expressions of coefficients of the series. In the sixth section, we establish a limit theorem for coefficients and specify the rate of convergence to the limiting distribution. In the seventh section, we prove a theorem for the error term of the series and consider the computational complexity of the method. The numerical aspects of the technique are discussed in the eighth part of the paper. The ninth section of the study is devoted to 3D visualizations of zeta functions and the Riemann hypothesis.

Throughout this paper, we denote by Φ(x) the cumulative distribution function of the standard normal distribution, and by Φ‾(x) we denote the corresponding tail distribution Φ‾(x)=1−Φ(x) . E(X) stands for the expected value of a random variable X. Γ(s) denotes the gamma function. Cnk are the binomial coefficients. ⌊x⌋ and ⌈x⌉ stand for the floor function and the ceiling functions respectively. A×B stands for the Cartesian product of two sets A and B. All limits, unless specified, are taken as n→∞ .

 Series with Binomial-Like Coefficients

Let s=σ+it be a complex variable. The Riemann zeta function is defined by ζ(s)= ∑n=1∞1ns=∏p∈P(1−1ps)−1, for σ>1 1$]]>, and by analytic continuation elsewhere except for a simple pole at s=1 . Let us consider a globally convergent series with binomial-like coefficients for the Riemann zeta function.

For s≠1+2iπm/log2 , m∈Z , ζ(s)=11−21−slimn→∞ ∑k=0n(−1)k(k+1)sqnk, here qnk=ank2n+1. Coefficients ank satisfy a class of triangular arrays and are defined by the following recurrent expression: for n,k∈N , ank=an−1,k−1+an−1,k, where a00=1 , ank=0 if n<k , and an0=2an−1,0+1 if n>0 0$]]> (cf. Table 1).

Numbers ank (3) form a triangular array.

0 1 2 3 4 5 …
0 1 0 0 0 0 0 …
1 3 1 0 0 0 0 …
2 7 4 1 0 0 0 …
3 15 11 5 1 0 0 …
4 31 26 16 6 1 0 …
5 63 57 42 22 7 1 …
… … … … … … … …

It is significant to note that though the partial difference equation (3) bears a close resemblance to the partial difference equation of the binomial coefficients, coefficients ank are not symmetric, because of different boundary conditions. Indeed, the coefficients are decreasing by k for any fixed n. In order to investigate the underlying structure of these numbers, we introduce the double ordinary and the double semi-exponential generating functions (this aspect will be dealt in detail in Section 3).

The series (1) is a modification of the globally convergent series for the Riemann zeta function, first given by Lerch (1897), ζ(s)=11−21−s ∑n=0∞12n+1 ∑k=0n(−1)kCnk(k+1)s. Borwein studied a similar series with binomial-like coefficients, rapidly convergent and well-suited for high precision numerical calculations (Borwein, 2000; Belovas and Sakalauskas, 2018) introduced an improved version of his algorithm. Vepštas (2008) and Coffey (2009) developed efficient algorithms for Hurwitz zeta function, while Šleževičienė (2004) adapted Borwein’s result for Dirichlet L-functions. These algorithms, using series with binomial-like coefficients, are nearly optimal in the sense that there is no sequence of n-term exponential polynomials that can converge to the Riemann zeta function on an interval much faster that those of the algorithms (cf. Theorem 3.1 in Borwein, 2000 and Theorem 5 in Šleževičienė, 2004). In present paper we extend the study of this class of algorithms. We begin by examining the generating functions of the coefficients of the series (1).

 Generating Functions of the Coefficients of the Series

In this section we establish the double ordinary generating function G(x,y)= ∑n=0∞ ∑k=0∞ankxnyk= ∑n=0∞ ∑k=0nankxnyk and the double semi-exponential generating function F(x,y)= ∑n=0∞ ∑k=0∞ankxnn!yk= ∑n=0∞ ∑k=0nankxnn!yk for the coefficients (3). First we will prove an auxiliary lemma for the boundary generating functions.

Coefficients ank have the boundary ordinary generating function G0(x) and the boundary exponential generating function F0(x) as follows: (i)G0(x):=G(x,0)=((1−2x)(1−x))−1, (ii)F0(x):=F(x,0)=2e2x−ex.

Let us consider the formal series (5) of the double ordinary generating function (note that ank=0 for n<k ), G(x,y)= ∑n=0∞ ∑k=0nankxnyk=a00︸=1+ ∑n=1∞an0xn+ ∑n=1∞ ∑k=1nankxnyk. Hence, the boundary ordinary generating function equals G0(x)=G(x,0)=1+ ∑n=1∞an0xn=1+ ∑n=1∞(2an−1,0+1)xn=2 ∑n=1∞an−1,0xn+ ∑n=0∞xn=2x ∑n=0∞an0xn︸=G0(x)+11−x. Thereby, statement (i) of the lemma follows.

Next, take a look at the formal series (6) of the double semi-exponential generating function, F(x,y)= ∑n=0∞ ∑k=0nankxnn!yk=a00︸=1+ ∑n=1∞an0xnn!+ ∑n=1∞ ∑k=1nankxnn!yk. Thus, the boundary exponential generating function equals F0(x)=F(x,0)=1+ ∑n=1∞an0xnn!=1+ ∑n=1∞(2an−1,0+1)xnn!=1+2 ∑n=1∞an−1,0xnn!+ ∑n=1∞xnn!=2 ∑n=0∞an,0xn+1(n+1)!+ex. Differentiating, we receive F0′(x)=2 ∑n=0∞an,0xnn!+ex, and, as a result, the linear differential equation F0′(x)−2F0(x)=ex with the initial condition F0(0)=1 . Indeed, F0(0)=a00+ ∑n=1∞ ∑k=1nankxnn!yk|(0,0)=1. Solving the ordinary differential equation, we obtain statement (ii) of the lemma.  □

Coefficients ank have the double ordinary generating function G(x,y) and the semi-exponential generating function F(x,y) as follows: (i)G(x,y)=((1−2x)(1−x(y+1)))−1, (ii)F(x,y)=(2e2x−(1+y)e(y+1)x)(1−y)−1.

Let us take the double ordinary generating function (5). Substituting the expression (3) into the generating function, we get G(x,y)=1+ ∑n=1∞an0xn+ ∑n=1∞ ∑k=1∞ankxnyk=G0(x)+ ∑n=1∞ ∑k=1∞an−1,k−1xnyk+ ∑n=1∞ ∑k=1∞an−1,kxnyk=G0(x)+xyG(x,y)+x(G(x,y)−G0(x)). Thus, G(x,y)=G0(x)(1−x)((1−x(y+1)))−1. By substituting (7) into (11) we obtain statement (i) of the lemma.

Next, turn to the double semi-exponential generating function (6). Substituting the recurrent expression (3) into the generating function, we have that F(x,y)= ∑n=1∞ ∑k=1nan−1,k−1xnn!yk+ ∑n=1∞ ∑k=1nan−1,kxnn!yk+ ∑n=0∞an0xnn!= ∑n=0∞ ∑k=0nankxn+1(n+1)!yk+1︸=y∫0xFdt+ ∑n=0∞ ∑k=0nankxn+1(n+1)!yk︸=∫0xFdt− ∑n=0∞an0xn+1(n+1)!︸=∫0xF0dt+ ∑n=0∞an0xnn!︸=F0. Therefore, by statement (ii) of Lemma 1, we obtain that F=(y+1)∫0xFdt−∫0xF0dt︸=e2x−ex+F0︸=2e2x−ex=(y+1)∫0xFdt+e2x. Calculating the partial derivative by x, we get the first-order linear partial differential equation, Fx′−(y+1)F=2e2x, with the initial condition F|x=0=1 . Indeed, by (6), we get F(0,y)=a00+ ∑n=1∞ ∑k=1nankxnn!yk|x=0=1. Solving the partial differential equation (note that we can solve it as an ordinary one), we obtain statement (ii) of the lemma.  □

 Analytic Expressions of the Coefficients of the Series

The recurrent formula (3) is often inconvenient in practical computations, thus it is purposeful to seek an analytic expression or an asymptotic expansion. We may view (3) as a partial difference equation with constant coefficients. In this section we will establish analytic expressions for the coefficients ank , using the double ordinary and the semi-exponential generating functions from (see Lemma 2). The coefficients ank can be expressed by linear combinations of binomial coefficients or by the incomplete beta function.

The following analytic expressions are equivalent: (i)ank= ∑j=kn2n−jCjk,(ii)ank= ∑j=k+1n+1Cn+1j,(iii)ank=2n+1I1/2(k+1,n−k+1). Here Ix(a,b) stands for the regularized incomplete beta function: Ix(a,b)=∫0xta−1(1−t)b−1dt∫01ta−1(1−t)b−1dt.

Formal Taylor series in two variables for generating functions G(x,y) and F(x,y) , defined by formulas (5) and (6), are equal to G(x,y)= ∑n=0∞ ∑k=0∞(∂n+k∂xn∂ykG(x,y)(0,0))xnykn!k! and F(x,y)= ∑n=0∞ ∑k=0∞(∂n+k∂xn∂ykF(x,y)(0,0))xnykn!k!, respectively. Statements (i) and (ii) of the lemma are obtained by partial differentiation of the generating functions at (0,0) . By proposition (i) of Lemma 2, we have that ank=1n!k!∂n+k∂xn∂ykG(x,y)(0,0)=1n!k!∂n∂xn11−2xk!xk(1−x(y+1))k+1(0,0)=1n! ∑j=0nCnj(xk)(j)(1(1−2x)(1−x(y+1))k+1)(n−j)(0,0)=1(n−k)! ∑j=0n−kCn−kj(11−2x)(j)(1(1−x(y+1))k+1)(n−k−j)(0,0)= ∑j=0n−k1j!(n−k−j)!j!2j(k+1)(k+2)…(n−j)= ∑j=0n−k2jCn−jk.

Next, using proposition (ii) of Lemma 2, we get ank=1k!∂n+k∂xn∂ykF(x,y)(0,0)=1k!∂k∂yk2n+1e2x−(1+y)n+1e(y+1)x1−y(0,0)=1k! ∑j=0kCkj∂j∂yj2n+1e2x−(1+y)n+1e(y+1)x(0,0)∂k−j∂yk−j11−y(0,0)=− ∑j=1k1j! ∑i=0jCji∂i∂yi(1+y)n+1(0,0)∂j−i∂yj−ie(y+1)x(0,0)+2n+1−1=− ∑j=1k(n+1)!j!(n+1−j)!+2n+1−1=2n+1− ∑j=0kCn+1j= ∑j=k+1n+1Cn+1j.

The third statement of the lemma follows from the formula for the distribution of the probabilities in the Bernoulli scheme (see (19) on page 27 in Bolshev and Smirnov, 1983), ∑k=mnCnkpk(1−p)n−k=Ip(m,n−m+1). Substituting p=1/2 , n=r+1 and m=j+1 , we have that I1/2(j+1,r−j+1)=2−r−1 ∑k=j+1r+1Cr+1k. Thus (cf. (13)), ∑j=k+1n+1Cn+1j︸=ank=2n+1I1/2(k+1,n−k+1), concluding the proof.  □

The first and the second statements of Lemma 3 also can be proved by the direct application of the definition (3) (or by induction). However, these approaches do not reveal the underlying nature of the numbers ank .

Now we can proceed with the proof of Theorem 1.

 Proof of Theorem <xref rid="j_infor434_stat_001">1</xref>

Consider Lerch’s series (4) for the Riemann zeta function ζ(s)=11−21−s ∑n=0∞ ∑k=0n12n+1(−1)kCnk(k+1)s. Changing the order of the summation in the double sum and applying the first statement of Lemma 3 (12), we obtain that ζ(s)=11−21−s ∑k=0∞ ∑n=k∞12n+1(−1)kCnk(k+1)s=11−21−slimN→∞ ∑k=0N(−1)k(k+1)s ∑n=kNCnk2n+1=11−21−slimN→∞ ∑k=0N(−1)k(k+1)s2−N−1 ∑n=kN2N−nCnk︸=aNk yielding us the statement of the theorem.  □

 Limit Theorem for the Coefficients of the Series and the Rate of Convergence to the Limiting Distribution

We will use Hwang’s result on convergence rates in central limit theorems for combinatorial structures (see Theorem 7 and Corollary 2 in Section 4 in Hwang, 1998) to establish the limit theorem for qnk coefficients and specify the rate of convergence to the limiting distribution.

(See Hwang, 1998.) Let Pn(z) be a probability generating function of the random variable Ωn , taking only non-negative integral values, with mean μn and variance σn2 . Suppose that, for each fixed n⩾1 , Pn(z) is a Hurwitz polynomial. If σn→∞ , then Ωn satisfies P(Ωn−μnσn<x)=Φ(x)+O(1σn),x∈R.

Let us formulate a central limit theorem.

Suppose that Fn(x) is the cumulative distribution function of the random variable An with the probability mass function P(An=k)=Cn+1k2n+1−1,k=0,…,n, mean μn=n2(1+1n+O(12n)) and variance σn2=n4(1+1n+O(n2n)). Then it holds Fn(σnx+μn)=Φ(x)+O(1n),x∈R.

Let us consider the moment generating function of the random variable An Mn(s)=E(esAn)=∑kP(An=k)eks= ∑k=0nCn+1keks2n+1−1=12n+1−1( ∑k=0n+1Cn+1keks−es(n+1))=(es+1)n+1−es(n+1)2n+1−1. Calculating the derivatives of Mn(s) , we obtain that Mn′(s)=n+12n+1−1((es+1)nes−es(n+1)),Mn″(s)=n+12n+1−1(n(es+1)n−1e2s+(es+1)nes−(n+1)es(n+1)). Hence, μn=Mn′(0)=(n+1)(2n−1)2n+1−1=n+12(1− ∑k=1∞(12n+1)k)=n+12(1+O(12n))=n2(1+1n+O(12n)) and σn2=Mn″(0)−(Mn′(0))2=n+1(2n+1−1)2(22n−2n−1(n+2))=n+14(1− ∑k=1∞kn+k−12kn+k)=n4(1+1n+O(n2n)), granting us with estimates for the mean and the variance (cf. (16) and (17)). Note that σn→∞ .

Next, let us examine the probability generating function of the random variable An Pn(z)=E(zAn)=∑kP(An=k)zk= ∑k=0nCn+1kzk2n+1−1=(z+1)n+1−zn+12n+1−1. Hurwitz polynomial is a polynomial whose zeros are located in the left half-plane of the complex plane or on the imaginary axis. Let us locate the roots of the polynomial (19). Since z≠0 , we have that (z+1z)n+1=1. Calculating an nth root of unity, we obtain that zk+1zk=cos2πkn+1+isin2πkn+1,k=0,…,n. Thus, the roots of (19) are equal to zk=−12−i2cotπkn+1. The real part of every root is negative. Hence, the probability generating function Pn(z) is a Hurwitz polynomial. By Theorem 2, we have that P(An−μnσn)−Φ(x)=O(1σn)=O(1n), yielding us the statement of the lemma.  □

Now we can formulate the central limit theorem for the coefficients qnk .

Under conditions of Lemma 4, the coefficients qnk satisfy qnk=Φ‾(k−μnσn)+O(1n).

Note that the cumulative distribution function of the random variable An equals Fn(σnx+μn)=∑j⩽σnx+μnCn+1j2n+1−1. Denoting k=⌊σnx+μn⌋ and taking into account Lemma 4 and the equality (cf. (13)) ank2n+1=1− ∑j=0kCn+1j2n+1, we obtain that 2n+12n+1−1 ∑j=0kCn+1j2n+1=2n+12n+1−1(1−ank2n+1)=Φ(k−μnσn)+O(1n). Thus, ank2n+1︸=qnk=1−2n+1−12n+1Φ(k−μnσn)+O(1n)︸=Φ‾(k−μnσn)+O(1n), and the statement of the theorem follows.  □

 Error Term of the Series and the Computational Complexity

Let us investigate the error term of the series (1).

For σ>0 0$]]> and s≠1+2iπm/log2 , m∈Z , ζ(s)=11−21−s ∑k=0n(−1)k(k+1)sqnk+ξn+1(s). Here |ξn+1(s)|⩽(1+|t|/σ)eπ|t|/22n+1|1−21−s|.

By (2.3) of Algorithm 1 in Borwein (2000), the error term of the series equals ξn+1(s)=1pn+1(−1)(1−21−s)1Γ(s)∫01pn+1(x)|logx|s−11+xdx. By the definition of the coefficients of the series and the expression of the generating function of the binomial coefficients, the polynomial pn+1(x) equals pn+1(x)= ∑k=0n+1(−1)kCn+1kxk=(1−x)n+1. Thus, pn+1(−1)=2n+1 and, for 0<x<1 , we have 0<pn+1(x)<1 . It follows |ξn+1(s)|⩽12n+1|1−21−s|1|Γ(s)|∫01(−logx)σ−11+xdx. For the integral factor, we have that ∫01(−logx)σ−11+xdx⩽∫01(−logx)σ−1dx=Γ(σ). The product representation of the gamma function equals |Γ(σ)Γ(s)|2= ∏n=0∞(1+t2(σ+n)2). The product representation of the hyperbolic sine equals sinhπt=πt ∏n=1∞(1+t2n2). Combining (26) and (27), we get |Γ(σ)Γ(s)|2=(1+t2σ2) ∏n=1∞(1+t2(σ+n)2)⩽(1+t2σ2)sinhπ|t|π|t|. Next, combining (24), (25) and (28), we obtain that |ξn+1(s)|⩽12n+1|1−21−s|(1+t2σ2)1/2(sinhπ|t|π|t|)1/2︸⩽(1+|t|/σ)eπ|t|/2, yielding us the statement of the theorem.  □

For σ>0 0$]]> and s≠1+2iπm/log2 , m∈Z , the inequality |ξn+1(s)|⩽12n+1|1−21−s|Γ(σ)|Γ(s)| holds. Hence, to compute the Riemann zeta function with d decimal digits of accuracy, the approach requires a number n of terms in the sum (21), n=⌊logΓ(σ)−log|Γ(s)|+dlog10−log|1−21−s|log2⌋=⌊logΓ(σ)−log|Γ(s)|+dlog10−log1−22−σcos(tlog2)+22−2σlog2⌋⩽⌊logΓ(σ)−log|Γ(σ+it)|+dlog10−log|1−21−σ|log2⌋. Noticing that (cf. 8.328.1, Gradshteyn and Ryzhik, 2014), limt→∞|Γ(σ+it)|t12−σeπ2t=2π, we obtain, for fixed d and σ, that n=πt2log2(1+O(logtt)). For numerical purposes (if t is large enough, e.g. t>0.2 0.2$]]>), we can choose n by the following approximate formula n=⌈πt2log2+(12−σ)logtlog2+Cσ+Cd⌉. Here Cσ=logΓ(σ)−log|1−21−σ|log2,Cd=log10log2d−log2π2log2.

Numerical Aspects

Apart from a certain theoretical interest, the representation (20) of qnk is pivotal in numerical applications of the series (1), providing a considerable speedup in calculations. It relieves the computer memory, circumvents the computation of the binomial coefficients (i.e. time-consuming calculations of factorials). It allows us to cut the tails of the series. We do not to calculate qnk while k→0 and qnk→1 , and do not include corresponding terms into the sum, while k→n and qnk→0 . The process is controlled by choosing the accuracy level ε. For the right tail, we have Φ((k−μn)/σn)>1−ε 1-\varepsilon $]]>. Solving the inequality we obtain that n1=⌈μn+z1−εσn⌉. For the left tail, we have n0=⌊μn−z1−εσn⌋ . Here zp is the quantile function of the standard normal distribution.

The algorithm requires a number of terms n (32) to compute ζ(s) with d decimal digits. However, the use of the asymptotic, proposed in Theorem 3, or analytic expression (iii) of Lemma 3, allows us to select n adaptively, depending on accuracy required and t given. This implies that there is no need to recalculate coefficients with every new s.

Naturally, a direct application of the recurrent equation (3) (if values are stored) is faster, compared to straightforward repeated recalculations. However, this storage-intensive approach is unattractive (note order of n2 storage).

Time required to achieve six decimal digits of accuracy, [sec].

Alg. mod. H=102 and (σ,t)∈Σi×T1 H=103 and (σ,t)∈Σi×T2
Σ0×T1 Σ1×T1 Σ2×T1 Σ0×T2 Σ1×T2 Σ2×T2
RE 0.38 0.36 0.31 293.55 298.09 290.84
BF 0.14 0.14 0.13 12.20 12.47 11.58
NA 0.05 0.05 0.05 4.16 4.41 4.07
RE n1 0.18 0.18 0.16 90.20 97.51 87.74
BF n1 0.10 0.10 0.09 6.63 6.96 6.57
NA n1 0.04 0.04 0.04 2.37 2.53 2.38

In Table 2 we compare the performance of six modifications of the algorithm for the computation of the Riemann zeta function, based on the series with binomial-like coefficients. Namely,

Recurrent Expression (RE). The coefficients of the series (21) are calculated using the recurrent expression (3). The number of terms n in the sum is obtained by applying the approximation (33).

Beta function (BF). The coefficients of the series (21) are calculated using the regularized incomplete beta function (see (iii) of Lemma 3). The number of terms n in the sum is obtained by applying the approximation (33).

Normal approximation (NA). The coefficients of the series (21) are calculated using the asymptotic (20). The standard normal cumulative distribution function is approximated with a 16-digit precision using the fast computation algorithm proposed by (West, 2005). The number of terms n in the sum is obtained by applying the approximation (33).

Recurrent Expression with n1 (RE n1 ). The coefficients of the series (21) are calculated using the recurrent expression (3) and the number of terms n1  (34) in the sum.

Beta function with n1 (BF n1 ). Coefficients of the series (21) are calculated using the regularized incomplete beta function (see (iii) of Lemma 3) and the number of terms n1  (34) in the sum.

Normal approximation with n1 (NA n1 ). Coefficients of the series (21) are calculated using the normal approximation (20) and the number of terms n1  (34) in the sum.

Let us denote the sets for σ: Σ0=(0.7) , Σ1=(1/2,1) , Σ2=(1,10) and for t: T1=(101,102) , T2=(102,103) . We generate uniformly distributed samples of arguments s=(σ,t) of the Riemann zeta function from the sets Σi×Tj . Here × stands for the Cartesian product. Number H in Table 2 stands for the length of a sample. All calculations were performed with i7-2600 CPU, 8 GB RAM (Python 3.6.3).

We can see that asymptotic-based modifications are 2–3 times faster than beta-based modifications, while recurrence-based modifications have proved too computationally demanding for intensive practical use. Note that C++ version of the program would bring more substantial benefits in speedup (cf. Belovas and Sakalauskas, 2018). Numerical experiments have shown that processing 3D visualizations of surfaces and curves, associated with zeta functions, the performance comparable with Zetafast algorithm (Fischer, 2017), while the specified accuracy level does not exceed six decimal digits of accuracy and the coefficients qnk for a mesh are precalculated.

 3D visualizations of Zeta Functios

Visualizing a complex function of a complex variable fC:C→C we have to bear in mind that actually we are faced with a mapping of a 2D subspace to another 2D subspace – an essentially four-dimensional problem fR:R2→R2 . There are a lot of ways of using three dimensions to handle a 4D-structure. We can colour a complex function of a complex variable according to the argument of the function, with the height representing the modulus of the function. We have tried this approach and have found the visualizations unappealing. Our goal is to provide a clear visualization of underlying structures of surfaces and 3D curves associated with zeta functions. Hence, we emphasize on the visualization of 3D intersection lines between the real ℜ(ζ˜) and the imaginary ℑ(ζ˜) surfaces of a zeta function, along with 2D lines of intersection between the real ℜ(ζ˜) and the imaginary ℑ(ζ˜) surfaces of a zeta function with the complex plane. We apply the aforementioned algorithms to obtain illustrations, disclosing the placement of zeroes of a zeta function under consideration, and as a result, enabling us to probe a hypothesis concerning the distribution of zeroes visually. All graphs of zeta functions presented in this work are generated using the normal approximation (NA k0 ) of the coefficients qnk . This algorithm is chosen as the fastest.

The first two figures present two examples of 3D visualizations of the pattern of the allocation of non-trivial zeros of the Riemann zeta function. Figure 1 illustrates the initial allocation (for (σ,t)∈(−1/2,3/2)×(10,35) ). The first five non-trivial zeros are indicated by cyan globes at the points of intersections of the real surface of the Riemann zeta function (yellow surface), the imaginary surface of the Riemann zeta function (red surface) and the critical line σ=1/2 .

Zeroes, ℜζ(s) and ℑζ(s) surfaces for (σ,t)∈(−1/2,3/2)×(10,35) .

Figure 2 visualizes a segment of the pattern, containing 91–100th non-trivial zeroes of the Riemann zeta function (for (σ,t)∈(−1/2,3/2)×(10,35) . The height is cut at the level of |ℜζ(s)|⩽5 and |ℑζ(s)|⩽5 ). Note close pairs of zeroes, resembling Lehmer’s phenomenon (a set of pairs of zeros of the Riemann zeta function, that are close to each other) (Stopple, 2017). It is an unsolved problem, whether there exist infinitely many Lehmer pairs.

Zeroes, ℜζ(s) and ℑζ(s) surfaces for (σ,t)∈(−1/2,3/2)×(10,35) .

ℜζ(s) and ℑζ(s) surfaces for (σ,t)∈(−40,10)×(−20,100) .

Figure 3 shows numerous ridges shaped by the real (yellow) and the imaginary (red) surfaces of the Riemann zeta function (for (σ,t)∈(−40,10)×(−20,100) ). The height is cut at the level of |ℜζ(s)|⩽5 and |ℑζ(s)|⩽5 . Note the peak at the point of a simple pole s=1 . Non-trivial zeros are indicated by small cyan balls. Trivial zeros are indicated by green ones.

Figure 4 extends the visualization of the previous illustration, demonstrating a 2D pattern of intersections of the real and the imaginary surfaces of the Riemann zeta function with the complex plane. Non-trivial zeros are indicated by small cyan balls. Trivial zeros are indicated by green ones.

2D curves corresponding the intersections of ℜζ(s) and ℑζ(s) surfaces with the complex plane, (σ,t)∈(−40,10)×(−20,100) .

Figure 5 shows a pattern of intersections of the real (yellow) and the imaginary (red) surfaces of the Riemann zeta function (for (σ,t)∈(−40,10)×(−20,100) ). The height was cut at the level of |ℜζ(s)|⩽50 and |ℑζ(s)|⩽50 . 3D intersection lines between the real and the imaginary surfaces of the Riemann zeta function are highlighted (bold black).

ℜζ(s) and ℑζ(s) surfaces for (σ,t)∈(−40,10)×(−20,100) .

Figure 6 extends the visualization of the previous illustration, revealing a knot of 3D curves corresponding the intersections of the surfaces, while the surfaces themselves are not shown.

3D curves corresponding the intersections of ℜζ(s) and ℑζ(s) surfaces for (σ,t)∈(−40,10)×(−20,100) .

Figure 7 presents a visualization of the Riemann hypothesis. We can see ribs, formed by the intersections of the real and the imaginary surfaces of the Riemann zeta function, going through the critical line (red) in the critical strip 0<σ<1 (green) at non-trivial zero points (cyan balls).

3D curves corresponding the intersections of ℜζ(s) and ℑζ(s) surfaces in the critical strip.

The next visualizations present a counterexample to a “generalized Riemann hypothesis”. Let us consider the linear combination of the Riemann zeta function and the Dirichlet L-function (cf. Garunkštis and Šimėnas, 2015), f(s,τ)=(1−τ)(1+55s)ζ(s)+τL(s,χ), where τ∈[0,1] and L(s,χ) is the Dirichlet L-function (χ mod 5, χ(2)=−1 ). Hence, L(s,χ)= ∑n=1∞χ(n)ns= ∑n=1∞(1(5n−4)s−1(5n−3)s−1(5n−2)s+1(5n−1)s). By the joint universality theorem for Dirichlet L-functions, it follows that for any 0<τ<1 there are infinitely many zeros of f(s,τ) in the strip 1/2<σ<1 (see Theorem 2 by Kaczorowski and Kulas, 2007).

Let us go with τ=3/4 . The 3D visualization in Fig. 8 shows the pattern of intersections of the real (yellow) and the imaginary (red) surfaces of the linear combination (35) (for (σ,t)∈(−1,2)×(1.65,1.95) ).

ℜf(s,3/4) and ℑf(s,3/4) surfaces for (σ,t)∈(−1,2)×(1.65,1.95) .

Figure 9 demonstrates the intersections of 3D curves (corresponding the intersections of the real and the imaginary surfaces of the linear combination (35)) with the complex plane in the critical strip 0<σ<1 (while the surfaces themselves are not shown). We can see two pairs of intersections in the critical strip (green) not belonging to the critical line (blue).

3D curves corresponding the intersections of ℜf(s,3/4) and ℑf(s,3/4) surfaces in the critical strip.

 References Arias de Reyna, J. (2011). High precision computation of Riemann’s zeta function by the Riemann–Siegel formula. Mathematics of Computation, 80(274), 995–1009. https://doi.org/10.1090/S0025-5718-2010-02426-3. Belovas, I. (2019a). A central limit theorem for coefficients of the modified Borwein method for the calculation of the Riemann zeta-function. Lithuanian Mathematical Journal, 59(1), 17–23. https://doi.org/10.1007/s10986-019-09421-4. Belovas, I. (2019b). A local limit theorem for coefficients of modified Borwein’s method. Glasnik Matematički, 54(74), 1–9. https://doi.org/10.3336/gm.54.1.01. Belovas, I., Sakalauskas, L. (2018). Limit theorems for the coefficients of the modified Borwein method for the calculation of the Riemann zeta-function values. Colloquium Mathematicum, 151(2), 217–227. https://doi.org/10.4064/cm7086-2-2017. Bolshev, L.N., Smirnov, N.V. (1983). Tablitsy matematicheskoj statistiki. Nauka, Moscow (in Russian). Borwein, P. (2000). An efficient algorithm for the Riemann zeta function. In: Constructive, Experimental, and Nonlinear Analysis, Limoges, 1999, CMS Conference Proceedings,Vol. 27. American Mathematical Society, Providence, RI, pp. 29–34. Borwein, J.M., Bradley, D.M., Crandall, R.E. (2000). Computational strategies for the Riemann zeta function. Journal of Computational and Applied Mathematics, 121(1–2), 247–296. https://doi.org/10.1016/S0377-0427(00)00336-8. Coffey, M.W. (2009). An efficient algorithm for the Hurwitz zeta and related functions. Journal of Computational and Applied Mathematics, 225(2), 338–346. https://doi.org/10.1016/j.cam.2008.07.040. Fischer, K. (2017). The Zetafast algorithm for computing zeta functions. arXiv:1703.01414. Garunkštis, R., Šimėnas, R. (2015). On the Speiser equivalent for the Riemann hypothesis. European Journal of Mathematics, 1(2), 337–350. https://doi.org/10.1007/s40879-014-0033-1. Gradshteyn, I.S., Ryzhik, I.M. (2014). Table of Integrals, Series, and Products (8th ed.). Academic Press. Hwang, H.-K. (1998). On convergence rates in the central limit theorems for combinatorial structures. European Journal of Combinatorics, 19(3), 329–343. https://doi.org/10.1006/eujc.1997.0179. Kaczorowski, J., Kulas, M. (2007). On the non-trivial zeros off the critical line for L-functions from the extended Selberg class. Monatshefte für Mathematik, 150(3), 217–232. https://doi.org/10.1007/s00605-006-0412-x. Lerch, M. (1897). Expressions nouvelles de la constante d’Euler, Věstník Královské české společnosti náuk. Tř. mathematicko-přírodovědecká, 42, 1–5. Šleževičienė, R. (2004). An efficient algorithm for computing Dirichlet L-functions. Journal Integral Transforms and Special Functions, 15(6), 513–522. https://doi.org/10.1080/1065246042000272072. Stopple, J. (2017). Lehmer pairs revisited. Experimental Mathematics, 26(1), 45–53. https://doi.org/10.1080/10586458.2015.1107870. Vepštas, L. (2008). An efficient algorithm for accelerating the convergence of oscillatory series, useful for computing the polylogarithm and Hurwitz zeta functions. Numerical Algorithms, 47(3), 211–252. https://doi.org/10.1007/s11075-007-9153-8. West, G. (2005). Better approximations to cumulative normal functions. Wilmott Magazine, 70–76.