Series with Binomial-Like Coefficients for Evaluation and 3D Visualization of Zeta Functions

Abstract

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.

Keywords

zeta functions asymptotic normality limit theorem efficient algorithms 3D visualization

1 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 $ζ (σ + i t)$ 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 $Φ (x)$ the cumulative distribution function of the standard normal distribution, and by $\overline{Φ} (x)$ we denote the corresponding tail distribution $\overline{Φ} (x) = 1 - Φ (x)$ . $E (X)$ stands for the expected value of a random variable X. $Γ (s)$ denotes the gamma function. $C_{n}^{k}$ are the binomial coefficients. $⌊ x ⌋$ and $⌈ x ⌉$ 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$ .

2 Series with Binomial-Like Coefficients

Let $s = σ + i t$ be a complex variable. The Riemann zeta function is defined by $\begin{matrix} ζ (s) = \sum_{n = 1}^{\infty} \frac{1}{n^{s}} = \prod_{p \in P} {(1 - \frac{1}{p^{s}})}^{- 1}, \end{matrix}$ for $σ > 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.

Theorem 1.
For $s \neq 1 + 2 i π m / log 2$ , $m \in Z$ , $\begin{matrix} ζ (s) = \frac{1}{1 - 2^{1 - s}} lim_{n \to \infty} \sum_{k = 0}^{n} \frac{{(- 1)}^{k}}{{(k + 1)}^{s}} q_{n k}, \end{matrix}$ (1) here $\begin{matrix} q_{n k} = \frac{a_{n k}}{2^{n + 1}} . \end{matrix}$ (2) Coefficients $a_{n k}$ satisfy a class of triangular arrays and are defined by the following recurrent expression: for $n, k \in N$ , $\begin{matrix} a_{n k} = a_{n - 1, k - 1} + a_{n - 1, k}, \end{matrix}$ (3) where $a_{00} = 1$ , $a_{n k} = 0$ if $n < k$ , and $a_{n 0} = 2 a_{n - 1, 0} + 1$ if $n > 0$ (cf. Table 1).

Table 1
Numbers $a_{n k}$ (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 $a_{n k}$ 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), $\begin{matrix} ζ (s) = \frac{1}{1 - 2^{1 - s}} \sum_{n = 0}^{\infty} \frac{1}{2^{n + 1}} \sum_{k = 0}^{n} \frac{{(- 1)}^{k} C_{n}^{k}}{{(k + 1)}^{s}} . \end{matrix}$ (4) 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).
3 Generating Functions of the Coefficients of the Series

	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	…
…	…	…	…	…	…	…	…

In this section we establish the double ordinary generating function $\begin{matrix} G (x, y) = \sum_{n = 0}^{\infty} \sum_{k = 0}^{\infty} a_{n k} x^{n} y^{k} = \sum_{n = 0}^{\infty} \sum_{k = 0}^{n} a_{n k} x^{n} y^{k} \end{matrix}$ (5) and the double semi-exponential generating function $\begin{matrix} F (x, y) = \sum_{n = 0}^{\infty} \sum_{k = 0}^{\infty} a_{n k} \frac{x^{n}}{n!} y^{k} = \sum_{n = 0}^{\infty} \sum_{k = 0}^{n} a_{n k} \frac{x^{n}}{n!} y^{k} \end{matrix}$ (6) for the coefficients (3). First we will prove an auxiliary lemma for the boundary generating functions. Lemma 1.

Coefficients $a_{n k}$ have the boundary ordinary generating function $G_{0} (x)$ and the boundary exponential generating function $F_{0} (x)$ as follows: $\begin{aligned} (i) G_{0} (x) : = & G (x, 0) = {((1 - 2 x) (1 - x))}^{- 1}, \end{aligned}$ (7) $\begin{aligned} (ii) F_{0} (x) : = & F (x, 0) = 2 e^{2 x} - e^{x} . \end{aligned}$ (8)

Proof.

Let us consider the formal series (5) of the double ordinary generating function (note that $a_{n k} = 0$ for $n < k$ ), $\begin{matrix} G (x, y) = \sum_{n = 0}^{\infty} \sum_{k = 0}^{n} a_{n k} x^{n} y^{k} = \underset{= 1}{\underset{︸}{a_{00}}} + \sum_{n = 1}^{\infty} a_{n 0} x^{n} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{n} a_{n k} x^{n} y^{k} . \end{matrix}$ Hence, the boundary ordinary generating function equals $\begin{aligned} G_{0} (x) & = G (x, 0) = 1 + \sum_{n = 1}^{\infty} a_{n 0} x^{n} = 1 + \sum_{n = 1}^{\infty} (2 a_{n - 1, 0} + 1) x^{n} \\ = 2 \sum_{n = 1}^{\infty} a_{n - 1, 0} x^{n} + \sum_{n = 0}^{\infty} x^{n} = 2 x \underset{= G_{0} (x)}{\underset{︸}{\sum_{n = 0}^{\infty} a_{n 0} x^{n}}} + \frac{1}{1 - x} . \end{aligned}$ Thereby, statement (i) of the lemma follows.

