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 $\zeta (\sigma +it)$ in the critical strip for large values of t with a number of terms proportional to $\sqrt{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 $\Phi (x)$ the cumulative distribution function of the standard normal distribution, and by $\overline{\Phi }(x)$ we denote the corresponding tail distribution $\overline{\Phi }(x)=1-\Phi (x)$. $E(X)$ stands for the expected value of a random variable X. $\Gamma (s)$ denotes the gamma function. ${C_{n}^{k}}$ are the binomial coefficients. $\lfloor x\rfloor $ and $\lceil x\rceil $ stand for the floor function and the ceiling functions respectively. $A\times B$ stands for the Cartesian product of two sets A and B. All limits, unless specified, are taken as $n\to \infty $.

Let $s=\sigma +it$ be a complex variable. The Riemann zeta function is defined by for $\sigma >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.

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 ${a_{nk}}$ 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), 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).

In this section we establish the double ordinary generating function and the double semi-exponential generating function for the coefficients (3). First we will prove an auxiliary lemma for the boundary generating functions.

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 ${a_{nk}}$, using the double ordinary and the semi-exponential generating functions from (see Lemma 2). The coefficients ${a_{nk}}$ can be expressed by linear combinations of binomial coefficients or by the incomplete beta function.

Now we can proceed with the proof of Theorem 1.

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 ${q_{nk}}$ coefficients and specify the rate of convergence to the limiting distribution.

Let us formulate a central limit theorem.

Now we can formulate the central limit theorem for the coefficients ${q_{nk}}$.

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

Apart from a certain theoretical interest, the representation (20) of ${q_{nk}}$ 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 ${q_{nk}}$ while $k\to 0$ and ${q_{nk}}\to 1$, and do not include corresponding terms into the sum, while $k\to n$ and ${q_{nk}}\to 0$. The process is controlled by choosing the accuracy level ε. For the right tail, we have $\Phi ((k-{\mu _{n}})/{\sigma _{n}})>1-\varepsilon $. Solving the inequality we obtain that For the left tail, we have ${n_{0}}=\lfloor {\mu _{n}}-{z_{1-\varepsilon }}{\sigma _{n}}\rfloor $. Here ${z_{p}}$ is the quantile function of the standard normal distribution.

The algorithm requires a number of terms n (32) to compute $\zeta (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 ${n^{2}}$ storage).

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,

Let us denote the sets for σ: ${\Sigma _{0}}=(0.7)$, ${\Sigma _{1}}=(1/2,1)$, ${\Sigma _{2}}=(1,10)$ and for t: ${T_{1}}=({10^{1}},{10^{2}})$, ${T_{2}}=({10^{2}},{10^{3}})$. We generate uniformly distributed samples of arguments $s=(\sigma ,t)$ of the Riemann zeta function from the sets ${\Sigma _{i}}\times {T_{j}}$. 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 ${q_{nk}}$ for a mesh are precalculated.

Visualizing a complex function of a complex variable ${f_{\mathbb{C}}}:\mathbb{C}\to \mathbb{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 ${f_{\mathbb{R}}}:{\mathbb{R}^{2}}\to {\mathbb{R}^{2}}$. 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 $\mathrm{\Re }(\tilde{\zeta })$ and the imaginary $\mathrm{\Im }(\tilde{\zeta })$ surfaces of a zeta function, along with 2D lines of intersection between the real $\mathrm{\Re }(\tilde{\zeta })$ and the imaginary $\mathrm{\Im }(\tilde{\zeta })$ 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${k_{0}}$) of the coefficients ${q_{nk}}$. 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 $(\sigma ,t)\in (-1/2,3/2)\times (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 $\sigma =1/2$.

Figure 2 visualizes a segment of the pattern, containing 91–100th non-trivial zeroes of the Riemann zeta function (for $(\sigma ,t)\in (-1/2,3/2)\times (10,35)$. The height is cut at the level of $|\mathrm{\Re }\zeta (s)|\leqslant 5$ and $|\mathrm{\Im }\zeta (s)|\leqslant 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.

Figure 3 shows numerous ridges shaped by the real (yellow) and the imaginary (red) surfaces of the Riemann zeta function (for $(\sigma ,t)\in (-40,10)\times (-20,100)$). The height is cut at the level of $|\mathrm{\Re }\zeta (s)|\leqslant 5$ and $|\mathrm{\Im }\zeta (s)|\leqslant 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.

Figure 5 shows a pattern of intersections of the real (yellow) and the imaginary (red) surfaces of the Riemann zeta function (for $(\sigma ,t)\in (-40,10)\times (-20,100)$). The height was cut at the level of $|\mathrm{\Re }\zeta (s)|\leqslant 50$ and $|\mathrm{\Im }\zeta (s)|\leqslant 50$. 3D intersection lines between the real and the imaginary surfaces of the Riemann zeta function are highlighted (bold black).

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.

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<\sigma <1$ (green) at non-trivial zero points (cyan balls).

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), where $\tau \in [0,1]$ and $L(s,\chi )$ is the Dirichlet L-function (χ mod 5, $\chi (2)=-1$). Hence, By the joint universality theorem for Dirichlet L-functions, it follows that for any $0<\tau <1$ there are infinitely many zeros of $f(s,\tau )$ in the strip $1/2<\sigma <1$ (see Theorem 2 by Kaczorowski and Kulas, 2007).

Let us go with $\tau =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 $(\sigma ,t)\in (-1,2)\times (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<\sigma <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).