Sharedwww / talks / bernoulli / current.texOpen in CoCalc
Author: William A. Stein
1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2%  (c)  William Stein
3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4
5\documentclass{beamer}
6\usepackage{beamerthemesplit}
7\usepackage{graphicx}
8\newcommand{\page}[1]{\frame{#1}}
9%\voffset=-0.1\textheight
10
11\def\Z{{\mathbb{Z}}}
12\def\Q{{\mathbb{Q}}}
13\def\N{{\mathbb{N}}}
14\def\F{{\mathbb{F}}}
15\def\R{{\mathbb{R}}}
16\def\C{{\mathbb{C}}}
17\newcommand{\eps}{\varepsilon}
18
19\usepackage{amsmath}
20\usepackage{amsfonts}
21\usepackage{amssymb}
22\usepackage{amsthm}
23\usepackage{graphicx}
24\usepackage{pstricks}
25\usepackage{color}
26\definecolor{dbluecolor}{rgb}{0,0,0.6}
27\definecolor{dredcolor}{rgb}{.5,0.0,0.0}
28\definecolor{dgreencolor}{rgb}{0,0.4,0}
29\definecolor{blue}{rgb}{.02, .02, .908}
30\newcommand{\dred}{\color{dredcolor}\bf}
31\newcommand{\dblue}{\color{dbluecolor}\bf}
32\newcommand{\dgreen}{\color{dgreencolor}\bf}
33%\newcommand{\dred}{\bf}
34%\newcommand{\dblue}{\bf}
35%\newcommand{\dgreen}{\bf}
36
37\usepackage{xspace}  % to allow for text macros that don't eat space
38\newcommand{\SAGE}{{\color{blue}\sf SAGE}\xspace}
39\newcommand{\sage}{\SAGE}
40
41\newcommand{\graph}[1]{\includegraphics[width=0.45\textwidth]{#1}}
42
43\title{Computing Bernoulli Numbers}
44\author{William Stein}
45\date{(joint work with Kevin McGown of UCSD)\vspace{4ex}\\February 16, 2006}
46\bibliographystyle{amsalpha}
47%\newcommand{\Q}{\mathbb{Q}}
48%\newcommand{\Z}{\mathbb{Z}}
49\newcommand{\defn}[1]{{\em #1}}
50% ---- SHA ----
51\DeclareFontEncoding{OT2}{}{} % to enable usage of cyrillic fonts
52  \newcommand{\textcyr}[1]{%
53    {\fontencoding{OT2}\fontfamily{wncyr}\fontseries{m}\fontshape{n}%
54     \selectfont #1}}
55\newcommand{\Sha}{{\mbox{\textcyr{Sh}}}}
56\DeclareMathOperator{\rank}{rank}
57\DeclareMathOperator{\tor}{tor}
58\DeclareMathOperator{\denom}{denom}
59\DeclareMathOperator{\cond}{cond}
60\DeclareMathOperator{\HH}{H}
61\DeclareMathOperator{\Ker}{Ker}
62\DeclareMathOperator{\ord}{ord}
63\renewcommand{\H}{\HH}
64\newcommand{\hra}{\hookrightarrow}
65\newcommand{\isom}{\cong}
66\newcommand{\ds}{\displaystyle}
67\renewcommand{\t}[1]{\begin{center}{\dblue \large\bf #1}\end{center}}
68
69%%%% Theoremstyles
70\theoremstyle{plain}
71\newtheorem{proposition}{Proposition}
72\newtheorem{claim}{Claim}
73\newtheorem{hypothesis}{Hypothesis}
74\newtheorem{conjecture}{Conjecture}
75
76\theoremstyle{definition}
77\newtheorem{question}[theorem]{Question}
78\newtheorem{alg}[theorem]{Algorithm}
79\newtheorem{openproblem}[theorem]{Open Problem}
80
81%\theoremstyle{remark}
82\newtheorem{goal}[theorem]{Goal}
83\newtheorem{remark}[theorem]{Remark}
84\newtheorem{remarks}[theorem]{Remarks}
85\newtheorem{exercise}[theorem]{Exercise}
86
87\begin{document}
88
89\page{
90\maketitle
91}
92
93\begin{frame}[fragile]
94\t{Bernoulli Numbers}
95Defined by Jacques Bernoulli in posthumous work {\em Ars conjectandi Bale, 1713}.
96
97$$98 \frac{x}{e^x-1}=\sum_{n=0}^\infty \frac{B_n}{n!}\,x^n 99$$
100
101\vfill
102
103$$104 B_0=1,\;\;\; 105 B_1=-\frac{1}{2}\,\;\;\; 106 B_2=\frac{1}{6},\;\;\; 107 B_3=0,\;\;\; 108 B_4=-\frac{1}{30},\;\;\; 109$$
110$$111 B_5=0,\;\;\; 112 B_6=\frac{1}{42},\;\;\; 113 B_7=0,\;\;\; 114 B_8=-\frac{1}{30},\;\;\; 115 B_9=0,\;\;\; 116$$
117
118
119\end{frame}
120
121\begin{frame}[fragile]
122\t{Connection with Riemann Zeta Function}
123For integers $n\geq 2$ we have
124\begin{align*}
125  \zeta(2n)&=\frac{(-1)^{n+1}(2\pi)^{2n}}{2\cdot (2n)!}\,B_{2n}\\
126  \zeta(1-n) &= -\frac{B_n}{n}
127\end{align*}
128So for $n\geq 2$ even:
129$$130 |B_n| =\frac{2n!}{(2\pi)^n}\,\zeta(n) = \pm \frac{n}{\zeta(1-n)}. 131$$
132
133%(Sign is wrong in second equality for $n=1$.)
134\end{frame}
135
136\begin{frame}[fragile]
137\t{Computing Bernoulli Numbers -- say $B_{500}$}
138\begin{verbatim}
139sage: a = maple('bernoulli(500)')     # Wall time: 1.35
140sage: a = maxima('bern(500)')         # Wall time: 0.81
141sage: a = maxima('burn(500)')         # broken...
142sage: a = magma('Bernoulli(500)')     # Wall time: 0.66
143sage: a = gap('Bernoulli(500)')       # Wall time: 0.53
144sage: a = mathematica('BernoulliB[500]')  #W time: 0.18
145  calcbn (http://www.bernoulli.org)   #      Time: 0.020
146sage: a = gp('bernfrac(500)')         # Wall time: 0.00 ?!
147\end{verbatim}
148\end{frame}
149
150\begin{frame}[fragile]
151\t{Computing Bernoulli Numbers -- say $B_{1000}$}
152\begin{verbatim}
153sage: a = maple('bernoulli(1000)')    # Wall time: 9.27
154sage: a = maxima('bern(1000)')        # Wall time: 5.49
155sage: a = magma('Bernoulli(1000)')    # Wall time: 2.58
156sage: a = gap('Bernoulli(1000)')      # Wall time: 5.92
157sage: a = mathematica('BernoulliB[1000]') #W time: 1.01
158  calcbn (http://www.bernoulli.org)   #      Time: 0.06
159sage: a = gp('bernfrac(1000)')        # Wall time: 0.00?!
160\end{verbatim}
161
162{\small NOTE: Mathematica 5.2 is much faster than Mathematica 5.1
163at computing Bernoulli numbers, and the timing is almost
164identical to PARI (for $n>1000$), though amusingly Mathematica 5.2 is
165{\em slow} for $n<1000$!}
166\end{frame}
167
168\begin{frame}[fragile]
169\t{World Records?}
170
171Largest one ever computed was $B_{5000000}$ by O. Pavlyk, which was
172done in Oct. 8, 2005, and whose numerator has $27332507$ digits.
173Computing $B_{10^7}$ is the next obvious challenge.
174\vspace{3ex}
175
176{\bf Bernoulli numbers are really big!}
177\vfill
178
179{\bf Sloane Sequence A103233:}
180\begin{center}
181\begin{tabular}{|l|l|l|l|l|l|l|l|l|}\hline
182n & 0 & 1 & 2   & 3    & 4    &  5     & 6  & 7 \\\hline
183a(n) & 1 & 1 & 83 & 1779 & 27691 & 376772 & 4767554 & ???\\ \hline
184\end{tabular}
185\end{center}
186\vfill
187Here
188$a(n) =$ Number of digits of numerator of $B_{10^n}$.
189\end{frame}
190
191\begin{frame}[fragile]
192\t{Number of Digits}
193Clausen and von Staudt: $\ds d_n = \denom(B_n) = \prod_{p-1\mid m} p$.
194
195Number of digits of numerator is
196$$\lceil \log_{10}(d_n \cdot |B_n|) \rceil$$
197But
198\begin{align*}
199\log(|B_n|) &= \log\left(\frac{2n!}{(2\pi)^n}\,\zeta(n)\right)\\
200     &= \log(2) + \sum_{m=1}^{n} \log(m) - \log(2) - n\log(\pi) + \log(\zeta(n)),
201   \end{align*}
202and $\zeta(n)\sim 1$.   In 10 minutes this gives {\em two
203new entries} for Sloane's sequence:
204   $$205 a(10^7) = 57675292 \qquad\text{and}\qquad a(10^8) = 676752609. 206$$
207\end{frame}
208
209\begin{frame}[fragile]
210\t{Stark's Observation (after talk)}
211Use Stirling's formula, which, ammusingly, involves
212small Bernoulli numbers:
213$$214 \log(\Gamma(z)) = 215 \frac{1}{\log(2\pi)} + \left(z - \frac{1}{2}\right)\log(z) 216 - z + \sum_{n=1}^{\infty} 217 \frac{B_{2n}}{2n(2n-1)z^{2n-1}}. 218$$
219This would make computation of the number of digits
220of the numerator of $B_n$ pretty easy.
221See
222\url{http://mathworld.wolfram.com/StirlingsSeries.html}
223\end{frame}
224
225\begin{frame}[fragile]
226\t{Tables?}
227
228I couldn't find {\em any} interesting tables at all!
229
230\vfill
231
232But from \url{http://mathworld.wolfram.com/BernoulliNumber.html}
233
234"The only known Bernoulli numbers $B_n$ having prime numerators occur for
235n=10, 12, 14, 16, 18, 36, and 42 (Sloane's A092132) [...] with
236no other primes for $n\leq 55274$ (E. W. Weisstein, Apr. 17, 2005)."
237
238\vfill
239
240So maybe $55274$ is the biggest enumeration of $B_k$'s ever? Not
241anymore... since I just used SAGE to script a bunch of PARI's on my
242new 64GB 16-core computer, and made a table of $B_k$ for {\dred $k\leq 24394000$}.  It's very compressed but takes over 3.4GB.
244\end{frame}
245
246\begin{frame}[fragile]
247\t{Buhler et al.}
248Basically, compute $B_k \pmod{p}$ for all $k\leq p$
249and $p$ up to $16\cdot 10^6$ using clever Newton
250iteration to find $1/(e^x-1)$.   In particular,
251if $g$ is an approximation to $f^{-1}$ then ...
252$h=2g-fg^2$'' is twice as good.  (They also gain
253a little using other tricks.)
254\end{frame}
255
256
257\begin{frame}[fragile]
258\t{Math 168 Student Project}
259
260{\em Figure out why PARI is vastly faster than anything else at
261computing $B_k$ and explain it to me.}
262
263{\dred Kevin McGown} rose to the challenge.
264
265{\tiny
266\begin{verbatim}
267/* assume n even > 0. Faster than standard bernfrac for n >= 6 */
268GEN
269bernfrac_using_zeta(long n)
270{
271  pari_sp av = avma;
272  GEN iz, a, d, D = divisors(utoipos( n/2 ));
273  long i, prec, l = lg(D);
274  double t, u;
275
276  d = utoipos(6); /* 2 * 3 */
277  for (i = 2; i < l; i++) /* skip 1 */
278  { /* Clausen - von Staudt */
279    ulong p = 2*itou(gel(D,i)) + 1;
280    if (isprime(utoipos(p))) d = muliu(d, p);
281  }
282  /* 1.712086 = ??? */
283  t = log( gtodouble(d) ) + (n + 0.5) * log(n) - n*(1+log2PI) + 1.712086;
284  u = t / (LOG2*BITS_IN_LONG); prec = (long)ceil(u);
285  prec += 3;
286  iz = inv_szeta_euler(n, t, prec);
287  a = roundr( mulir(d, bernreal_using_zeta(n, iz, prec)) );
288  return gerepilecopy(av, mkfrac(a, d));
289}
290\end{verbatim}
291}
292
293\end{frame}
294
295
296\begin{frame}[fragile]
297\t{Compute $1/\zeta(n)$ to VERY high precision}
298
299{\tiny
300\begin{verbatim}
301/* 1/zeta(n) using Euler product. Assume n > 0.
302 * if (lba != 0) it is log(bit_accuracy) we _really_ require */
303GEN
304inv_szeta_euler(long n, double lba, long prec)
305{
306  GEN z, res = cgetr(prec);
307  pari_sp av0 = avma;
308  byteptr d =  diffptr + 2;
309  double A = n / (LOG2*BITS_IN_LONG), D;
310  long p, lim;
311
312  if (!lba) lba = bit_accuracy_mul(prec, LOG2);
313  D = exp((lba - log(n-1)) / (n-1));
314  lim = 1 + (long)ceil(D);
315  maxprime_check((ulong)lim);
316
317  prec++;
318  z = gsub(gen_1, real2n(-n, prec));
319  for (p = 3; p <= lim;)
320  {
321    long l = prec + 1 - (long)floor(A * log(p));
322    GEN h;
323
324    if (l < 3)         l = 3;
325    else if (l > prec) l = prec;
326    h = divrr(z, rpowuu((ulong)p, (ulong)n, l));
327    z = subrr(z, h);
329  }
330  affrr(z, res); avma = av0; return res;
331}
332\end{verbatim}
333}
334\end{frame}
335
336\begin{frame}[fragile]
337\t{What Does PARI Do?}
338Use
339$$340 |B_n| =\frac{2n!}{(2\pi)^n}\,\zeta(n) 341$$
342and {\em tightly bound} precisions needed to compute each
343quantity.
344{\small
345\begin{verbatim}
346>  (1) Do you know who came up with or implemented the idea
347>      in PARI for computing Bernoulli numbers quickly by
348>      approximating the zeta function and using Classen
349>      and von Staudt's identification of the denominator
350>      of the Bernoulli number?
351
352Henri did, and wrote the initial implementation.
353I wrote the current one (same idea, faster details).
354
355The idea independently came up (Bill Daly) on pari-dev
356as a speed up to Euler-Mac Laurin formulae for zeta or
357gamma/loggamma (that specific one has not been tested/
358implemented so far).
359\end{verbatim}
360}
361\end{frame}
362
363\begin{frame}
364\t{http://www.bernoulli.org/}
365
366Bernd C. Kellner's program at \url{http://www.bernoulli.org/}
367(2002-2004)
368also appears to uses
369$$370 |B_n| =\frac{2n!}{(2\pi)^n}\,\zeta(n) 371$$
372but Kellner's program is closed source and noticeably slower  than PARI
373(2.2.10.alpha).  He claims his program calculates Bernoulli numbers up to
374index $n=10^6$ extremely quickly.''
375\vfill
376Also: {\bf Maxima's} documentation claims to have a function
377{\tt burn} that uses zeta, but it doesn't work (for me).
378\vfill
379%Mathematica 5.2 claims to be able to  but every Mathematica
380%I have access to is version 5.1
381\end{frame}
382
383\begin{frame}[fragile]
384\t{Kevin McGown Project}
385
386\textbf{The Algorithm:}
387Suppose $n\geq 2$ is even.
388\begin{enumerate}
389  \item
390  $\displaystyle 391 K:=\frac{2n!}{(2\pi)^n} 392$
393  \item
394  $\displaystyle 395 d:=\prod_{p-1 | n} p 396$
397  \item
398  $\displaystyle 399 N := \left\lceil (Kd)^{1/(n-1)} \right\rceil 400$
401  \item
402  $\displaystyle 403 z := \prod_{p\leq N} (1-p^{-n})^{-1} 404$
405  \item
406  $\displaystyle 407 a := (-1)^{n/2+1}\,\lceil dKz\rceil 408$
409  \item
410  $\displaystyle 411 B_n = \frac{a}{d} 412$
413\end{enumerate}
414
415
416\end{frame}
417
418\begin{frame}[fragile]
419\t{What About Generalized Bernoulli Numbers?}
420\vfill
421\begin{verbatim}
422>     (2) Has a generalization to generalized
423>         Bernoulli numbers attached to an integer
424>         and Dirichlet character been written
425>         down or implemented?
426
427Not to my knowledge.
428
429Cheers,
430    Karim.
431\end{verbatim}
432\vfill
433\end{frame}
434
435
436
437\begin{frame}[fragile]
438\t{Generalized Bernoulli Numbers}
439Defined in 1958 by H.\thinspace{}W. Leopoldt.
440$$441 \sum_{r=1}^{f-1} \chi(r)\;\frac{t e^{rt}}{e^{ft}-1} 442 =\sum_{n=0}^\infty B_{n,\chi}\;\frac{t^n}{n!} 443$$
444Here $\chi:(\Z/m\Z)\to \C$ is a Dirichlet character.
445
446\vfill
447These give {\bf values at negative integers} of
448associated Dirichlet $L$-functions:
449$$450 L(1-n,\chi) = -\frac{B_{n,\chi}}{n} 451$$
452
453\vfill
454Kubota-Leopoldt $p$-adic $L$-function ($p$-adic interpolation)...
455%$$L_p(1-n, \chi) = -\frac{B_n,{\chi}}{n} (1-\chi(p) p^{n-1}).$$
456
457%where~$\omega$ is the Teichmuller character.
458
459\end{frame}
460
461\begin{frame}[fragile]
462\t{$B_{n,\psi}$ Very Important to Computing Modular Forms}
463
464$$465 E_{k,\chi,\psi}(q) = c_0 + \sum_{m \geq 1} \left( 466 \sum_{n|m} \psi(n) \cdot \chi(m/n) \cdot n^{k-1}\right) q^{m} 467\in \Q(\chi, \psi)[[q]], 468$$
469where
470$$471c_0 = \begin{cases} 0 & \text{ if } L=\cond(\chi) >1, \\ 472\ds- \frac{B_{k,\psi}}{2k} & \text{ if } L=1. 473\end{cases} 474$$
475\begin{theorem}\label{thm:eisgen}
476The (images of) the Eisenstein series above generate
477the Eisenstein subspace $E_k(N,\eps)$, where
478$N=L\cdot \cond(\psi)$ and $\eps=\chi/\psi$.
479\end{theorem}
480
481\end{frame}
482
483\begin{frame}[fragile]
484\t{The Torsion Subgroup of $J_1(p)$}
485
486Let $J_1(p)$ be the Jacobian of the modular curve $X_1(p)$.
487
488\begin{conjecture}[Stein]
489$$\# J_1(p)(\Q)_{\tor} = \frac{p}{2^{p-3}} \cdot\prod_{\chi\neq 1} B_{2,\chi},$$
490where the $\chi$ have modulus $p$.
491(Equivalently, the torsion subgroup is generated by the rational
492cuspidal subgroup---see Kubert-Lang.)
493\end{conjecture}
494(This is a generalization of Ogg's conjecture for $J_0(p)$,
495which Mazur proved.)
496\end{frame}
497
498\begin{frame}[fragile]
499\t{Compute $B_{n,\chi}$? One way.}
500
501Let $N$=modulus of $\chi$, assumed $>1$.
502
503\begin{enumerate}
504\item{} Compute $g = x/(e^{Nx}-1) \in \Q[[x]]$ to
505  precision $O(x^{n+1})$ by computing $e^{Nx}-1 = \sum_{m\geq 1} N^m x^m/m!$
506to precision $O(x^{n+2})$, and computing the inverse
507 $1/(e^{Nx}-1)$, e.g., using Newton iteration as in Buhler et al.
508\item{} For each $a=1,\ldots, N-1$, compute $f_a = g 509 \cdot e^{ax}\in\Q[[x]]$, to precision $O(x^{k+1})$.  This requires
510computing $e^{ax}=\sum_{m\geq 0} a^m x^m/m!$ to precision $O(x^{k+1})$.
511
512\item{} Then for $j\leq n$, we have
513$\displaystyle B_{j,\eps} = j!\cdot 514 \sum_{a=1}^{N-1} \eps(a) \cdot c_j(f_a), 515$ where $c_j(f_a)$ is the
516  coefficient of $x^j$ in~$f_a$.
517\end{enumerate}
518This requires arithmetic {\bf only in $\Q$}, except in the last easy step.
519
520\end{frame}
521
522\begin{frame}[fragile]
523\t{Analytic Method}
524
525Is there an analytic method to compute $B_{n,\chi}$
526that is impressively fast in practice like the one Cohen/Kellner/etc.
527invented for $B_n$?
528\vfill
529
530{\bf YES.}
531\vfill
532\end{frame}
533
534\begin{frame}[fragile]
535\t{Analytic Method}
536
537Assume $\chi$ primitive now.
538
539If $$540 K_{n,\chi}:=(-1)^{n-1}\, 2 n!\left(\frac{N}{2 i}\right)^n 541$$
542then
543$$544 B_{n,\chi}=\frac{K_{n,\chi}}{\pi^n\,\tau(\chi)}\; 545 L(n,\overline{\chi}) 546$$
547
548There is a simple formula for a $d$ such that $d \cdot B_{n,\chi}$ is
549an algebraic integer (analogue of Clausen and von Staudt).
550
551For $n$ large we can compute $L(n,\overline{\chi})$
552{\em very quickly} to high precision; hence  we can compute
553$B_{n,\chi}$ (at least if $\Q(\chi)$ isn't too big,
554e.g., $\Q(\chi)=\Q$ wouldn't be a problem).  (Note, for
555small $n$ that $L(n,\overline{\chi})$ converges slowly; but
556then just use the power series algorithm.)
557
558Compute the conjugates of $d\cdot B_{n,\chi}$ approximately;
559compute minimal polynomial over $\Z$; factor that over $\Q(\chi)$,
560then recognize the right root from the numerical approximation
561to $d\cdot B_{n,\chi}$.
562\end{frame}
563\end{document}
564
565\begin{frame}[fragile]
566\vfill
567\t{Ideas:}
568\begin{enumerate}
569\item Compute $\text{denom} \cdot B_{n,\chi}$ mod $p$ using method
570  similar to Buhler (but with arithmetic as explained in my algorithm
571  above) then use Chinese Remainder Theorem.  Note that we can throw
572  away all but the $n$th coefficient for each $p$, so as to save
573  memory.
574\end{enumerate}
575
576\vfill
577\end{frame}
578
579
580
581
582