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 2*n* bits of precision, or 2*n*/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.

Another option is to use Python’s “decimal” module. The translation is straight-forward, with the main new complexity being “with localcontext() as cxt:” / “cxt.prec = n” instead of “scale = n”, and using a Decimal(“2.5”) and Decimal(i).

Though it looks like Python’s “prec” is more like bc’s “length” than “scale”, so make that “cxt.prec = n+1”.