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}
95
Defined 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}
123
For 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*}
128
So 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}
139
sage: a = maple('bernoulli(500)') # Wall time: 1.35
140
sage: a = maxima('bern(500)') # Wall time: 0.81
141
sage: a = maxima('burn(500)') # broken...
142
sage: a = magma('Bernoulli(500)') # Wall time: 0.66
143
sage: a = gap('Bernoulli(500)') # Wall time: 0.53
144
sage: a = mathematica('BernoulliB[500]') #W time: 0.18
145
calcbn (http://www.bernoulli.org) # Time: 0.020
146
sage: 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}
153
sage: a = maple('bernoulli(1000)') # Wall time: 9.27
154
sage: a = maxima('bern(1000)') # Wall time: 5.49
155
sage: a = magma('Bernoulli(1000)') # Wall time: 2.58
156
sage: a = gap('Bernoulli(1000)') # Wall time: 5.92
157
sage: a = mathematica('BernoulliB[1000]') #W time: 1.01
158
calcbn (http://www.bernoulli.org) # Time: 0.06
159
sage: 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
163
at computing Bernoulli numbers, and the timing is almost
164
identical 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
171
Largest one ever computed was $B_{5000000}$ by O. Pavlyk, which was
172
done in Oct. 8, 2005, and whose numerator has $27332507$ digits.
173
Computing $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
182
n & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\\hline
183
a(n) & 1 & 1 & 83 & 1779 & 27691 & 376772 & 4767554 & ???\\ \hline
184
\end{tabular}
185
\end{center}
186
\vfill
187
Here
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}
193
Clausen and von Staudt: $\ds d_n = \denom(B_n) = \prod_{p-1\mid m} p$.
194
195
Number of digits of numerator is
196
$$\lceil \log_{10}(d_n \cdot |B_n|) \rceil$$
197
But
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*}
202
and $\zeta(n)\sim 1$. In 10 minutes this gives {\em two
203
new 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)}
211
Use Stirling's formula, which, ammusingly, involves
212
small 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
$$
219
This would make computation of the number of digits
220
of the numerator of $B_n$ pretty easy.
221
See
222
\url{http://mathworld.wolfram.com/StirlingsSeries.html}
223
\end{frame}
224
225
\begin{frame}[fragile]
226
\t{Tables?}
227
228
I couldn't find {\em any} interesting tables at all!
229
230
\vfill
231
232
But from \url{http://mathworld.wolfram.com/BernoulliNumber.html}
233
234
"The only known Bernoulli numbers $B_n$ having prime numerators occur for
235
n=10, 12, 14, 16, 18, 36, and 42 (Sloane's A092132) [...] with
236
no other primes for $n\leq 55274$ (E. W. Weisstein, Apr. 17, 2005)."
237
238
\vfill
239
240
So maybe $55274$ is the biggest enumeration of $B_k$'s ever? Not
241
anymore... since I just used SAGE to script a bunch of PARI's on my
242
new 64GB 16-core computer, and made a table of $B_k$ for {\dred $k\leq
243
94000$}. It's very compressed but takes over 3.4GB.
244
\end{frame}
245
246
\begin{frame}[fragile]
247
\t{Buhler et al.}
248
Basically, compute $B_k \pmod{p}$ for all $k\leq p$
249
and $p$ up to $16\cdot 10^6$ using clever Newton
250
iteration to find $1/(e^x-1)$. In particular,
251
``if $g$ is an approximation to $f^{-1}$ then ...
252
$h=2g-fg^2$'' is twice as good. (They also gain
253
a 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
261
computing $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 */
268
GEN
269
bernfrac_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 */
303
GEN
304
inv_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);
328
NEXT_PRIME_VIADIFF(p,d);
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?}
338
Use
339
$$
340
|B_n| =\frac{2n!}{(2\pi)^n}\,\zeta(n)
341
$$
342
and {\em tightly bound} precisions needed to compute each
343
quantity.
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
352
Henri did, and wrote the initial implementation.
353
I wrote the current one (same idea, faster details).
354
355
The idea independently came up (Bill Daly) on pari-dev
356
as a speed up to Euler-Mac Laurin formulae for zeta or
357
gamma/loggamma (that specific one has not been tested/
358
implemented so far).
359
\end{verbatim}
360
}
361
\end{frame}
362
363
\begin{frame}
364
\t{http://www.bernoulli.org/}
365
366
Bernd C. Kellner's program at \url{http://www.bernoulli.org/}
367
(2002-2004)
368
also appears to uses
369
$$
370
|B_n| =\frac{2n!}{(2\pi)^n}\,\zeta(n)
371
$$
372
but 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
374
index $n=10^6$ extremely quickly.''
375
\vfill
376
Also: {\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:}
387
Suppose $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
427
Not to my knowledge.
428
429
Cheers,
430
Karim.
431
\end{verbatim}
432
\vfill
433
\end{frame}
434
435
436
437
\begin{frame}[fragile]
438
\t{Generalized Bernoulli Numbers}
439
Defined 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
$$
444
Here $\chi:(\Z/m\Z)\to \C$ is a Dirichlet character.
445
446
\vfill
447
These give {\bf values at negative integers} of
448
associated Dirichlet $L$-functions:
449
$$
450
L(1-n,\chi) = -\frac{B_{n,\chi}}{n}
451
$$
452
453
\vfill
454
Kubota-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
$$
469
where
470
$$
471
c_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}
476
The (images of) the Eisenstein series above generate
477
the 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
486
Let $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},$$
490
where the $\chi$ have modulus $p$.
491
(Equivalently, the torsion subgroup is generated by the rational
492
cuspidal subgroup---see Kubert-Lang.)
493
\end{conjecture}
494
(This is a generalization of Ogg's conjecture for $J_0(p)$,
495
which Mazur proved.)
496
\end{frame}
497
498
\begin{frame}[fragile]
499
\t{Compute $B_{n,\chi}$? One way.}
500
501
Let $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!$
506
to 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
510
computing $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}
518
This 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
525
Is there an analytic method to compute $B_{n,\chi}$
526
that is impressively fast in practice like the one Cohen/Kellner/etc.
527
invented for $B_n$?
528
\vfill
529
530
{\bf YES.}
531
\vfill
532
\end{frame}
533
534
\begin{frame}[fragile]
535
\t{Analytic Method}
536
537
Assume $\chi$ primitive now.
538
539
If $$
540
K_{n,\chi}:=(-1)^{n-1}\, 2 n!\left(\frac{N}{2 i}\right)^n
541
$$
542
then
543
$$
544
B_{n,\chi}=\frac{K_{n,\chi}}{\pi^n\,\tau(\chi)}\;
545
L(n,\overline{\chi})
546
$$
547
548
There is a simple formula for a $d$ such that $d \cdot B_{n,\chi}$ is
549
an algebraic integer (analogue of Clausen and von Staudt).
550
551
For $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,
554
e.g., $\Q(\chi)=\Q$ wouldn't be a problem). (Note, for
555
small $n$ that $L(n,\overline{\chi})$ converges slowly; but
556
then just use the power series algorithm.)
557
558
Compute the conjugates of $d\cdot B_{n,\chi}$ approximately;
559
compute minimal polynomial over $\Z$; factor that over $\Q(\chi)$,
560
then recognize the right root from the numerical approximation
561
to $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