Numerical calculation of the Riemann zeta function at odd integer arguments: A direct formula method

In this article, we introduce a recurrence formula which only involves two adjacent values of the Riemann zeta function at integer arguments. Based on the formula, an algorithm to evaluate $\zeta$-values(i.e. the values of Riemann zeta function) at odd-integers from the two nearest $\zeta$-values at even-integers is posed and proved. The behavior of the error bound is $O(10^{-n})$ approximately where $n$ is the argument. Our method is especially powerful for the calculation of Riemann zeta function at large argument, while for smaller ones it can also reach spectacular accuracies such as more than ten decimal places.


Introduction.
Zeta functions of various kinds, such as Hurwitz zeta function, Epstein zeta function and Dirichlet L-function, are all-pervasive objects in modern mathematics, especially in analytical number theory, and among which the prototype zeta function is the famous Riemann zeta function. It is classically defined as the sum of the infinite series [1,2]  The Riemann zeta function, as a function of a real argument, was studied by Euler and Bernoulli brothers in the first half of the 18th century before Riemann extended the definition to a complex variable in 1859 [3]. Historically, people tend to study the closed form of the Riemann zeta function in that those special values seem to dictate the properties of the objects they associated. In 1735, Euler gave the solution of the so-called Basel problem [4] (1. 3) ζ(2) = π 2 6 .
It is well known that for positive even integer arguments the Riemann zeta function can be expressed in a brief form: [5] (1.4) ζ(2n) = (−1) n+1 (2π) 2n 2(2n)! B 2n in terms of the Bernoulli numbers B n . On the contrary, for positive odd integer arguments the Riemann zeta function only has an integral expression: [6] (1.5) The celebrated Goldbach-Euler theorem states that the sum over the set of perfect powers(excluding 1 and omitting repetitions) converges to 1, namely [7,8] For the Riemann zeta function, Goldbach-Eulers theorem (1.6) assumes the elegant form denotes the fractional part of the real number x. In addition, we can show that frac(ζ(2k + 1)) = 1 4 .
Anyway, the formula (1.4) and (1.5), along with (1.8) do reveal somewhat similarity for the values of Riemann zeta function at even and odd arguments. Even up to now, the explicit formula for Riemann zeta function at odd values is difficult if not fundamentally impossible to obtain. Besides, The calculation of Riemann zeta function and related series is a hot topic in computational mathematics. The traditional methods are Euler-Maclaurin formula and Riemann-Siegel formula. Algorithms are still being developed in earnest ever since [9,10,11].
In this paper we mainly obtain a recurrence formula (3.2) relating to the Riemann zeta function and based on which we construct an algorithm for the calculation of the Riemann zeta function at odd integers. In addition, numerical calculation implies that the algorithm can reach considerable accuracies with small odd integer arguments, not to speak of larger ones. Quantificationally, the behavior of the error bound is O(10 −n ) where n is the argument.

Notations and Preliminaries.
We begin by recalling the definition of the Bernoulli polynomials B n (x) and their basic properties in a nutshell to render the paper essentially self-contained. The generating function of the Bernoulli polynomials B n (x) is [12,13] Taking a derivative with respect to x on both sides of (2.1), we find that . Bernoulli polynomials can also be expressed explicitly from Bernoulli numbers For convenience, we introduce two kind of reduced Bernoulli numbers(RBNs), one relates to the even-labeled Bernoulli numbers (denoted by +) , and another relates to the odd-labeled Bernoulli polynomials (denoted by -) In this section we will demonstrate the asymptotic representation of the two kinds of RBNs in a uniform way and establish their integral representation subsequently.
2.1. Asymptotic representations of RBNs. The asymptotic expressions of Bernoulli polynomials at even and odd subscript are respectively [14] (2.6a) Moreover, the Bernoulli polynomials can also be expressed in a stronger form basing on Fourier sine and cosine series expansion [15] (2.7a) Obviously (2.6a)(or (2.6b)) is just the corollary of (2.7a)(or (2.7b)). From the analysis above we find the asymptotic behavior of RBNs are where we have used the fact that 1 0 sin(2πx) cot(πx)dx = 1. The asymptotic representations clearly reveal somewhat similarity of the two kinds of reduced Bernoulli numbers.