Next, take a look at the formal series (6) of the double semi-exponential generating function, $\begin{matrix} F (x, y) = \sum_{n = 0}^{\infty} \sum_{k = 0}^{n} a_{n k} \frac{x^{n}}{n!} y^{k} = \underset{= 1}{\underset{︸}{a_{00}}} + \sum_{n = 1}^{\infty} a_{n 0} \frac{x^{n}}{n!} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{n} a_{n k} \frac{x^{n}}{n!} y^{k} . \end{matrix}$ Thus, the boundary exponential generating function equals $\begin{aligned} F_{0} (x) & = F (x, 0) = 1 + \sum_{n = 1}^{\infty} a_{n 0} \frac{x^{n}}{n!} = 1 + \sum_{n = 1}^{\infty} (2 a_{n - 1, 0} + 1) \frac{x^{n}}{n!} \\ = 1 + 2 \sum_{n = 1}^{\infty} a_{n - 1, 0} \frac{x^{n}}{n!} + \sum_{n = 1}^{\infty} \frac{x^{n}}{n!} = 2 \sum_{n = 0}^{\infty} a_{n, 0} \frac{x^{n + 1}}{(n + 1)!} + e^{x} . \end{aligned}$ Differentiating, we receive $\begin{matrix} F_{0}^{'} (x) = 2 \sum_{n = 0}^{\infty} a_{n, 0} \frac{x^{n}}{n!} + e^{x}, \end{matrix}$ and, as a result, the linear differential equation $\begin{matrix} F_{0}^{'} (x) - 2 F_{0} (x) = e^{x} \end{matrix}$ with the initial condition $F_{0} (0) = 1$ . Indeed, $\begin{matrix} F_{0} (0) = a_{00} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{n} a_{n k} \frac{x^{n}}{n!} y^{k} |_{(0, 0)} = 1 . \end{matrix}$ Solving the ordinary differential equation, we obtain statement (ii) of the lemma. □

Lemma 2.

Coefficients $a_{n k}$ have the double ordinary generating function $G (x, y)$ and the semi-exponential generating function $F (x, y)$ as follows: $\begin{aligned} (i) G (x, y) & = {((1 - 2 x) (1 - x (y + 1)))}^{- 1}, \end{aligned}$ (9) $\begin{aligned} (ii) F (x, y) & = (2 e^{2 x} - (1 + y) e^{(y + 1) x}) {(1 - y)}^{- 1} . \end{aligned}$ (10)

Proof.

