In the previous post we were introduced to *metamorphisms*, which consist of an unfold after a fold—typically on lists, and the fold part typically a . A canonical example is the conversion of a fraction from one base to another. For simplicity, let’s consider here only infinite fractions, so we don’t have to deal with the end of the input and flushing the state:

So for example, we can convert an infinite fraction in base 3 to one in base 7 with

where

In this post, we’ll see another number conversion problem, which will deliver the digits of . For more details, see my paper—although the presentation here is simpler now.

## Series for pi

Leibniz showed that

From this, using Euler’s convergence-accelerating transformation, one may derive

or equivalently

This can be seen as the number in a funny mixed-radix base , just as the usual decimal expansion

is represented by the number in the fixed-radix base . Computing the decimal digits of is then a matter of conversion from the mixed-radix base to the fixed-radix base.

## Conversion from a fixed base

Let’s remind ourselves of how it should work, using a simpler example: conversion from one fixed base to another. We are given an infinite-precision fraction in the unit interval

in base , in which for each digit . We are to convert it to a similar representation

in base , in which for each output digit . The streaming process maintains a state , a pair of rationals; the invariant is that after consuming input digits and producing output digits, we have

so that represents a linear function that should be applied to the value represented by the remaining input.

We can initialize the process with . At each step, we first try to produce another output digit. The remaining input digits represent a value in the unit interval; so if and have the same integer part, then that must be the next output digit, whatever the remaining input digits are. Let be that integer. Now we need to find such that

for any remainder ; then we can increment and set to and the invariant is maintained. A little algebra shows that we should take and .

If and have different integer parts, we cannot yet tell what the next output digit should be, so we must consume the next input digit instead. Now we need to find such that

for any remainder ; then we can increment and set to and the invariant is again maintained. Again, algebraic manipulation leads us to and .

For example, , and the conversion starts as follows:

That is, the initial state is . This state does not yet determine the first output digit, so we consume the first input digit 0 to yield the next state . This state still does not determine the first output, and nor will the next; so we consume the next two input digits 2 and 0, yielding state . This state does determine the next digit: and both start with a 1 in base 7. So we can produce a 1 as the first output digit, yielding state . And so on.

The process tends to converge. Each production step widens the non-empty window by a factor of , so it will eventually contain multiple integers; therefore we cannot produce indefinitely. Each consumption step narrows the window by a factor of , so it will tend towards eventually producing the next output digit. However, this doesn’t always work. For example, consider converting to base 3:

The first output digit is never determined: if the first non-3 in the input is less than 3, the value is less than a third, and the first output digit should be a 0; if the first non-3 is greater than 3, then the value is definitely greater than a third, and it is safe to produce a 1 as the first output digit; but because the input is all 3s, we never get to make this decision. This problem will happen whenever the value being represented has a finite representation in the output base.

## Conversion from a mixed base

Let’s return now to computing the digits of . We have the input

which we want to convert to decimal. The streaming process maintains a pair of rationals—but this time representing the linear function , since this time our expression starts with a sum rather than a product. The invariant is similar: after consuming input digits and producing output digits, we have

Note that the output base is fixed at 10; but more importantly, the input *digits* are all fixed at 2, and it is the input *base* that varies from digit to digit.

We can initialize the process with . At each step, we first try to produce an output digit. What value might the remaining input

represent? Each of the bases is at least , so it is clear that , where

which has unique solution . Similarly, each of the bases is less than , so it is clear that , where

which has unique solution . So we consider the bounds and ; if these have the same integer part , then that is the next output digit. Now we need to find such that

for any remainder , so we pick and . Then we can increment and set to , and the invariant is maintained.

If the two bounds have different integer parts, we must consume the next input digit instead. Now we need to find such that

for all , so we pick and . Then we can increment and set to , and again the invariant is maintained.

The conversion starts as follows:

Happily, non-termination ceases to be a problem: the value being represented does not have a finite representation in the output base, being irrational.

## Code

We can plug these definitions straight into the function above:

where

and

(The s make rational numbers in Haskell, and force the ambiguous fractional type to be rather than .)

In fact, this program can be considerably simplified, by inlining the definitions. In particular, the input digits are all 2, so we need not supply them. Moreover, the component of the state is never used, because we treat each output digit in the same way (in contrast to the input digits); so that may be eliminated. Finally, we can eliminate some of the numeric coercions if we represent the component as a rational in the first place:

Then we have

Weirdly this didn’t work correctly until I compiled it. When I tried it just in ghci I got errors after 16 digits:

> take 20 $ piDigits

[3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,1,3,3,4]

But once I compiled it, the error went away.

Oops, yes. The numeric type of the non-integer calculations is ambiguous, and defaults to Double unless you do something different. That’s what the type declaration for

gois for in the more concise program. I’ve added two percents in the longer program, which have the same effect.This is one of my favourite programs! Glad to see a blog post about it 🙂

Pingback: Resumen de lecturas compartidas durante noviembre de 2017 | Vestigium