2.2.
Integral representations of RBNs. Let us consider two auxiliary integrals [5] (2.9a) with m and n are integers and n ≥ 1. Specially, when n = 1, direct computation shows that .. By virtue of (2.2) and integrating by parts twice, readily yields Combining (2.10a),(2.10b) and (2.11a),(2.11b) we find that hold if m is even. Immediately, the integral representations of the two RBNs B + n and B − n are (2.13a) 3. Algorithm to calculate the Riemann zeta function.
Mathematically, Riemann zeta function is said to be monotonically decreasing since its values are only falling and never rising with increasing values of s with For brevity we denote the reciprocal function as below Now that what we concerned most is the values of the Riemann zeta function at integers for the moment, the asymptotic behavior of the ratio of the reciprocal function (3.1) at odd integers and even integers interests us. Motivated by (1.7) and (1.8), we manage to demonstrate a formula, the so-called recurrence formula(not in a strict sense, though), on condition that the argument is a positive integer.
Motivated by the recurrence, we manage to construct an algorithm to compute the Riemann zeta function.
3.1. Demonstration of the recurrence formula.
Theorem 1. If n is a positive integer such that n ≥ 1, the recurrence relation holds Proof. Making use of (1.4) and (1.5) and the definition of reciprocal function (3.1), we have Ihe limitation of the rightmost term is exactly equal to 1/2 if n is large enough, thus what we want to prove is From the asymptotic formulae (2.8a) and (2.8b) of the two kinds of RBNs we find that, without any difficulty, we have finished the demonstration of the recurrence formula of the Riemann zeta function.

3.2.
Basic ideas for the algorithm. The explicit formula for odd-argument ζvalues remains to be an open problem though some results shed light on it. Since there exits a closed form expression for even-arguments Riemann zeta function, Naturally, one way to construct the algorithm to evaluate the odd-arguments Riemann zeta function is by using of the even-argument ζ-values of the nearest. We know that when the argument is large enough, (3.2) can be rewritten as where ρ l (2n+ 1) and ρ r (2n+ 1) represent two different representation of the asymptotic behavior of ρ(2n + 1). Judging by appearance, One can use any of the formula above to calculate the Riemann zeta function at odd integers. When considering that those two formulae give the upper and lower bound of the zeta-values at odd integers, we come up with the idea that we can combine them together by a special method. Let us give a proposition relating to Dirichlet eta function firstly before we move forward to give another theorem.
Lemma 1. If n is a positive integer such that n ≥ 1, the two inequalities hold Those two inequalities are quite new to the authors because we haven't seen them in any literature or monograph. However, we are not intended to give the details here since the demonstration is rather elementary. The theorem below holds once we take advantage of the lemma.
The basic steps of the calculation of ζ(2n + 1) are presented as follow. Firstly, we calculate ρ(2n) and ρ(2n + 2) from (1.4), (1.2) and (3.1)) in sequence. Secondly, we obtain the value of ρ(2n + 1) from (3.10). Lastly, we calculate ζ(2n + 1) from (3.1)) and (1.2), reversely. In order to start our method, we need to know some ζ-values at even integers. For example, ζ(2) and ζ(4) should be available to get ζ(3). We can obtain ζ(2n) through (1.4) systematically for small argument. However, it is almost impossible to obtain Bernoulli numbers by the ordinary recursive methods thus we can't know the values of ζ(2n) if the argument is large enough. Many methods for computing Bernoulli numbers have been invented. David Harvey introduced an efficient multimodular algorithm [16] and a subquadratic algorithm [17] which ensure us to obtain the Bernoulli numbers B n at n = 10 8 . What's more, he implenmented the algorithms in a C++ package named bernmm(Bernoulli multi-modular). In addition, one can also use the inline function Zeta[s] in Mathematica since it is also based on an efficient algorithm. For convenience, our computation platform is main on Mathematica and we regard the values obtained by it as accuracy values.

