\documentclass[12pt]{article}
\renewcommand{\baselinestretch}{1.07}
\usepackage{amsmath}
\usepackage{amssymb}
\DeclareMathOperator{\tor}{tor}
\DeclareMathOperator{\im}{im}
\DeclareMathOperator{\lcm}{lcm}
\DeclareMathOperator{\denom}{denom}
\def\Z{{\mathbb{Z}}}
\def\Q{{\mathbb{Q}}}
\def\N{{\mathbb{N}}}
\def\F{{\mathbb{F}}}
\def\R{{\mathbb{R}}}
\def\C{{\mathbb{C}}}
\def\beqn{\begin{eqnarray*}}
\def\eeqn{\end{eqnarray*}}
\begin{document}
\title{Computing Bernoulli Numbers Quickly}
\author{Kevin J. McGown}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The Bernoulli numbers are defined via the coefficients of
the power series expansion of $t/(e^t-1)$.
Namely, for integers $m\geq 0$ we define $B_m$ so that
$$
\frac{t}{e^t-1}=\sum_{m=0}^\infty \frac{B_m}{m!}\,t^m
.
$$
Multiplying both sides above by $e^t-1$ and equating
coefficients of $t^{m+1}$ yields:
$$
B_0=1\text{,}
\qquad
(m+1)B_m=-\sum_{k=0}^{m-1}\binom{m+1}{k}B_{k}
$$
Some authors take the above recurrence to be the
definition of the Bernoulli numbers.
This recurrence provides a straightforward method
for calculating $B_m$ and is especially convenient
for computing $B_m$ for all $m$ up to some bound.
The first few Bernoulli numbers are:
$$
B_0=1,\;\;\;
B_1=-\frac{1}{2}\,\;\;\;
B_2=\frac{1}{6},\;\;\;
B_3=0,\;\;\;
B_4=-\frac{1}{30},\;\;\;
$$
$$
B_5=0,\;\;\;
B_6=\frac{1}{42},\;\;\;
B_7=0,\;\;\;
B_8=-\frac{1}{30},\;\;\;
B_9=0,\;\;\;
$$
$$
B_{10}=\frac{55}{66}\;\;,\;
B_{11}=0,\;\;\;
B_{12}=-\frac{691}{2730}\;\;
B_{13}=0,\;\;\;
B_{14}=-\frac{7}{6}\;\;
$$
The values above provide evidence for two basic results
regarding Bernoulli numbers.
First, $B_m=0$ for odd $m\geq 3$, and secondly, for even $m\geq 2$,
$B_m=(-1)^{m/2+1}\,|B_m|$.
Henceforth we will denote
$$
B_m=\frac{a}{d}
$$
where $a,d\in\Z$, $d\geq 1$, and $(a,d)=1$.
From the properties mentioned above, it is clear that
$a=(-1)^{m/2+1}\;|a|$
for even $m\geq 2$.
The goal of this paper is to compute $B_m$ rapidly,
when $m$ is potentially very large.
Computing $B_m$ via the recurrence
is slow; it requires us to sum
over $m(m+1)/2$ terms. In addition, this method
requires storing the numbers $B_0,\dots,B_{m-1}$ in memory.
In order to speed up this computation,
we will describe an important connection between
the Bernoulli numbers and the Riemann Zeta Function.
Much of what we will describe was gleaned
from the PARI-2.2.11.alpha source code. The algorithm this version of
PARI uses to compute Bernoulli numbers was written by Henri Cohen and
later refined by Karim Belabas; it was originally designed to
speed up the computation of zeta values.
For real $s> 1$, Euler defined the function
$$
\zeta(s)=\sum_{n=1}^\infty n^{-s}
\,,
$$
for which he proved the product formula
$$
\zeta(s)=\prod_{p} (1-p^{-s})^{-1}
\,.
$$
(Riemann showed that $\zeta(s)$ has an analytic continuation to
the entire complex plane with a simple pole at $s=1$,
and hence the function bears his name.)
In addition to giving the product formula,
Euler was able to evaluate the
zeta function at the even positive integers.~\cite{B:1}
For any integer $m\geq 1$,
$$
2\zeta(2m)=\frac{(-1)^{m+1}(2\pi)^{2m}}{(2m)!}\,B_{2m}
\,.
$$
It follows that for any even integer $m\geq 4$,
$$
|B_m|=\frac{2m!}{(2\pi)^m}\,\zeta(m)
\,.
$$
It is now clear how to compute decimal approximations to $B_m$;
we merely approximate $\zeta(m)$
using the Euler product and plug the result into the above equation.
A priori, this is not enough to compute $B_m$ as
a ratio of integers, but fortunately a theorem of
Clausen and von Staudt precisely describes the denominator of $B_m$
in terms of the divisors of $m$.~\cite{B:1} For even $m\geq 2$,
$$
d:=\denom(B_m)=\prod_{p-1 | m} p
\,.
$$
Now we show how to compute $a$.
First define
$$
K:=\frac{2m!}{(2\pi)^m}
$$
so that $|B_m|=K\zeta(m)$.
Using the Euler product, we may approximate
$\zeta(m)$ from below with arbitrary precision.
Suppose that we have computed a number $z$ such that
$$
0\leq \zeta(m)-z<(Kd)^{-1}
\,;
$$
then we have
$$
0\leq |B_m|-zK1$ and
$\varepsilon>0$, find $N\in\Z^+$ so that
when we set
$$
z:=\prod_{p\leq N} (1-p^{-s})^{-1}
\,,
$$
we are guaranteed that $0\leq\zeta(s)-z<\varepsilon$.
We always have $0\leq\zeta(s)-z$. Further,
it is not hard to see that
$$
\sum_{n\leq N} n^{-s}\leq\prod_{p\leq N} (1-p^{-s})^{-1}
$$
and therefore
\begin{eqnarray*}
\zeta(s)-z
&\leq&
\sum_{n=N+1}^\infty n^{-s}\\
&\leq&
\int_N^\infty x^{-s}\;dx\\
&=&
\frac{1}{(s-1)N^{s-1}}
\,.
\end{eqnarray*}
If we choose $N>\varepsilon^{-1/(s-1)}$, then we have
$$
\frac{1}{(s-1)N^{s-1}}
\leq
\frac{1}{N^{s-1}}
<
\varepsilon
\,,
$$
which implies
$\zeta(s)-z<\varepsilon$, as required.
For our purposes, we have $s=m$ and $\varepsilon=(Kd)^{-1}$
and therefore it suffices to choose
$N>(Kd)^{1/(m-1)}$.
\pagebreak
\vspace{24pt}
\noindent
\textbf{The Algorithm:}
Suppose $m\geq 2$ is even.
\begin{enumerate}
\item
$$
K:=\frac{2m!}{(2\pi)^m}
$$
\item
$$
d:=\prod_{p-1 | m} p
$$
\item
$$
N := \left\lceil (Kd)^{1/(m-1)} \right\rceil
$$
\item
$$
z := \prod_{p\leq N} (1-p^{-m})^{-1}
$$
\item
$$
a := (-1)^{m/2+1}\,\lceil dKz\rceil
$$
\item
$$
B_m = \frac{a}{d}
$$
\end{enumerate}
Some remarks are in order. In step (1), we must be careful to
compute $K$ to sufficient precision so that the calculation
in (5) gives the desired result.
In order to compute (4), it is useful to first compute all primes
$p\leq N$; this may be done quickly using the Sieve of
Eratosthenes. One may also compute the product in (2)
via a sieving process. Finally, for the value of $N$ we may choose
any integer greater than or equal to the one specified in (3), so we
need not worry about computing $(Kd)^{1/(m-1)}$ to much precision.
It is interesting to note that the algorithm above also gives
a way of approximating $\zeta(m)$ quickly for even $m$. Namely, compute
$B_{m}$ as a rational number using this algorithm and plug it into
Euler's formula for $\zeta(m)$ along with an approximation of $\pi$
sufficiently many decimal places.
\pagebreak
\vspace{24pt}
\noindent
\textbf{Example:}
We use the modest size example of $m=50$
for the sake of readability.
Using 50 digits of precision,
we compute
\begin{eqnarray*}
K=7500866746076957704747736.7155247316456403804367604
\,.
\end{eqnarray*}
The divisors of $m$ are ${1, 2, 5, 10, 25, 50}$ and hence
$$
d = (2)(3)(11)=66
\,.
$$
We find $N=4$ and compute
$$
z =
1.0000000000000008881784210930815902983501390146827
\,.
$$
Finally we compute
$$
dKz=
495057205241079648212477524.99999999442615111210652
$$
and therefore
$$
B_{50}=\frac{495057205241079648212477525}{66}
\,.
$$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{1}
\bibitem{B:1}
Ireland, Kenneth. Rosen, Michael.
\emph{A Classical Introduction to Modern Number Theory.}
New York: Springer-Verlag, 1990.
\end{thebibliography}
\end{document}