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 ζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta $$\end{document}-values (i.e., the values of Riemann zeta function) at odd integers from the two nearest ζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta $$\end{document}-values at even integers is posed and proved. The behavior of the error bound is O(10-n)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(10^{-n})$$\end{document} approximately where n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document} 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][3] fðsÞ ¼ with the complex variable s ¼ r þ it. Specially, the series converges if r ¼ Re s [ 1. We can extend fðsÞ from s with Re s [ 1 to s with Re ðsÞ [ 0; s 6 ¼ 1 by the following formula where gðsÞ is the Dirichlet eta function or alternating eta function.
Historically, people prefer to study the closed form of the Riemann zeta function at positive integer arguments in that those special values seem to dictate the properties of the objects they associated. In condensed matter physics for instance, the famous Sommerfeld expansion, which is useful for the calculation of particle number and internal energy of electrons, involves Riemann zeta function at even integers [4], while the spin-spin correlation functions of isotropic spin-1=2 Heisenberg model are expressed by ln 2 and Riemann zeta functions with odd-integer arguments [5]. It was, without doubt, a profound discovery of Euler in 1736 to work out the prolonged Basel problem [6] fð2Þ ¼ superbly. It is well known that for positive even-integer arguments, the Riemann zeta function can be expressed explicitly as [7] fð2nÞ ¼ ðÀ1Þ nþ1 ð2pÞ 2n 2ð2nÞ! B 2n ð4Þ in terms of the Bernoulli numbers B n . On the contrary, however, the explicit formula for Riemann zeta function at odd values is difficult if not fundamentally impossible to obtain. Euler himself once conjectured that fð2n þ 1Þ ¼ cðnÞp 2nþ1 and c involve the irrational constant gð1Þ ¼ ln 2 [8]. This suggests that Riemann zeta function at odd integers produces a recurrence relation that is self-recursive. Even up to now, for positive odd-integer arguments, the Riemann zeta function can only be expressed by series and integral [see (35) and (36) for detail]. One possible integral expression is [9] fð2n þ 1Þ ¼ ðÀ1Þ nþ1 ð2pÞ 2nþ1 2ð2n þ 1Þ!
where B 2nþ1 ðxÞ are Bernoulli polynomials. A relevant aspect is that, for Riemann zeta function, the celebrated Goldbach-Euler theorem [10] assumes the elegant form where fracðxÞ ¼ x À ½x denotes the fractional part of the real number x. It turns out that Indeed, the formulas (4) and (5), along with (7) do reveal somewhat similarity for the values of Riemann zeta function at even and odd arguments. Meanwhile, 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, and algorithms are still being developed in earnest ever since [11][12][13][14]. Typically, a particular numerical method is limited to a special domain. Therefore, when concentrating on Riemann zeta function at odd integers, a special method should be constructed in view of the connection of Riemann zeta function values between odd and even integers. In this paper, we mainly obtain a recurrence formula (22) 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 [1][2][3] te tx e t À 1 Taking a derivative with respect to x on both sides of (8), we find that Bernoulli polynomials can also be expressed explicitly from Bernoulli numbers For convenience, we introduce two kinds 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 framework and establish their integral representation subsequently.

Asymptotic representations of RBNs
The asymptotic expressions of Bernoulli polynomials at even and odd subscript are, respectively [17] ðÀ1Þ nþ1 ð2pÞ 2n 2ð2nÞ! B 2n ðxÞ $ cosð2pxÞ; ð13aÞ Moreover, the Bernoulli polynomials can also be expressed in a stronger form based on Fourier sine and cosine series expansion [18] B 2n ðxÞ ¼ ðÀ1Þ nþ1 2ð2nÞ! ð2pÞ 2n Obviously, (13a, b) is just the corollary of (14a, b). Joining together we, therefore, obtain the asymptotic behavior of RBNs where we use the fact that R 1 0 sinð2pxÞ cotðpxÞdx ¼ 1.

Integral representations of RBNs
Let us consider two auxiliary integrals [7] I c ðn; mÞ ¼ with m and n are integers and n ! 1. Specially, when n ¼ 1, direct computation shows that  By virtue of (9) and integrating by parts twice, readily yields 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 (21) at odd integers and even integers interests us. Motivated by (6) and (7), 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.

Demonstration of the recurrence formula
Theorem 1 If n is a positive integer such that n ! 1, the recurrence relation holds Proof Using (4) and (5) and the definition of reciprocal function (21), we have Since the limitation of the rightmost term is exactly equal to 1=2 if n is large enough, what we want to prove is or equivalently From the asymptotic formulae (15a) and (b) 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. h As a matter of fact, we can extend the validity of (22) from positive integers to positive real numbers straightforward. We, therefore, can obtain the asymptotic behavior of Riemann zeta function as using (22). The application of (26) can be diverse, here we just pick a example relating to prime number theorem. The positive integer x is s-free if and only if in the prime factorization of x, no prime number occurs more than s À 1. Indeed, if Qðx; sÞ denotes the number of s-free integers (e.g., 2-free integers being square-free integers) between 1 and x, one can show that [19] Qðx; sÞ ¼ therefore, we find that the asymptotic density of s-free integers Qðx; sÞ=x $ 1 fðsÞ is nothing but (26). Another intriguing issue is to what degree can (26) reveals its ability to obtain the f-values. Figure 1 is thus plotted as follow.
The fact that all the stars ''*'' lie on the solid curve indicates that (26) may be a suitable candidate for the calculation of Riemann zeta function. The emergence of the abnormal slope between s ¼ 3 and s ¼ 4 in the inserted figure, however, implies that any f-value obtained from its nearest neighbors should be much more accuracy. We, therefore, come up with a satisfactory proposal which is postponed until next subsection.
Basic ideas for the algorithm Abundant methods to evaluate the fð2nÞ have appeared in the mathematical literatures from now and then ever since Euler's seminal work. In contract, the explicit formula for odd-argument f-values remains to be an open problem though some results shed light on it [15,16]. By analogy to fð2nÞ, several authors have established the series and integral representations of fð2n þ 1Þ, which, to some degree, provides some perspectives on the difficulty of evaluating fð2n þ 1Þ as opposed to fð2nÞ. From the viewpoint of numerical method, one natural way to construct the corresponding algorithm to evaluate the oddargument Riemann zeta function is by virtue of the evenargument f-values near to them. In the current paper, only the two nearest f-values are taken into consideration currently for simplicity. When n is large enough, (22) can be rewritten as where q l ð2n þ 1Þ and q r ð2n þ 1Þ represent two different representations of the asymptotic behavior of qð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 (see Theorem 2), we come up with the idea that we can combine them together by a special method. It happens to us that there may exist a somewhat mysterious map from fð2nÞ and fð2n þ 2Þ to fð2n þ 1Þ, which will ensure us to obtain the approximation values of fð2n þ 1Þ with higher precision. Let us give a proposition relating to Dirichlet eta function first before we move forward to give another theorem.
Lemma 1 If n is a positive integer such that n ! 1, the two inequalities hold gð2nÞ Those two inequalities are quite new to the authors because we have not seen them in any literature or monograph before. However, we are not intended to give the details here since the demonstration is rather Theorem 2 If n is a positive integer such that n ! 1, the inequality holds where f l ð2n þ 1Þ and f r ð2n þ 1Þ correspond to q l ð2n þ 1Þ and q r ð2n þ 1Þ, respectively.
Since Riemann zeta function is a monotonic decreasing function, the exact value fð2n þ 1Þ is just between f l ð2n þ 1Þ and f r ð2n þ 1Þ for any given positive integer n. For the benefit of accuracy, we regard the geometric mean values of (28a) and (b) as the approximate values of the reciprocal function qð2n þ 1Þ, namely which is the most valuable ingredient of our algorithm. The basic steps for the calculation of fð2n þ 1Þ are presented as follow. First, qð2nÞ and qð2n þ 2Þ should be calculated from (4), (2) and (21) in sequence. Second, the value of qð2n þ 1Þ is ready to be obtained in light of (32). Lastly, the ultimate aim, i.e., fð2n þ 1Þ is just at hand from (21) and (2), reversely. Our algorithm does not bother circulation of any kind, it just looks like a formula, therefore, we refer it as the direct formula method.
To start our method, we need to know some f-values at even integers. For example, fð2Þ and fð4Þ should be available to get fð3Þ. We can obtain fð2nÞ through (4) systematically for small argument. However, it is almost impossible to obtain Bernoulli numbers by the ordinary recursive methods thus we hardly know the values of fð2nÞ, if the argument is large enough. Many methods for computing Bernoulli numbers have been invented. David Harvey introduced an efficient multimodular algorithm [20] which ensures us to obtain the Bernoulli numbers B n at n ¼ 10 8 . However, one can also use the intrinsic function Zeta[s] in Mathematica since it is also based on an efficient algorithm. Therefore, for convenience, our computation platform is mainly on Mathematica and we regard those values as benchmarks.

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 condensed matter physics. Various approaches to accomplish this task have been proposed [11][12][13]21], 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 fðsÞ to some appropriate integral forms and introduce a method to compute the Riemann zeta function based on Gauss-Hermite and Gauss-Laguerre quadratures [11]. Numerical result shows that 20 points are capable of producing an accuracy of seven-decimal place for small arguments. Besides, many rapidly converging series for fð2n þ 1Þ have been introduced by Srivastava in a review article [13] and by other authors [12,14]. In this section, we first 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 f-values at large oddinteger arguments.

Numerical test and error bound of the algorithm
We regard the f-values obtained by Mathematica as benchmarks. 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 values and the absolute errors are also presented at the same time.
Table 1 tells us 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. It is interesting to find that only the Apéry's constant fð3Þ sightly larger than the approximate value obtained by our method. It is also funny to see the errors present an upsidedown stair configuration, which implies that the error declines about ten times as long as the argument n increase 1.
In Table 2, we present the absolute errors ðnÞ versus n, for the purpose of exploring the error bound when the argument n is large enough. It is clear that, from Table 2, the error is of the order Oð10 Àn Þ approximately. Using least square method, we notice that lgððnÞÞ ¼ À0:9542n À 1:6884: This formula suggests that when the argument of Riemann zeta function is large enough, our algorithm should be powerful enough to obtain the f-values at odd integers.

Comparison with the existed methods
In this subsection, we aim to compare our algorithm with the already existed ones, namely the Gauss-Hermite quadrature (Integral method, see [11], Corollary 3.1) and rapid converging series (Series method, see [13, eq.(3.30)]). The Gauss-Hermite quadrature formula has the form [11] Z where x k is one of the zeros of H N ðxÞ, the Hermite polynomial of degree N, and w k ¼ À 2 Nþ1 N! ffiffi 2 N ð2NÞ! f 2N ðgÞ,g 2 ðÀ1; 1Þ is, obviously, the error bound of the above integral. Riemann zeta function is such an amazing function that it can be transformed into [11] fðsÞ ¼ whose numerator and denominator are of the form presented in (34). Among all the series representations of Riemann zeta function, the series below converges most rapidly as pointed out by Srivastava [13].
When n ¼ 1 for instance, the error bound R ðsÞ N of the N-th partial sum of the infinite series in (36) satisfies jR ðsÞ N j ¼ where we have used the fact that fðsÞ\ 1 1À2 1Às since gðsÞ\1. If N ¼ 25, the error bound is jR ðsÞ 25 j\1:0 Â 10 À20 , which is superior to other rapid series jR 25 j\0:9 Â 10 À18 as noted in [8,12]. Specially, when N is larger than some typical numbers, the asymptotic behavior of (37) reads lgðjR ðsÞ N jÞ $ À 2 lg 2ðN þ 1Þ À 3 lg N: The accuracy of the latter two methods relies on the number of zeros (denoted as N 1 ) of the associated polynomial (in this occasion, it is Hermite polynomial) and the terms (denoted as N 2 ) of partial sum of the infinite series, respectively. We set two integers to be the same value, i.e., N 1 ¼ N 2 ¼ 25 since the corresponding methods are both efficient as have been declared by many Mathematicians. The behaviors of the error bound of integral 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 arguments, but the errors decrease dramatically with argument increasing. It outstrips integration method and series method before n ¼ 12 and n ¼ 21, respectively. Our method  superior to them absolutely afterwards. To reach the accuracy obtained by our method, the number of nodes and terms in the above two methods should be augmented largely. In the series method for instance, the terms of the order n in the infinity series should be included according to (33) and (38). Obviously, it is almost impossible to carry on within the limited CPU time when n is an astronomical number.

Conclusion
In summary, we first introduce two kinds of reduced Bernoulli numbers (RBNs) and prove their asymptotic behaviors in an uniform framework, and their series and integral representations are available at the same time. What is more, we discover and prove a recurrence formula (22) of the Riemann zeta function original and construct an algorithm to evaluate the Riemann zeta function at odd integers based on it. The idea of our method is quiet simple, but it turns out to be a competent algorithm. The behavior of the error bound ðnÞ is governed by lgððnÞÞ ¼ À0:9542n À 1:6884 or ðnÞ ¼ Oð10 Àn Þ approximately, which, of course, suggests that our method is especially suit for the calculation of f-values at large odd-integer arguments. Therefore, our results can also work as benchmarks to test the accuracy of other related algorithms. However, more works should be carried on to improve the accuracy at small arguments in future. Remarkably, the recurrence formula (22) is likely to act as a touchstone to explore the closed form of the Riemann zeta function at positive integers since it witnesses the connection between f-values at odd integers and even integers.