Let us take the double ordinary generating function (5). Substituting the expression (3) into the generating function, we get $\begin{aligned} G (x, y) & = 1 + \sum_{n = 1}^{\infty} a_{n 0} x^{n} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} a_{n k} x^{n} y^{k} \\ = G_{0} (x) + \sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} a_{n - 1, k - 1} x^{n} y^{k} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} a_{n - 1, k} x^{n} y^{k} \\ = G_{0} (x) + x y G (x, y) + x (G (x, y) - G_{0} (x)) . \end{aligned}$ Thus, $\begin{matrix} G (x, y) = G_{0} (x) (1 - x) {((1 - x (y + 1)))}^{- 1} . \end{matrix}$ (11) 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 $\begin{aligned} F (x, y) & = \sum_{n = 1}^{\infty} \sum_{k = 1}^{n} a_{n - 1, k - 1} \frac{x^{n}}{n!} y^{k} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{n} a_{n - 1, k} \frac{x^{n}}{n!} y^{k} + \sum_{n = 0}^{\infty} a_{n 0} \frac{x^{n}}{n!} \\ = \underset{= y \int_{0}^{x} F d t}{\underset{︸}{\sum_{n = 0}^{\infty} \sum_{k = 0}^{n} a_{n k} \frac{x^{n + 1}}{(n + 1)!} y^{k + 1}}} + \underset{= \int_{0}^{x} F d t}{\underset{︸}{\sum_{n = 0}^{\infty} \sum_{k = 0}^{n} a_{n k} \frac{x^{n + 1}}{(n + 1)!} y^{k}}} \\ - \underset{= \int_{0}^{x} F_{0} d t}{\underset{︸}{\sum_{n = 0}^{\infty} a_{n 0} \frac{x^{n + 1}}{(n + 1)!}}} + \underset{= F_{0}}{\underset{︸}{\sum_{n = 0}^{\infty} a_{n 0} \frac{x^{n}}{n!}}} . \end{aligned}$ Therefore, by statement (ii) of Lemma 1, we obtain that $\begin{matrix} F = (y + 1) \int_{0}^{x} F d t - \underset{= e^{2 x} - e^{x}}{\underset{︸}{\int_{0}^{x} F_{0} d t}} + \underset{= 2 e^{2 x} - e^{x}}{\underset{︸}{F_{0}}} = (y + 1) \int_{0}^{x} F d t + e^{2 x} . \end{matrix}$ Calculating the partial derivative by x, we get the first-order linear partial differential equation, $\begin{matrix} F_{x}^{'} - (y + 1) F = 2 e^{2 x}, \end{matrix}$ with the initial condition $F |_{x = 0} = 1$ . Indeed, by (6), we get $\begin{matrix} F (0, y) = a_{00} + \sum_{n = 1}^{\infty} \sum_{k = 1}^{n} a_{n k} \frac{x^{n}}{n!} y^{k} |_{x = 0} = 1 . \end{matrix}$ Solving the partial differential equation (note that we can solve it as an ordinary one), we obtain statement (ii) of the lemma. □

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

Lemma 3.
The following analytic expressions are equivalent: $\begin{aligned} (i) a_{n k} & = \sum_{j = k}^{n} 2^{n - j} C_{j}^{k}, \\ (ii) a_{n k} & = \sum_{j = k + 1}^{n + 1} C_{n + 1}^{j}, \\ (iii) a_{n k} & = 2^{n + 1} I_{1 / 2} (k + 1, n - k + 1) . \end{aligned}$ (12) Here $I_{x} (a, b)$ stands for the regularized incomplete beta function: $\begin{matrix} I_{x} (a, b) = \frac{\int_{0}^{x} t^{a - 1} {(1 - t)}^{b - 1} d t}{\int_{0}^{1} t^{a - 1} {(1 - t)}^{b - 1} d t} . \end{matrix}$
Proof.
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 $\begin{matrix} G (x, y) = \sum_{n = 0}^{\infty} \sum_{k = 0}^{\infty} ({\frac{\partial^{n + k}}{\partial x^{n} \partial y^{k}} G (x, y)|}_{(0, 0)}) \frac{x^{n} y^{k}}{n! k!} \end{matrix}$ and $\begin{matrix} F (x, y) = \sum_{n = 0}^{\infty} \sum_{k = 0}^{\infty} ({\frac{\partial^{n + k}}{\partial x^{n} \partial y^{k}} F (x, y)|}_{(0, 0)}) \frac{x^{n} y^{k}}{n! k!}, \end{matrix}$ 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 $\begin{aligned} a_{n k} & = {\frac{1}{n! k!} \frac{\partial^{n + k}}{\partial x^{n} \partial y^{k}} G (x, y)|}_{(0, 0)} = \frac{1}{n! k!} \frac{\partial^{n}}{\partial x^{n}} {\frac{1}{1 - 2 x} \frac{k! x^{k}}{{(1 - x (y + 1))}^{k + 1}}|}_{(0, 0)} \\ = \frac{1}{n!} \sum_{j = 0}^{n} C_{n}^{j} {(x^{k})}^{(j)} {{(\frac{1}{(1 - 2 x) {(1 - x (y + 1))}^{k + 1}})}^{(n - j)}|}_{(0, 0)} \\ = \frac{1}{(n - k)!} \sum_{j = 0}^{n - k} C_{n - k}^{j} {(\frac{1}{1 - 2 x})}^{(j)} {{(\frac{1}{{(1 - x (y + 1))}^{k + 1}})}^{(n - k - j)}|}_{(0, 0)} \\ = \sum_{j = 0}^{n - k} \frac{1}{j! (n - k - j)!} j! 2^{j} (k + 1) (k + 2) \dots (n - j) = \sum_{j = 0}^{n - k} 2^{j} C_{n - j}^{k} . \end{aligned}$

Next, using proposition (ii) of Lemma 2, we get $\begin{aligned} a_{n k} & = {\frac{1}{k!} \frac{\partial^{n + k}}{\partial x^{n} \partial y^{k}} F (x, y)|}_{(0, 0)} = \frac{1}{k!} \frac{\partial^{k}}{\partial y^{k}} {\frac{2^{n + 1} e^{2 x} - {(1 + y)}^{n + 1} e^{(y + 1) x}}{1 - y}|}_{(0, 0)} \\ = \frac{1}{k!} \sum_{j = 0}^{k} C_{k}^{j} {\frac{\partial^{j}}{\partial y^{j}} 2^{n + 1} e^{2 x} - {(1 + y)}^{n + 1} e^{(y + 1) x}|}_{(0, 0)} {\frac{\partial^{k - j}}{\partial y^{k - j}} \frac{1}{1 - y}|}_{(0, 0)} \\ = - \sum_{j = 1}^{k} \frac{1}{j!} {\sum_{i = 0}^{j} C_{j}^{i} \frac{\partial^{i}}{\partial y^{i}} {(1 + y)}^{n + 1}|}_{(0, 0)} \frac{\partial^{j - i}}{\partial y^{j - i}} {e^{(y + 1) x}|}_{(0, 0)} + 2^{n + 1} - 1 \\ = - \sum_{j = 1}^{k} \frac{(n + 1)!}{j! (n + 1 - j)!} + 2^{n + 1} - 1 = 2^{n + 1} - \sum_{j = 0}^{k} C_{n + 1}^{j} = \sum_{j = k + 1}^{n + 1} C_{n + 1}^{j} . \end{aligned}$ (13)

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), $\begin{matrix} \sum_{k = m}^{n} C_{n}^{k} p^{k} {(1 - p)}^{n - k} = I_{p} (m, n - m + 1) . \end{matrix}$ Substituting $p = 1 / 2$ , $n = r + 1$ and $m = j + 1$ , we have that $\begin{matrix} I_{1 / 2} (j + 1, r - j + 1) = 2^{- r - 1} \sum_{k = j + 1}^{r + 1} C_{r + 1}^{k} . \end{matrix}$ Thus (cf. (13)), $\begin{matrix} \underset{= a_{n k}}{\underset{︸}{\sum_{j = k + 1}^{n + 1} C_{n + 1}^{j}}} = 2^{n + 1} I_{1 / 2} (k + 1, n - k + 1), \end{matrix}$ concluding the proof. □
Remark 1.
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 $a_{n k}$ .

