Computing Bernoulli Numbers

William Stein

(joint work with Kevin McGown of UCSD) April 14, 2006

Bernoulli Numbers

Defined by Jacques Bernoulli in posthumous work Ars conjectandi Bale, 1713.

\frac{x}{e^{x} − 1} = ∑_{n=0}^{∞} \frac{B_{n}}{n!} x^{n}

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,

Connection with Riemann Zeta Function

For integers n ≥ 2 we have

ζ(2n) = \frac{(−1)^{n+1}(2π)^{2n}}{2 · (2n)!} B_{2n}
ζ(1 − n) = −\frac{B_{n}}{n}

So for n ≥ 2 even:

|B_{n}| = \frac{2 · n!}{(2π)^{n}} ζ(n) = ±\frac{n}{ζ(1 − n)}.

Computing Bernoulli Numbers – say B_{500}

sage: a = maple('bernoulli(500)')     # Wall time: 1.35
sage: a = maxima('bern(500)')         # Wall time: 0.81
sage: a = maxima('burn(500)')         # broken...
sage: a = magma('Bernoulli(500)')     # Wall time: 0.66
sage: a = gap('Bernoulli(500)')       # Wall time: 0.53
sage: a = mathematica('BernoulliB[500]')  #W time: 0.18
  calcbn (http://www.bernoulli.org)   #      Time: 0.020
sage: a = gp('bernfrac(500)')         # Wall time: 0.00 ?!

Computing Bernoulli Numbers – say B_{1000}

sage: a = maple('bernoulli(1000)')    # Wall time: 9.27
sage: a = maxima('bern(1000)')        # Wall time: 5.49
sage: a = magma('Bernoulli(1000)')    # Wall time: 2.58
sage: a = gap('Bernoulli(1000)')      # Wall time: 5.92
sage: a = mathematica('BernoulliB[1000]') #W time: 1.01
  calcbn (http://www.bernoulli.org)   #      Time: 0.06
sage: a = gp('bernfrac(1000)')        # Wall time: 0.00?!

NOTE: Mathematica 5.2 is much faster than Mathematica 5.1 at computing Bernoulli numbers; it takes only about twice as long as PARI (for n > 1000), though amusingly Mathematica 5.2 is slow for n ≤ 1000!

World Records?

Largest one ever computed was B_{5000000} by O. Pavlyk, which was done in Oct. 8, 2005, and whose numerator has 27332507 digits. Computing B_{10^{7}} is the next obvious challenge.

Bernoulli numbers are really big!

Sloane Sequence A103233:

n  0 1 2 3  4  5  6  7
a(n) 1 1 83 1779 27691 376772 4767554 ???

Here a(n) = Number of digits of numerator of B_{10^{n}}.

Number of Digits

Clausen and von Staudt: d_{n} = denom(B_{n}) = ∏_{p−1|m}p. Number of digits of numerator is

⌈log_{10}(d_{n} · |B_{n}|)⌉

But

log(|B_{n}|) = log (\frac{2n!}{(2π)^{n}} ζ(n))
= log(2) + ∑_{m=1}^{n}log(m) − n log(2) − n log(π) + log(ζ(n)),

and ζ(n) ∼ 1. This quickly gives new entries for Sloane’s sequence:

a(10^{7}) = 57675292  and  a(10^{8}) = 676752609.
William Stein  Computing Bernoulli Numbers

Stark’s Observation (after talk)

Use Stirling’s formula, which, ammusingly, involves small Bernoulli numbers:

log(Γ(z)) = \frac{1}{log(2π)} + (z − \frac{1}{2})log(z) − z + ∑_{n=1}^{∞} \frac{B_{2n}}{2n(2n − 1)z^{2n−1}}.

This would make computation of the number of digits of the numerator of B_{n} pretty easy. See

http://mathworld.wolfram.com/StirlingsSeries.html

Tables?

I couldn’t find any interesting tables at all!

But from

http://mathworld.wolfram.com/BernoulliNumber.html

”The only known Bernoulli numbers B_{n} having prime numerators occur for n=10, 12, 14, 16, 18, 36, and 42 (Sloane’s A092132) [...] with no other primes for n ≤ 55274 (E. W. Weisstein, Apr. 17, 2005).”

So maybe 55274 is the biggest enumeration of B_{k}’s ever? Not anymore... since I used SAGE to script a bunch of PARI’s on my new 64GB 16-core computer, and made a table of B_{k} for k ≤ 100000. It’s very compressed but takes over 3.4GB (and is “stuck” in that broken computer.)

Buhler et al.

Basically, compute B_{k} (mod p) for all k ≤ p and p up to 16 · 10^{6} using clever Newton iteration to find 1/(e^{x} − 1). In particular, “if g is an approximation to f^{−1} then ... h = 2g − fg^{2}” is twice as good. (They also use a few other tricks.)

Math 168 Student Project

Figure out why PARI is vastly faster than anything else at computing B_{k} and explain it to me. Kevin McGown rose to the challenge.

/* assume n even > 0. Faster than standard bernfrac for n >= 6 */
GEN
bernfrac_using_zeta(long n)
{
pari_sp av = avma;
GEN iz, a, d, D = divisors(utoipos( n/2 ));
long i, prec, l = lg(D);
double t, u;
d = utoipos(6); /* 2 * 3 */
for (i = 2; i < l; i++) /* skip 1 */
{ /* Clausen - von Staudt */
  ulong p = 2*itou(gel(D,i)) + 1;
  if (isprime(utoipos(p))) d = muliu(d, p);
}
/* 1.712086 = ??? */
t = log( gtodouble(d) ) + (n + 0.5) * log(n) - n*(1+log2PI) + 1.712086;
u = t / (LOG2*BITS_IN_LONG); prec = (long)ceil(u);
prec += 3;
iz = inv_szeta_euler(n, t, prec);
a = roundr( mulir(d, bernreal_using_zeta(n, iz, prec)) );
return gerepilecopy(av, mkfrac(a, d));
}

Compute 1/ζ(n) to VERY high precision

/* 1/zeta(n) using Euler product. Assume n > 0.
* if (lba != 0) it is log(bit_accuracy) we _really_require */
GEN
inv_szeta_euler(long n, double lba, long prec)
{
GEN z, res = cgetr(prec);
pari_sp av0 = avma;
byteptr d =  diffptr + 2;
double A = n / (LOG2*BITS_IN_LONG), D;
long p, lim;
if (!lba) lba = bit_accuracy_mul(prec, LOG2);
D = exp((lba - log(n-1)) / (n-1));
lim = 1 + (long)ceil(D);
maxprime_check((ulong)lim);
prec++;
z = gsub(gen_1, real2n(-n, prec));
for (p = 3; p <= lim;)
{
  long l = prec + 1 - (long)floor(A * log(p));
  GEN h;
  if (l < 3)         l = 3;
  else if (l > prec) l = prec;
  h = divrr(z, rpowuu((ulong)p, (ulong)n, l));
  z = subrr(z, h);
  NEXT_PRIME_VIADIFF(p,d);
}
affrr(z, res); avma = av0; return res;
}

What Does PARI Do?

Use

|B_{n}| = \frac{2n!}{(2π)^{n}} ζ(n)

and tightly bound precisions needed to compute each quantity.

>  (1) Do you know who came up with or implemented the idea
>      in PARI for computing Bernoulli numbers quickly by
>      approximating the zeta function and using Classen
>      and von Staudt's identification of the denominator
>      of the Bernoulli number?
Henri did, and wrote the initial implementation.
I wrote the current one (same idea, faster details).
The idea independently came up (Bill Daly) on pari-dev
as a speed up to Euler-Mac Laurin formulae for zeta or
gamma/loggamma (that specific one has not been tested/
implemented so far).

http://www.bernoulli.org/

Bernd C. Kellner’s program at http://www.bernoulli.org/

(2002-2004) also appears to uses

|B_{n}| = \frac{2n!}{(2π)^{n}} ζ(n)

but Kellner’s program is closed source and noticeably slower than PARI (2.2.10.alpha). He claims his program “calculates Bernoulli numbers up to index n = 10^{6} extremely quickly.”

Also: Maxima’s documentation claims to have a function burn that uses zeta, but it doesn’t work (for me).

Kevin McGown Project

The Algorithm: Suppose n ≥ 2 is even.

1. K = \frac{2n!}{(2π)^{n}}

2. d = ∏_{p−1|n}p

3. N = ⌈(Kd)^{1/(n−1)}⌉

4. z = ∏_{p≤N}(1 − p^{−n})^{−1}

5. a = (−1)^{n/2+1} ⌈dKz⌉

6. B_{n} = \frac{a}{d}

What About Generalized Bernoulli Numbers?

>     (2) Has a generalization to generalized
>         Bernoulli numbers attached to an integer
>         and Dirichlet character been written
>         down or implemented?
Not to my knowledge.
Cheers,
    Karim.

Generalized Bernoulli Numbers

Defined in 1958 by H. W. Leopoldt.

∑_{r=1}^{f −1}χ(r) \frac{te^{rt}}{e^{ft} − 1} = ∑_{n=0}^{∞}B_{n,χ} \frac{t^{n}}{n!}

Here χ : (ℤ/mℤ) → ℂ is a Dirichlet character. These give values at negative integers of associated Dirichlet L-functions:

L(1 − n, χ) = −\frac{B_{n,χ}}{n}

Kubota-Leopoldt p-adic L-function (p-adic interpolation)...

William Stein  Computing Bernoulli Numbers

B_{n,ψ} Very Important to Computing Modular Forms

E_{k,χ,ψ}(q) = c_{0}+∑_{m≥1} (∑_{n|m}ψ(n) · χ(m/n) · n^{k−1})q^{m} ∈ ℚ(χ, ψ)[[q]],

where

c_{0} = \begin{cases}0 & if L = cond(χ) > 1, \\ −\frac{B_{k,ψ}}{2k} & if L = 1. \end{cases}

Theorem

The (images of) the Eisenstein series above generate the Eisenstein subspace E_{k}(N, ε), where N = L · cond(ψ) and ε = χ/ψ.

The Torsion Subgroup of J_{1}(p)

Let J_{1}(p) be the Jacobian of the modular curve X_{1}(p).

Conjecture (Stein)

#J_{1}(p)(ℚ)_{tor} = \frac{p}{2^{p−3}} · ∏_{χ≠1} B_{2,χ},

where the χ have modulus p. (Equivalently, the torsion subgroup is generated by the rational cuspidal subgroup—see Kubert-Lang.) (This is a generalization of Ogg’s conjecture for J_{0}(p), which Mazur proved.)

Compute B_{n,χ}? One way.

Let N=modulus of χ, assumed > 1.

1. Compute g = x/(e^{Nx} − 1) ∈ ℚ[[x]] to precision O(x^{n+1}) by computing e^{Nx} − 1 = ∑_{m≥1} N^{m}x^{m}/m! to precision O(x^{n+2}), and computing the inverse 1/(e^{Nx} − 1), e.g., using Newton iteration as in Buhler et al.

2. For each a = 1, …, N − 1, compute f_{a} = g · e^{ax} ∈ ℚ[[x]], to precision O(x^{k+1}). This requires computing e^{ax} = ∑_{m≥0} a^{m}x^{m}/m! to precision O(x^{k+1}).

3. Then for j ≤ n, we have B_{j,ε} = j! · ∑_{a=1}^{N−1}ε(a) · c_{j}(f_{a}), where

c_{j}(f_{a}) is the coefficient of x^{j} in f_{a}.

This requires arithmetic only in ℚ, except in the last easy step.

Analytic Method

Is there an analytic method to compute B_{n,χ} that is impressively fast in practice like the one Cohen/Kellner/etc. invented for B_{n}?

YES.

Analytic Method

Assume χ primitive now. If

K_{n,χ} = (−1)^{n−1} 2n! (\frac{N}{2i} )^{n}

then

B_{n,χ} = \frac{K_{n,χ}}{π^{n} τ(χ)} L(n, \overline{χ})

There is a simple formula for a d such that d · B_{n,χ} is an algebraic integer (analogue of Clausen and von Staudt). For n large we can compute L(n, \overline{χ}) very quickly to high precision; hence we can compute B_{n,χ} (at least if ℚ(χ) isn’t too big, e.g., ℚ(χ) = ℚ wouldn’t be a problem). (Note, for small n that L(n, \overline{χ}) converges slowly; but then just use the power series algorithm.) Compute the conjugates of d · B_{n,χ} approximately; compute minimal polynomial over ℤ; factor that over ℚ(χ), then recognize the right root from the numerical approximation to d · B_{n,χ}.

William Stein  Computing Bernoulli Numbers