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 go is 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