Now we can proceed with the proof of Theorem 1.
5 Proof of Theorem 1

Proof.
Consider Lerch’s series (4) for the Riemann zeta function $\begin{matrix} ζ (s) = \frac{1}{1 - 2^{1 - s}} \sum_{n = 0}^{\infty} \sum_{k = 0}^{n} \frac{1}{2^{n + 1}} \frac{{(- 1)}^{k} C_{n}^{k}}{{(k + 1)}^{s}} . \end{matrix}$ Changing the order of the summation in the double sum and applying the first statement of Lemma 3 (12), we obtain that $\begin{aligned} ζ (s) & = \frac{1}{1 - 2^{1 - s}} \sum_{k = 0}^{\infty} \sum_{n = k}^{\infty} \frac{1}{2^{n + 1}} \frac{{(- 1)}^{k} C_{n}^{k}}{{(k + 1)}^{s}} \\ = \frac{1}{1 - 2^{1 - s}} lim_{N \to \infty} \sum_{k = 0}^{N} \frac{{(- 1)}^{k}}{{(k + 1)}^{s}} \sum_{n = k}^{N} \frac{C_{n}^{k}}{2^{n + 1}} \\ = \frac{1}{1 - 2^{1 - s}} lim_{N \to \infty} \sum_{k = 0}^{N} \frac{{(- 1)}^{k}}{{(k + 1)}^{s}} 2^{- N - 1} \underset{= a_{N k}}{\underset{︸}{\sum_{n = k}^{N} 2^{N - n} C_{n}^{k}}} \end{aligned}$ yielding us the statement of the theorem. □
6 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 $q_{n k}$ coefficients and specify the rate of convergence to the limiting distribution.

Theorem 2.
(See Hwang , 1998 .) Let $P_{n} (z)$ be a probability generating function of the random variable $Ω_{n}$ , taking only non-negative integral values, with mean $μ_{n}$ and variance $σ_{n}^{2}$ . Suppose that, for each fixed $n ⩾ 1$ , $P_{n} (z)$ is a Hurwitz polynomial. If $σ_{n} \to \infty$ , then $Ω_{n}$ satisfies $\begin{matrix} P (\frac{Ω_{n} - μ_{n}}{σ_{n}} < x) = Φ (x) + O (\frac{1}{σ_{n}}), x \in R . \end{matrix}$ (14)

