\documentstyle[11pt]{article}
\title{\bf Zeta Function Expansions\\ of Classical Constants}
\author{Philippe Flajolet and Ilan Vardi}
\date{February 18, 1996}
\newcommand{\ds}{\displaystyle}
\begin{document}
\maketitle
Many mathematical constants are expressed as
slowly convergent sums of the form
\begin{equation}\label{z1}
C=\sum_{n\in A} f(\frac{1}{n}),
\end{equation}
for some well-behaved analytic function $f$ and
some ``reasonable'' subset $A$ of the integers.
The convergence of such sums can be accelerated easily
once values at the integers of the zeta function
\begin{equation}\label{z2}
\zeta_A(s):=\sum_{n\in A} \frac{1}{n^s}
\end{equation}
are known: we have, formally at least,
\begin{equation}\label{z3}
C=\sum f_n\zeta_A(n)\quad\hbox{where}\quad
f(z)=\sum_m f_m z^m
\end{equation}
is the Taylor expansion of $f$ at~0.
This scheme is especially effective
in the context of high-precision evaluation of
mathematical constants as it normally exhibits geometric convergence.
It is also very easy to implement in current symbolic manipulation
systems that have built in many mechanisms for the fast computation of
zeta values.
In this note, we show the application of the rearrangement
summarized by (\ref{z1}-\ref{z3}) in the case where $A$
is either the whole set of integers, the prime numbers,
or a congruence subset of the primes. This yields
fast numerical schemes for the evaluation of many constants
like: the number $\pi$, Euler's constant~$\gamma$,
Khinchin's constant $K$, the Hafner-Sarnak-McCurley constant~$\sigma$,
Hardy-Littlewood's twin prime constant~$H$, or the Landau-Ramanujan
constant~$\lambda$. The reader is directed to Finch's beautiful
pages~\cite{Finch95b} for background information on these and other
classical constants.
\medskip
{\em \S1. The number $\pi$ and Gregory's formula.}
Many classical constants can be defined by
a series of the form
\begin{equation}\label{a1}
S=\sum_{n=1}^\infty f(\frac{1}{n}),
\end{equation}
where $f(z)$ is analytic at $0$ and $f(z)=O(z^2)$ there.
If $f(z)$ is analytic the closed unit disk,
then rearranging the series gives
\begin{equation}\label{a2}
S=\sum_{m=2}^\infty f_m \zeta(m),
\end{equation}
where
\[f(z)=\sum_{m=2}^\infty f_m z^m\]
is the Taylor expansion of $f$ and $\zeta(s)=\sum_{n\ge1}n^{-s}$
is the Riemann zeta function.
While the original series~(\ref{a1}) is slowly convergent,
the transformed series~(\ref{a2}) has geometric convergence.
Similarly, if $f(z)$ is analytic in $|z|< 1/m_0$
for some positive integer $m_0$, then
\begin{equation}\label{a3}
S=\sum_{n=1}^{m_0} f(\frac{1}{n})+\sum_{m=2}^\infty f_m
\left[\zeta(m)-\frac{1}{1^m}-\cdots-\frac{1}{m_0}\right],
\end{equation}
a formula that again converges geometrically.
Vardi's book~\cite[p.~156]{Vardi91} gives as an example
Gregory's series for $\pi$,
\[\frac{\pi}{4}
=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots
=\sum_{n=1}^\infty \left(\frac{1}{4n-3}-\frac{1}{4n-1}\right),
\]
whose rearrangement is found to be
\[\pi=\sum_{m=1}^\infty \frac{3^m-1}{4^m}\zeta(m+1),\]
The same technique will apply to Catalan's constant,
\[G=1-\frac{1}{3^2}+\frac{1}{5^2}-\frac{1}{7^2}+\cdots
=\sum_{n=1}^\infty \left(\frac{1}{(4n-3)^2}-\frac{1}{(4n-1)^2}\right).\]
Given the formul{\ae}~(\ref{a2},\ref{a3}),
a sum like~(\ref{a1}) is ``easily''
computable,
once values of the zeta function at the integers
have been tabulated, assuming that
the Taylor expansion of $f$ is itself easily
computable. The table of zeta values can be built from the
Euler-Maclaurin formula~\cite{Edwards74} while
zetas of even argument are
computable directly from Bernoulli numbers.
In this way, Stieltjes~\cite[vol.~II, p.~100]{Stieltjes93} determined
in 1887 the values $\zeta(2),\ldots,\zeta(70)$ to 30 digits of accuracy.
\bigskip
{\em \S2. Euler's constant.}
By definition, Euler's constant is
$$\gamma:=\lim_{n\to\infty} \left(H_n-\log n\right).$$
The limit definition transforms into a sum,
$$\gamma=1+\sum_{n=1}^\infty \left(
\frac{1}{n+1}+\log\frac{n}{n+1}\right),
$$
whose general term converges like $O(n^{-2})$.
The series rearrangement now applies upon taking $m_0=1$,
$$\gamma=\frac{3}{2}-\log 2-\sum_{m=2}^\infty (-1)^m\frac{m-1}{m}
\left[\zeta(m)-1\right].$$
This is one of the many ways to find a geometrically
convergent scheme for the computation of Euler's constant.
An equivalent formula
was already known to Euler
and Stieltjes used precisely the formula corresponding to $m_0=2$
to check his computation of zeta values.
\bigskip
{\em \S3. Khinchin's constant.}
By taking logarithms, infinite products can also be computed.
For instance, Vardi's book~\cite[p.~163]{Vardi91}
discusses the computation of
Khinchin's constant
$$K=\prod_{n=1}^\infty \left(1+\frac{1}{n(n+2)}\right)^{\log n/\log
2}
$$
along these lines. There, because of the logarithms in the exponent,
the rearranged series of $\log K$
involves the values $\zeta'(m)$.
Other rearrangements, like
\[\log K=
\frac{1}{\log2}\sum_{m=1}^\infty \frac{h_{m-1}}{m}(\zeta(2m)-1),\quad
h_m=\sum_{j=1}^m\frac{(-1)^{j-1}}{j}
\]
(due to Shanks and Wrench, 1959)
are known that do not require a special
computation of the derivatives. However, Vardi's method has the advantage
of complete generality as it
applies to sums and products of the form
\[\sum_{n\ge2}\log n f(\frac{1}{n}),\quad
\prod_{n\ge2} \left(1+f(\frac{1}{n})\right)^{\log n},
\]
of which Khinchin's constant is a particular instance.
\bigskip
{\em \S4. Sums over primes and the Hafner-Sarnak-McCurley constant.} A sum of the form
\begin{equation}\label{b1}
T=\sum_{p\ge2} f(\frac{1}{p}),
\end{equation}
where by usual conventions an index $p$ ranges over the prime numbers,
can be rearranged into
\begin{equation}\label{b2}
T=\sum_{m=2}^\infty f_m \Pi(m),
\end{equation}
provided $f$ is analytic in $|z|\le1$. There,
\begin{equation}\label{b3}
\Pi(s):=\sum_{p\ge2} \frac{1}{p^s}.
\end{equation}
If $f(z)$ is analytic in $|z|< 1/m_0$, one has
\begin{equation}\label{b4}
T=\sum_{n=1}^{m_0} f(\frac{1}{n})+\sum_{m=2}^\infty f_m
\left[\Pi(m)-\frac{1}{1^m}-\cdots-\frac{1}{m_0}\right],
\end{equation}
again with geometric convergence.
This technique makes it possible to evaluate $T$ efficiently.
In effect, $\Pi(s)$ is related to the zeta function,
as results from the Eulerian product,
\[\zeta(s)=\prod_p \left(1-\frac{1}{p^{-s}}\right)^{-1},\]
by taking logarithms,
\[\begin{array}{lll}
\log\zeta(s)&=&\ds{\sum_{k\ge1} \frac{1}{k}\sum_{p\ge2}
\frac{1}{p^ks}}\\
&=&\ds{\sum_{k\ge1}\frac{1}{k}\Pi(ks)}.
\end{array}
\]
The last formula is a clear case of application of Moebius inversion
and one finds
\begin{equation}\label{b5}
\Pi(s)=\sum_{k\ge1}\frac{\mu(k)}{k}\log\zeta(ks).
\end{equation}
In summary, $T$ is efficiently computable from the formula~(\ref{b2})
(or its variant~(\ref{b4})) with $\Pi(s)$ given by~(\ref{b5})
which permits to build a table of values of $\Pi(m)$.
Alternatively, one can compute $T$ directly
from the zeta values by combining the Moebius
formula for $\Pi(m)$ and the expansion of $f$:
\begin{equation}\label{b6}
T=\sum_{\nu\ge2} g_\nu \log\zeta(\nu),\qquad
g_\nu =\sum_{k=1}^\nu \frac{\mu(k)}{k} f_{\nu/k}.
\end{equation}
Hafner, Sarnak, and McCurley~\cite{HaSaMc93}
have shown that the probability that two
$m\times m$ and $n\times n$ matrices, $m,n$ large, have relatively
prime determinants is
\[
\sigma=\prod_{p}\left(1 - \left[1 - \prod_{n=1}^\infty
(1-1/p^n)\right]^2\right)
\]
(it is well-known that this is $6/\pi^2$ when
$m=n=1$).
An application of the method to $\log\sigma$~\cite[p.~174]{Vardi91}
gives an equivalent form
\[\sigma=\prod_{m\ge2}\zeta(m)^{-a_m}\]
and, upon taking 100 factors, we get
\[\sigma\approx 0.35323637185499598454\ldots,\]
The speed of convergence is determined by singularity
analysis~\cite[p.~258-261]{Vardi91} and is roughly $0.57^n$.
\bigskip
{\em \S5. The Twin-primes constant.}
Hardy and Littlewood~\cite[p.~371]{HaWr79}
developed in 1923 a heuristic model
for the distribution of prime pairs according to which the number
of prime pairs $p,p+2$ with $p\le x$ must be asymptotic to
\[ 2H {x \over (\log x)^2},\quad
H =\prod_{p\ge3} \left(1-{1\over (p-1)^2}\right).
\]
The same constant also occurs in heuristic models of the Goldbach
problem~\cite{Finch95b}.
The twin prime constant is a direct case of application of the method
of~(\ref{b6}).
It is amusing to note
that the coefficients $g_\nu$ that appear
there are related to finite fields.
Let $I_n$ be the number of monic irreducible
polynomials over the coefficient field $GF(2)$ with
degree $n$. We have~\cite[sec.~3.3]{Berlekamp68}
$$I_n=\frac{1}{n}\sum_{d\;|\;n} \mu(d)2^{n/d}.$$
On the other hand,
\[\log(1-\frac{1}{(p-1)^2})=
-\sum_{m=1}^\infty (2^m-2)\frac{1}{mp^m}.\]
Thus,
\[H=\prod_{n=2}^\infty \left(\zeta(n)(1-2^{-n})\right)^{-I_n},
\]
which converges like $(\frac{2}{3})^n$.
Using a cutpoint $m_0$ like in~(\ref{a3}) gives a whole
collection of formul{\ae} of which
\[H=
\frac{3}{4}\,\frac{15}{16}\,\frac{35}{36}
\prod_{n=2}^\infty \left(\zeta(n)(1-2^{-n})(1-3^{-n})
(1-5^{-n})(1-7^{-n})\right)^{-I_n}
\]
has rate of convergence $\approx(\frac{11}{2})^{-n}$, thus giving
approximately $3/4$ of a digit per iteration.
A closely related example is Mertens's constant~\cite[p.~351]{HaWr79},
\[B_1=\gamma+\sum_p \left(\log(1-p^{-1})-\frac{1}{p}\right),\]
that intervenes in the formula
\[\sum_{p\le x}=\log\log x+B_1+o(1).\]
By the same device, we have
\[
e^{B_1}=\prod_{m=2}^\infty \zeta(m)^{1-\mu(m)/m}.
\]
\bigskip
{\em \S6. The Landau Ramanujan constant.}
This constant is defined by
$$\lambda=\left({1\over2}\prod_{r}{1\over 1-r^{-2}}\right)^{1/2},$$
in which $\prod_r$ indicates a product where $r$ ranges over the prime numbers
that are congruent to $3$ modulo $4$.
It is known from work by Landau in 1908 rediscovered later by
Ramanujan and mentioned in his first letter to Hardy~\cite[Ch.~23]{Berndt94},
\cite{Hardy78} that
the number of integers less than $x$ that are sums of two squares is
asymptotic to
$$\lambda{x\over\sqrt{\log x}}.$$
Ramanujan gave the (correct) approximate value
$\lambda\doteq0.764$.
The general problem behind this
example is to estimate sums
over primes that satisfy congruence restrictions.
Here, we let $p$ range over all primes, $q$ range over primes
$\equiv 1 \bmod4$ and $r$ range over primes $\equiv 3 \bmod 4$.
Euler's product formula is
\begin{equation}\label{e1}
\zeta(s)=(1-2^{-s})^{-1}\prod_q (1-q^{-s})^{-1} \prod_r
(1-r^{-s})^{-1}.
\end{equation}
The function $L(s)$ defined by
$$L(s)=\sum_{n=0}^\infty {(-1)^n\over (2n+1)^s}$$
also admits an Eulerian decomposition
\begin{equation}\label{e2}
L(s)=\prod_q (1-q^{-s})^{-1} \prod_r
(1+r^{-s})^{-1}.
\end{equation}
Thus, comparison of~(\ref{e1}) and~(\ref{e2}) yields
\begin{equation}\label{e3}
(1-2^{-s}){\zeta(s)\over L(s)}=\prod_r{1+r^{-s}\over1-r^{-s}}.
\end{equation}
Taking logarithms in~(\ref{e3}), we get
\begin{equation}\label{e4}
\log \prod_r{1+r^{-s}\over1-r^{-s}}
=\ell(s)\qquad\hbox{where}\qquad \ell(s)=\log
\left((1-2^{-s}){\zeta(s)\over L(s)}\right).
\end{equation}
Introduce the ``base'' function
$$R(s)=\sum_{r} {1\over r^s}.$$
A Taylor expansion of the left hand side of~(\ref{e4}) leads to the
relation
\begin{equation}\label{e5}
R(s)+{1\over3}R(3s)+{1\over5}R(5s)+\cdots = {1\over2}\ell(s),
\end{equation}
which is easily inverted by a variant of Moebius inversion
\begin{equation}\label{e5a}
(s)={1\over2}\sum_{n=1}^\infty {\widehat\mu(n) \over n} \ell(ns),
\end{equation}
where
$$\widehat\mu(n)
=\left\{ \begin{array}{ll}
\mu(n) & \hbox{if $n\equiv 1 \pmod 2$}\\
0 & \hbox{if $n\equiv 0 \pmod 2$}.
\end{array}\right.
$$
Let $f(z)$ be any function analytic at the origin with radius of
convergence larger than $1/3$:
$$f(z)=\sum_{n=2}^\infty f_n z^n.$$
Then interchange of summations provides the identity
\begin{equation}\label{e6}
\sum_r f({1\over r})=\sum_{n=2}^\infty f_n R(n),
\end{equation}
which expresses a sum over primes $r\equiv 3 \pmod 4$ in terms of
values of the base function $R(s)$ at the positive integers
that is computable from~(\ref{e4},\ref{e5a}):
\[R(s)=\frac{1}{2}\sum\frac{\widehat \mu(n)}{n}\log
\left((1-2^{-s})\frac{\zeta(s)}{L(s)}\right).
\]
The values of $L(s)$ themselves can be obtained either from the relation
to the Hurwitz zeta function,
\[
L(s) = \frac{1}{4^s} (\zeta(s, 1/4)-\zeta(s, 3/4)),
\]
(this is computable directly by the Euler-Maclaurin summation formula)
or by reduction to Riemann zeta values in accordance with
the scheme employed
for Gregory's series and mentioned for Catalan's constant.
In the case of the Landau-Ramanujan constant,
applying the general algorithm leads to a wonderful
formula,
\begin{equation}\label{e7}
\lambda
= \frac{1}{\sqrt{2}}
\prod_{n=1}^\infty \left[\left(1-\frac{1}{2^{2^n}}\right)\, \zeta(2^n)
/L(2^n)\right]^{1/2^{n+1}}\,,
\end{equation}
as results from elementary Moebius function identities.
An alternative direct proof runs as follows.
Let
\[
g(s) = (1-2^{-s})\, \zeta(s)/L(s), \qquad
f(s) = \prod_r \frac{1}{1-r^{-s}}
\]
then it is seen that
\[
g(s) = \frac{f(s)^2}{f(s)}
\]
from which it follows that
\[
\frac{f(s)^2}{f(2^{k+1} s)^{1/2^k}}
=
g(s) \, g(2s)^{1/2}\, g(4s)^{1/4}\cdots g(2^k s)^{1/2^k}
\]
proving the identity along with an estimate of the rate of convergence.
Because of the lacunary character of the expression~(\ref{e7}),
the computation
is extremely fast and it takes only $6\cdot 10^8$ machine cycles
(6 seconds of CPU time of 1996!)
to get 200 digits of $\lambda$ by~(\ref{e7}) in pari/gp,
for instance,
\begin{small}
\[\lambda=
0.76422365358922066299069873125009232811679054139340951472\cdots\,.
\]
\end{small}
% \bibliographystyle{acm}
% \bibliography{/home/nuits/algo1/flajolet/lib/bib/algo}
\begin{thebibliography}{1}
\bibitem{Berlekamp68}
{\sc Berlekamp, E.~R.}
\newblock {\em Algebraic Coding Theory}.
\newblock Mc Graw-Hill, 1968.
\newblock Revised edition, 1984.
\bibitem{Berndt94}
{\sc Berndt, B.~C.}
\newblock {\em Ramanujan's Notebooks, Part {IV}}.
\newblock Springer Verlag, 1994.
\bibitem{Edwards74}
{\sc Edwards, H.~M.}
\newblock {\em {R}iemann's {Z}eta Function}.
\newblock Academic Press, 1974.
\bibitem{Finch95b}
{\sc Finch, S.}
\newblock Favorite mathematical constants.
\newblock Available under World Wide Web at {\tt http:
//www.mathsoft.com/asolve/constant.html}, 1995.
\bibitem{HaSaMc93}
{\sc Hafner, J.~L., Sarnak, P., and McCurley, K.}
\newblock Relatively prime values of polynomials.
\newblock In {\em Contemporary Mathematics\/} (1993), M.~Knopp and
M.~Sheingorn, Eds., vol.~143.
\bibitem{Hardy78}
{\sc Hardy, G.~H.}
\newblock {\em Ramanujan: {T}welve Lectures on Subjects Suggested by his Life
and Work}, third~ed.
\newblock Chelsea Publishing Company, New-York, 1978.
\newblock Reprinted and Corrected from the First Edition, Cambridge, 1940.
\bibitem{HaWr79}
{\sc Hardy, G.~H., and Wright, E.~M.}
\newblock {\em An Introduction to the Theory of Numbers}, fifth~ed.
\newblock Oxford University Press, 1979.
\bibitem{Stieltjes93}
{\sc Stieltjes, T.~J.}
\newblock {\em {\OE}uvres Compl{\`e}tes}.
\newblock Springer Verlag, 1993.
\newblock Edited by Gerrit van Dijk.
\bibitem{Vardi91}
{\sc Vardi, I.}
\newblock {\em Computational Recreations in {M}athematica}.
\newblock Addison Wesley, 1991.
\end{thebibliography}
\end{document}