In Mathematica, a core principle is that everything should be scalable. So in my job of creating algorithms for Mathematica I have to make sure that everything I produce is scalable.

Last week I decided to test this on one particular example. The problem I chose happens to be a classic. In fact, the very first nontrivial computer program ever written—by Ada Lovelace in 1842—was solving the same problem.

The problem is to compute Bernoulli numbers.

Bernoulli numbers have a long history, dating back at least to Jakob Bernoulli’s 1713 Ars Conjectandi.

Bernoulli’s specific problem was to find formulas for sums like .
Before Bernoulli, people had just made tables of results for specific n and m. But in a Mathematica-like way, Bernoulli pointed out that there was an algorithm that could automate this.

For any given n, the answer is a polynomial in m, and the coefficients are constants that Bernoulli showed could be computed by a simple recurrence formula.

It could have been that Bernoulli numbers would be useful only for solving this particular problem. But in fact over the past 300 years they have found their way into a remarkable range of areas of mathematics. They appear in the formula for the Riemann zeta function at even integers. They are in the coefficients in the series expansion for tan(x). They appear in the Euler-Maclaurin formula for approximating integrals by sums and in the Stirling series for the gamma function. They even relate to Fermat’s last theorem—which was first proved by Ernst Kummer for “regular primes” characterized by Bernoulli numbers.

In 1713 Bernoulli was rather proud of being able to compute the first ten Bernoulli numbers in “a quarter of an hour”.

Of course, in Mathematica, it’s now instantaneous:

But just how far can we go—295 years after Bernoulli, and now with Mathematica?

I decided to try scaling up Bernoulli’s computation—by a factor of a million—and computing the 10 millionth Bernoulli number.

Until recently, doing this would have been impractical, even in Mathematica. Because Mathematica used essentially the same classic recurrence relation for computing Bernoulli numbers that Jakob Bernoulli himself used, and that Ada Lovelace described programming on Charles Babbage’s Analytical Engine.

This algorithm has the feature (already recognized by Lovelace) that it takes about n^2 steps to compute the nth Bernoulli number. So even if one could compute the 10th Bernoulli number in a millisecond, it’d take several thousand years to compute the 10 millionth Bernoulli number.

But a few years ago I programmed a quite different algorithm into Mathematica. Instead of directly computing the Bernoulli numbers using a recurrence relation, I instead used a trick recently suggested by Bernd Kellner: computing Bernoulli numbers by computing the Riemann zeta function.

It’s the integrated nature of Mathematica that makes things like this practical. Without Mathematica, one has to use the simplest building blocks to make efficient algorithms. But with Mathematica, one can take for granted access to efficient very-high-level operations—like computing Riemann zeta functions.

Bernoulli numbers are related to the zeta function by:

The denominator of a Bernoulli number can be computed using a corollary to the von Staudt-Clausen theorem as:

To get the numerator, one then just has to use the relation to the zeta function. The right-hand side is, then, evaluated approximately and multiplied by the already-known denominator of the Bernoulli number. Provided the approximation is done with enough significant digits, the numerator of the Bernoulli number is a mere integer part of the product.

But there’s still a problem: to get all the digits in a Bernoulli number, one has to compute powers of pi to extreme precision.

Of course, Mathematica can do that—to billions of digits if necessary.

A week ago, I took our latest development version of Mathematica, and I typed BernoulliB[10^7].

And then I waited.

Yesterday—5 days, 23 hours, 51 minutes and 37 seconds later—I got the result! [Download 24MB bzip2 file]

The denominator is simple; it’s just 9601480183016524970884020224910.

But the numerator is not so simple; in fact it’s 57,675,292 digits long. It took less than 1 gigabyte of memory to compute, and required computing pi to about 66 million digits.

The numerator is negative, it begins with -47845869850733698144899338333210878162030638218660 and ends with 57164275665935124168181176013725629647185402960697.

Its digits seem almost random. I counted occurrences of all possible k-digit subsequences in the numerator of the 10 millionth Bernoulli number for k=1, 2, 3 and 4. The ratio of the standard deviation to the mean stays low, though grows with k.

So how can we tell if it’s correct?

Bernoulli numbers have lots of interesting properties. One that’s particularly revealing was discovered by Kummer in 1843, and goes by the name of p-adic continuity. In particular, for two integers n and m, and a prime number p such that neither n nor m is divisible by p-1, and a natural number r, such that:

the p-adic continuity asserts that:

This property is completely independent of the way we computed the Bernoulli number.

So let’s start checking. Obviously n=10^7. I checked the p-adic continuity for p=43, r=3, m=59776; p=59, r=2, m=916; p=7919, r=1, m=7484; and p=27449, r=1, m=8928. Every one of these congruences is 0, as it should be.

We’ve successfully found the 10 millionth Bernoulli number.

We’ve done what Ada Lovelace believed should be possible: we’ve mechanized the computation of Bernoulli numbers—so well, in fact, that 295 years after Jakob Bernoulli, its taken us only 500 times longer to compute a million-times-as-large Bernoulli number.

And all with a single line of Mathematica input.