The previous post discussed how you would plan an attempt to set a record in computing ζ(3), also known as Apéry’s constant. Specifically that post looked at how to choose your algorithm and how to anticipate the number of terms to use.

Now suppose you wanted to actually do the calculation. This post will carry out a calculation using the Unix utility bc. I often use bc for extended precision calculation, but mostly for simple things [1].

Calculating Apéry’s constant will require a function to compute binomial coefficients, something not built into bc, and so this post will illustrate how to write custom functions in bc.

(If you wanted to make a serious attempt at setting a record computing Apéry’s constant, or anything else, you would probably use an extended precision library written in C MPFR or write something even lower level; you would not use bc.)

Apéry’s series

The series we want to evaluate is

To compute this we need to compute the binomial coefficients in the denominator, and to do that we need to compute factorials.

From calculations in the previous post we estimate that summing n terms of this series will give us 2n bits of precision, or 2n/3 decimal places. So if we carry out our calculations to n decimal places, that gives us more precision than we need: truncation error is much greater than numerical error.

bc code

Here’s the code, which I saved in the file zeta3.b. I went for simplicity over efficiency. See [2] for a way to make the code much more efficient.

    # inefficient but simple factorial
    define fact(x) {
        if (x <= 1)
            return (1);
        return (x*fact(x-1));
    }

    # binomial coefficient n choose k
    define binom(n, k) {
        return (fact(n) / (fact(k)*fact(n-k)));
    }

    define term(n) {
        return ((-1)^(n-1)/(n^3 * binom(2*n, n)))
    }

    define zeta3(n) {
        scale = n
        sum = 0
        for (i = 1; i <= n; ++i)
            sum += term(i);
        return (2.5*sum);
    }

Now say we want 100 decimal places of ζ(3). Then we should need to sum about 150 terms of the series above. Let’s sum 160 terms just to be safe. I run the code above as follows [3].

    $ bc -lq zeta3.b
    zeta3(160)

This returns

    1.202056903159594285399738161511449990764986292340498881792271555341\
83820578631309018645587360933525814493279797483589388315482364276014\
64239236562218818483914927

How well did we do?

I tested this by computing ζ(3) to 120 decimals in Mathematica with

    N[RiemannZeta[3], 120]

and subtracting the value returned by bc. This shows that the error in our calculation above is approximately 10^-102. We wanted at least 100 decimal places of precision and we got 102.

Notes

[1] I like bc because it’s simple. It’s a little too simple, but given that almost all software errs on the side of being too complicated, I’m OK with bc being a little too simple. See this post where I used (coined?) the phrase controversially simple.

[2] Not only is the recursive implementation of factorial inefficient, computing factorials from scratch each time, even by a more efficient algorithm, is not optimal. The more efficient thing to do is compute each new coefficient by starting with the previous one. For example, once we’ve already computed the binomial coefficient (200, 100), then we can multiply by 202*201/101² in order to get the binomial coefficient (202, 101).

Along the same lines, computing (-1)^(n-1) is wasteful. When bootstrapping each binomial coefficient from the previous one, multiply by -1 as well.

[3] Why the options -lq? The -l option does two things: it loads the math library and it sets the precision variable scale to 20. I always use the -l flag because the default scale is zero, which lops off the decimal part of all floating point numbers. Strange default behavior! I also often need the math library. Turns out -l wasn’t needed here because we explicitly set scale in the function zeta3, and we don’t use the math library.

I also use the -q flag out of habit. It starts bc in quiet mode, suppressing the version and copyright announcement.

Before carrying out a big calculation, you want to have an idea how long various approaches would take. This post will illustrate this by planning an attempt to calculate Apéry’s constant

to enormous precision. This constant has been computed to many decimal places, in part because it’s an open question whether the number has a closed form. Maybe you’d like to set a new record for computing it.

This post was motivated by a question someone asked me in response to this post: What the advantage is in having an approximation for binomial coefficients when we can just calculate them exactly? This post will show these approximations in action.

Apéry’s series

First of all, calculating ζ(3) directly from the definition above is not the way to go. The series converges too slowly for that. If you compute the sum of the first N terms, you get about 2 log₁₀ N decimal places. This post explains why. We can do much better, with the number of decimal places achieved being proportional to N rather than log N.

A couple days ago I wrote about a sequence rapidly converging series for ζ(3), starting with Apéry’s series.

Suppose you want to compute ζ(3) to a million decimal places using this series. How would you do this?

Since this is an alternating series, the truncation error from summing only the first N terms is bounded by the N+1 term. So you could sum terms until you get a term less than

10^-106.

Great, but how long might that take?

Without some estimate of binomial coefficients, there’s no way to know how long this would take without computing a bunch of binomial coefficients. But with a little theory, you can do a back-of-the-envelope estimate.

Estimating how many terms to use

For large n, the binomial coefficient in the denominator of the nth term is approximately 4ⁿ/√(πn). So when you include the n³ term, you can see that the denominator of the nth term in the series is greater than 4ⁿ = 2²ⁿ. In other words, each term gives you at least 2 bits of precision.

If you want 1,000,000 digits, that’s approximately 3,000,000 bits. So to get a million digits of ζ(3) you’d need to sum about 1,500,000 terms of Apéry’s series. By comparison, summing 1,500,000 terms in the definition of ζ(3) would give you between 12 and 13 digits precision.

More efficient series

Apéry’s series is only the first of a sequence of rapidly converging series for ζ(3). The series get more complicated but more efficient as a parameter s increases. Is it worthwhile to use a larger value of s? If you were planning to devote, say, a month of computing resources to try to set a record, it would be a good idea to spend a few hours estimating whether another approach would get more work done in that month.

By combining information here and here you could find out that the series corresponding to s = 2 has terms whose denominator is on the order of 3³ⁿ. To estimate how many terms you’d need for a million digits of precision, you’d need to solve the equation

3³ⁿ = 10¹⁰⁶,

and by taking logs you find this is 1,047,951. In other words, about 1/3 fewer terms than using Apéry’s series. That was a win, so maybe you go on to estimate whether you should go up to s = 3. The key information you need in these calculations is the estimate on binomial coefficients given here.

Incidentally, computing ζ(3) to a million decimal places would have given you the world record in 1997. At the moment, the record is on the order of a trillion digits.

Update: The next post shows how to actually carry out (on a smaller scale) the calculations described here using the Unix utility bc.

I ran across an interesting footnote in Wolfram Koepf’s book Computer Algebra.

Gosper’s algorithm [1] was probably the first algorithm which would not have been found without computer algebra. Gosper writes in his paper: “Without the support of MACSYMA and its developer, I could not have collected the experiences necessary to provoke the conjectures that lead to the algorithm.”

When I first glanced at the footnote I misread it, coming away with the impression that Gosper said he could not have proved that his algorithm was correct without computer algebra. But he said he would not have formulated his algorithm without his experience using MACSYMA.

If a sum of hypergeometric terms is itself hypergeometric, Gosper’s algorithm gives a way to find the sum. Gosper’s algorithm, discovered using computer algebra, is itself an important computer algebra algorithm.

Related posts

• Brief notes on hypergeometric functions (4-page PDF)
• Hypergeometric distribution and series
• Hypergeometric functions are key

[1] Koepf cites Gosper as R. W. Gosper, Jr. Decision procedure for indefinite hypergeometric summation. Proc. Natl. Acad. Sci. USA, 75, 1978, 40–42.