Calculation of Riemann zeta function at odd integers.
The calculation of Riemann zeta function plays an essential role in the study of number theory and associated subjects such as statistical physics and string theory. Various approaches to accomplish this task have been proposed [18,19], especially for the evaluation of zeta function at integer arguments or in the critical strip (for the computation of Riemann's zeros). Most of the methods available consist of using integral forms of some particular functions or recursive series forms. Quite recently, Babolian et al transform ζ(s) to appropriate integral forms and introduce a method to compute the Riemann zeta function based on Gauss-Hermite and Gauss-Laguerre quadratures [9]. Besides, many rapidly converging series for ζ(2n+ 1) have been introduced by Srivastava in a review article [20]. In this section we firstly give some numerical examples according to our method to illustrate its accuracy, then we compare our method to two selected ones to show that our method is especially powerful to calculate the ζ-values at large odd integer arguments.

4.1.
Numerical test and error bound of the algorithm. We regard the zetavalues obtained by Mathematica as a benchmark. The result of the Riemann zeta function at odd integers with n = 1, 2, ..., 10 obtained by our method(approximate value) is presented in table (1). The accuracy value and the absolute error are also presented at the same time.  (1) we find that the idea that making the geometric mean instead of any of the "upper" or "lower" bound(see Theorem 2) be the best estimate of the Riemann zeta function dramatically reduces errors and satisfactory accuracy such as twelve decimal places in the tenth odd-argument of the Riemann zeta function can be achieved.
Indeed, what we concerned most is the error bound when the argument n is large enough. In table (2) we present the errors ǫ(n) versus n.  From table (2) we find that the error is of the order O(10 −n ) approximately. By using of least square method, we know that log(ǫ(n)) = −0.9542n − 1.5563. This formula suggests that when the argument of Riemann zeta function is large enough, our algorithm should be very powerful to obtain the zeta-values at odd integers.

4.2.
Compare with the other methods. In this subsection, we aim to compare our algorithm with the already existing ones, namely the Gauss-Hermite quadrature(Integration method, see [9], Corollary 3.1) and rapid converging series(Series method, see [20], eq.(3.30)). The latter two methods rely on the number of zeros of the associated polynomial(in this occasion it is Hermite polynomial) and the terms of partial sum of the infinite series respectively. Since those two method are both efficient, we only make the two integers be the same value 25. The behaviors of the error bound of integration method and series method, as can be seen from table 3, are totally different. When the argument increases, the errors of the former decrease exponential from a high level, while the latter maintain at a nearly constant low level despite of the variation of n. Our method exhibits the worst results for small argument, but the errors decrease dramatically with argument increasing. It outstrips integration method and series method before n = 12 and n = 21 respectively. To reach the accuracy obtained by our method, the number of nodes and items in the above two methods should be augmented largely. For example in the series method, the first cn(c is a constant) terms in the infinity series should be included. Obviously, it is impossible to carry on within the limited CPU time when n is very big.

Summary.
In summary we firstly introduce two kinds of reduced Bernoulli numbers and prove their asymptotic behaviors in a uniform method. Besides, we find a recurrence formula (3.2) of the Riemann zeta function original and construct an algorithm to evaluate the Riemann zeta function at odd integers based on the recurrence formula. The behavior of the error bound ǫ(n) is governed by log(ǫ(n)) = −0.9542n − 1.5563 or ǫ(n) = O(10 −n ) approximately, which, of course, suggests that our method is especially suit for the calculation of ζ-values at large odd integer arguments. However, once we can extend the validity of (3.2) from positive integers to positive real numbers, i.e., where the argument s is an arbitrary positive real number, a new approach to the calculation of the Riemann zeta function will be discovered. The problem above is precisely the direction where some subsequent work shall be done. Remarkably, the recurrence formula (3.2) is likely to act as a touchstone to explore the closed form of the Riemann zeta function at positive integers since it eclipses the mystery of the ζ-values at odd integers into shade.