Let us formulate a central limit theorem.
Lemma 4.
Suppose that $F_{n} (x)$ is the cumulative distribution function of the random variable $A_{n}$ with the probability mass function $\begin{matrix} P (A_{n} = k) = \frac{C_{n + 1}^{k}}{2^{n + 1} - 1}, k = 0, \dots, n, \end{matrix}$ (15) mean $\begin{matrix} μ_{n} = \frac{n}{2} (1 + \frac{1}{n} + O (\frac{1}{2^{n}})) \end{matrix}$ (16) and variance $\begin{matrix} σ_{n}^{2} = \frac{n}{4} (1 + \frac{1}{n} + O (\frac{n}{2^{n}})) . \end{matrix}$ (17) Then it holds $\begin{matrix} F_{n} (σ_{n} x + μ_{n}) = Φ (x) + O (\frac{1}{\sqrt{n}}), x \in R . \end{matrix}$ (18)
Proof.
Let us consider the moment generating function of the random variable $A_{n}$ $\begin{aligned} M_{n} (s) & = E (e^{s A_{n}}) = \sum_{k} P (A_{n} = k) e^{k s} = \sum_{k = 0}^{n} \frac{C_{n + 1}^{k} e^{k s}}{2^{n + 1} - 1} \\ = \frac{1}{2^{n + 1} - 1} (\sum_{k = 0}^{n + 1} C_{n + 1}^{k} e^{k s} - e^{s (n + 1)}) = \frac{{(e^{s} + 1)}^{n + 1} - e^{s (n + 1)}}{2^{n + 1} - 1} . \end{aligned}$ Calculating the derivatives of $M_{n} (s)$ , we obtain that $\begin{aligned} M_{n}^{'} (s) & = \frac{n + 1}{2^{n + 1} - 1} ({(e^{s} + 1)}^{n} e^{s} - e^{s (n + 1)}), \\ M_{n}^{″} (s) & = \frac{n + 1}{2^{n + 1} - 1} (n {(e^{s} + 1)}^{n - 1} e^{2 s} + {(e^{s} + 1)}^{n} e^{s} - (n + 1) e^{s (n + 1)}) . \end{aligned}$ Hence, $\begin{aligned} μ_{n} & = M_{n}^{^{'}} (0) = \frac{(n + 1) (2^{n} - 1)}{2^{n + 1} - 1} = \frac{n + 1}{2} (1 - \sum_{k = 1}^{\infty} {(\frac{1}{2^{n + 1}})}^{k}) \\ = \frac{n + 1}{2} (1 + O (\frac{1}{2^{n}})) = \frac{n}{2} (1 + \frac{1}{n} + O (\frac{1}{2^{n}})) \end{aligned}$ and $\begin{aligned} σ_{n}^{2} & = M_{n}^{″} (0) - {(M_{n}^{'} (0))}^{2} = \frac{n + 1}{{(2^{n + 1} - 1)}^{2}} (2^{2 n} - 2^{n - 1} (n + 2)) \\ = \frac{n + 1}{4} (1 - \sum_{k = 1}^{\infty} \frac{k n + k - 1}{2^{k n + k}}) = \frac{n}{4} (1 + \frac{1}{n} + O (\frac{n}{2^{n}})), \end{aligned}$ granting us with estimates for the mean and the variance (cf. (16) and (17)). Note that $σ_{n} \to \infty$ .

Next, let us examine the probability generating function of the random variable $A_{n}$ $\begin{matrix} P_{n} (z) = E (z^{A_{n}}) = \sum_{k} P (A_{n} = k) z^{k} = \sum_{k = 0}^{n} \frac{C_{n + 1}^{k} z^{k}}{2^{n + 1} - 1} = \frac{{(z + 1)}^{n + 1} - z^{n + 1}}{2^{n + 1} - 1} . \end{matrix}$ (19) 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 \neq 0$ , we have that $\begin{matrix} {(\frac{z + 1}{z})}^{n + 1} = 1 . \end{matrix}$ Calculating an nth root of unity, we obtain that $\begin{array}{r} \frac{z_{k} + 1}{z_{k}} = cos \frac{2 π k}{n + 1} + i sin \frac{2 π k}{n + 1}, k = 0, \dots, n . \end{array}$ Thus, the roots of (19) are equal to $\begin{matrix} z_{k} = - \frac{1}{2} - \frac{i}{2} cot \frac{π k}{n + 1} . \end{matrix}$ The real part of every root is negative. Hence, the probability generating function $P_{n} (z)$ is a Hurwitz polynomial. By Theorem 2, we have that $\begin{matrix} P (\frac{A_{n} - μ_{n}}{σ_{n}}) - Φ (x) = O (\frac{1}{σ_{n}}) = O (\frac{1}{\sqrt{n}}), \end{matrix}$ yielding us the statement of the lemma. □

Now we can formulate the central limit theorem for the coefficients $q_{n k}$ . Theorem 3.
Under conditions of Lemma 4 , the coefficients $q_{n k}$ satisfy $\begin{matrix} q_{n k} = \overline{Φ} (\frac{k - μ_{n}}{σ_{n}}) + O (\frac{1}{\sqrt{n}}) . \end{matrix}$ (20)
Proof.
Note that the cumulative distribution function of the random variable $A_{n}$ equals $\begin{matrix} F_{n} (σ_{n} x + μ_{n}) = \sum_{j ⩽ σ_{n} x + μ_{n}} \frac{C_{n + 1}^{j}}{2^{n + 1} - 1} . \end{matrix}$ Denoting $k = ⌊ σ_{n} x + μ_{n} ⌋$ and taking into account Lemma 4 and the equality (cf. (13)) $\begin{matrix} \frac{a_{n k}}{2^{n + 1}} = 1 - \sum_{j = 0}^{k} \frac{C_{n + 1}^{j}}{2^{n + 1}}, \end{matrix}$ we obtain that $\begin{matrix} \frac{2^{n + 1}}{2^{n + 1} - 1} \sum_{j = 0}^{k} \frac{C_{n + 1}^{j}}{2^{n + 1}} = \frac{2^{n + 1}}{2^{n + 1} - 1} (1 - \frac{a_{n k}}{2^{n + 1}}) = Φ (\frac{k - μ_{n}}{σ_{n}}) + O (\frac{1}{\sqrt{n}}) . \end{matrix}$ Thus, $\begin{matrix} \underset{= q_{n k}}{\underset{︸}{\frac{a_{n k}}{2^{n + 1}}}} = \underset{= \overline{Φ} (\frac{k - μ_{n}}{σ_{n}}) + O (\frac{1}{\sqrt{n}})}{\underset{︸}{1 - \frac{2^{n + 1} - 1}{2^{n + 1}} Φ (\frac{k - μ_{n}}{σ_{n}}) + O (\frac{1}{\sqrt{n}})}}, \end{matrix}$ and the statement of the theorem follows. □

7 Error Term of the Series and the Computational Complexity

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

Theorem 4.
For $σ > 0$ and $s \neq 1 + 2 i π m / log 2$ , $m \in Z$ , $\begin{matrix} ζ (s) = \frac{1}{1 - 2^{1 - s}} \sum_{k = 0}^{n} \frac{{(- 1)}^{k}}{{(k + 1)}^{s}} q_{n k} + ξ_{n + 1} (s) . \end{matrix}$ (21) Here $\begin{matrix} | ξ_{n + 1} (s) | ⩽ \frac{(1 + | t | / σ) e^{π | t | / 2}}{2^{n + 1} | 1 - 2^{1 - s} |} . \end{matrix}$ (22)
Proof.
By (2.3) of Algorithm 1 in Borwein (2000), the error term of the series equals $\begin{matrix} ξ_{n + 1} (s) = \frac{1}{p_{n + 1} (- 1) (1 - 2^{1 - s})} \frac{1}{Γ (s)} \int_{0}^{1} \frac{p_{n + 1} (x) | log x |^{s - 1}}{1 + x} d x . \end{matrix}$ (23) By the definition of the coefficients of the series and the expression of the generating function of the binomial coefficients, the polynomial $p_{n + 1} (x)$ equals $\begin{matrix} p_{n + 1} (x) = \sum_{k = 0}^{n + 1} {(- 1)}^{k} C_{n + 1}^{k} x^{k} = {(1 - x)}^{n + 1} . \end{matrix}$ Thus, $p_{n + 1} (- 1) = 2^{n + 1}$ and, for $0 < x < 1$ , we have $0 < p_{n + 1} (x) < 1$ . It follows $\begin{matrix} | ξ_{n + 1} (s) | ⩽ \frac{1}{2^{n + 1} | 1 - 2^{1 - s} |} \frac{1}{| Γ (s) |} \int_{0}^{1} \frac{{(- log x)}^{σ - 1}}{1 + x} d x . \end{matrix}$ (24) For the integral factor, we have that $\begin{matrix} \int_{0}^{1} \frac{{(- log x)}^{σ - 1}}{1 + x} d x ⩽ \int_{0}^{1} {(- log x)}^{σ - 1} d x = Γ (σ) . \end{matrix}$ (25) The product representation of the gamma function equals $\begin{matrix} {| \frac{Γ (σ)}{Γ (s)} |}^{2} = \prod_{n = 0}^{\infty} (1 + \frac{t^{2}}{{(σ + n)}^{2}}) . \end{matrix}$ (26) The product representation of the hyperbolic sine equals $\begin{matrix} sinh π t = π t \prod_{n = 1}^{\infty} (1 + \frac{t^{2}}{n^{2}}) . \end{matrix}$ (27) Combining (26) and (27), we get $\begin{matrix} {| \frac{Γ (σ)}{Γ (s)} |}^{2} = (1 + \frac{t^{2}}{σ^{2}}) \prod_{n = 1}^{\infty} (1 + \frac{t^{2}}{{(σ + n)}^{2}}) ⩽ (1 + \frac{t^{2}}{σ^{2}}) \frac{sinh π | t |}{π | t |} . \end{matrix}$ (28) Next, combining (24), (25) and (28), we obtain that $\begin{matrix} | ξ_{n + 1} (s) | ⩽ \frac{1}{2^{n + 1} | 1 - 2^{1 - s} |} \underset{⩽ (1 + | t | / σ) e^{π | t | / 2}}{\underset{︸}{{(1 + \frac{t^{2}}{σ^{2}})}^{1 / 2} {(\frac{sinh π | t |}{π | t |})}^{1 / 2}}}, \end{matrix}$ (29) yielding us the statement of the theorem. □
Remark 2.
For $σ > 0$ and $s \neq 1 + 2 i π m / log 2$ , $m \in Z$ , the inequality $\begin{matrix} | ξ_{n + 1} (s) | ⩽ \frac{1}{2^{n + 1} | 1 - 2^{1 - s} |} \frac{Γ (σ)}{| Γ (s) |} \end{matrix}$ (30) 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), $\begin{aligned} n & = ⌊ \frac{log Γ (σ) - log | Γ (s) | + d log 10 - log | 1 - 2^{1 - s} |}{log 2} ⌋ \\ = ⌊ \frac{log Γ (σ) - log | Γ (s) | + d log 10 - log \sqrt{1 - 2^{2 - σ} cos (t log 2) + 2^{2 - 2 σ}}}{log 2} ⌋ \\ ⩽ ⌊ \frac{log Γ (σ) - log | Γ (σ + i t) | + d log 10 - log | 1 - 2^{1 - σ} |}{log 2} ⌋ . \end{aligned}$ (31) Noticing that (cf. 8.328.1, Gradshteyn and Ryzhik, 2014), $\begin{matrix} lim_{t \to \infty} | Γ (σ + i t) | t^{\frac{1}{2} - σ} e^{\frac{π}{2} t} = \sqrt{2 π}, \end{matrix}$ we obtain, for fixed d and σ, that $\begin{matrix} n = \frac{π t}{2 log 2} (1 + O (\frac{log t}{t})) . \end{matrix}$ (32) For numerical purposes (if t is large enough, e.g. $t > 0.2$ ), we can choose n by the following approximate formula $\begin{matrix} n = ⌈ \frac{π t}{2 log 2} + (\frac{1}{2} - σ) \frac{log t}{log 2} + C_{σ} + C_{d} ⌉ . \end{matrix}$ (33) Here $\begin{matrix} C_{σ} = \frac{log Γ (σ) - log | 1 - 2^{1 - σ} |}{log 2}, C_{d} = \frac{log 10}{log 2} d - \frac{log 2 π}{2 log 2} . \end{matrix}$

8 Numerical Aspects

Apart from a certain theoretical interest, the representation (20) of $q_{n k}$ 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_{n k}$ while $k \to 0$ and $q_{n k} \to 1$ , and do not include corresponding terms into the sum, while $k \to n$ and $q_{n k} \to 0$ . The process is controlled by choosing the accuracy level ε. For the right tail, we have $Φ ((k - μ_{n}) / σ_{n}) > 1 - ε$ . Solving the inequality we obtain that $\begin{matrix} n_{1} = ⌈ μ_{n} + z_{1 - ε} σ_{n} ⌉ . \end{matrix}$ (34) For the left tail, we have $n_{0} = ⌊ μ_{n} - z_{1 - ε} σ_{n} ⌋$ . Here $z_{p}$ 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 $n^{2}$ storage).

Table 2

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

Alg. mod.	$H = 10^{2}$ and $(σ, t) \in Σ_{i} \times T_{1}$			$H = 10^{3}$ and $(σ, t) \in Σ_{i} \times T_{2}$
	$Σ_{0} \times T_{1}$	$Σ_{1} \times T_{1}$	$Σ_{2} \times T_{1}$	$Σ_{0} \times T_{2}$	$Σ_{1} \times T_{2}$	$Σ_{2} \times T_{2}$
	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 $n_{1}$	0.18	0.18	0.16	90.20	97.51	87.74
BF $n_{1}$	0.10	0.10	0.09	6.63	6.96	6.57
NA $n_{1}$	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 $n_{1}$ (RE $n_{1}$ ). The coefficients of the series (21) are calculated using the recurrent expression (3) and the number of terms $n_{1}$ (34) in the sum.

•

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

•

Normal approximation with $n_{1}$ (NA $n_{1}$ ). Coefficients of the series (21) are calculated using the normal approximation (20) and the number of terms $n_{1}$ (34) in the sum.

Let us denote the sets for σ: $Σ_{0} = (0.7)$ , $Σ_{1} = (1 / 2, 1)$ , $Σ_{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 = (σ, t)$ of the Riemann zeta function from the sets $Σ_{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_{n k}$ for a mesh are precalculated.

9 3D visualizations of Zeta Functios

Visualizing a complex function of a complex variable $f_{C} : C \to 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_{R} : R^{2} \to 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 $ℜ (\tilde{ζ})$ and the imaginary $ℑ (\tilde{ζ})$ surfaces of a zeta function, along with 2D lines of intersection between the real $ℜ (\tilde{ζ})$ and the imaginary $ℑ (\tilde{ζ})$ 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_{n k}$ . 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) \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 $σ = 1 / 2$ .

Fig. 1

Zeroes, $ℜ ζ (s)$ and $ℑ ζ (s)$ surfaces for $(σ, t) \in (- 1 / 2, 3 / 2) \times (10, 35)$ .

Figure 2 visualizes a segment of the pattern, containing 91–100th non-trivial zeroes of the Riemann zeta function (for $(σ, t) \in (- 1 / 2, 3 / 2) \times (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.

Fig. 2

Zeroes, $ℜ ζ (s)$ and $ℑ ζ (s)$ surfaces for $(σ, t) \in (- 1 / 2, 3 / 2) \times (10, 35)$ .

Fig. 3

$ℜ ζ (s)$ and $ℑ ζ (s)$ surfaces for $(σ, t) \in (- 40, 10) \times (- 20, 100)$ .

Figure 3 shows numerous ridges shaped by the real (yellow) and the imaginary (red) surfaces of the Riemann zeta function (for $(σ, t) \in (- 40, 10) \times (- 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.

Fig. 4

2D curves corresponding the intersections of $ℜ ζ (s)$ and $ℑ ζ (s)$ surfaces with the complex plane, $(σ, t) \in (- 40, 10) \times (- 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) \in (- 40, 10) \times (- 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).

Fig. 5

$ℜ ζ (s)$ and $ℑ ζ (s)$ surfaces for $(σ, t) \in (- 40, 10) \times (- 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.

Fig. 6

3D curves corresponding the intersections of $ℜ ζ (s)$ and $ℑ ζ (s)$ surfaces for $(σ, t) \in (- 40, 10) \times (- 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).

Fig. 7

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), $\begin{matrix} f (s, τ) = (1 - τ) (1 + \frac{\sqrt{5}}{5^{s}}) ζ (s) + τ L (s, χ), \end{matrix}$ (35) where $τ \in [0, 1]$ and $L (s, χ)$ is the Dirichlet L-function (χ mod 5, $χ (2) = - 1$ ). Hence, $\begin{matrix} L (s, χ) = \sum_{n = 1}^{\infty} \frac{χ (n)}{n^{s}} = \sum_{n = 1}^{\infty} (\frac{1}{{(5 n - 4)}^{s}} - \frac{1}{{(5 n - 3)}^{s}} - \frac{1}{{(5 n - 2)}^{s}} + \frac{1}{{(5 n - 1)}^{s}}) . \end{matrix}$ 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) \in (- 1, 2) \times (1.65, 1.95)$ ).

Fig. 8

$ℜ f (s, 3 / 4)$ and $ℑ f (s, 3 / 4)$ surfaces for $(σ, 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 < σ < 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).

Fig. 9

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

References

Arias de Reyna

(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

(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

(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

Sakalauskas

(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

(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

(2017). The Zetafast algorithm for computing zeta functions. arXiv:1703.01414.

10.

Garunkštis

Šimėnas

(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.

11.

Gradshteyn

I.S.

Ryzhik

I.M.

(2014). Table of Integrals, Series, and Products (8th ed.). Academic Press.

12.

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.

13.

Kaczorowski

Kulas

(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.

14.

Lerch

(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.

15.

Šleževičienė

(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.

16.

Stopple

(2017). Lehmer pairs revisited. Experimental Mathematics, 26(1), 45–53. https://doi.org/10.1080/10586458.2015.1107870.

17.

Vepštas

(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.

18.

West

(2005). Better approximations to cumulative normal functions. Wilmott Magazine, 70–76.

	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	…
…	…	…	…	…	…	…	…

	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	…
…	…	…	…	…	…	…	…

	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	…
…	…	…	…	…	…	…	…