The Digits of Pi

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 ${\mathit{foldl}}$ . 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:

$\displaystyle \begin{array}{@{}lcl@{}} \multicolumn{3}{@{}l}{\mathit{stream} :: (\beta \rightarrow \mathsf{Maybe}\;(\gamma,\beta)) \rightarrow (\beta \rightarrow \alpha \rightarrow \beta) \rightarrow \beta \rightarrow [\alpha] \rightarrow [\gamma]} \\ \mathit{stream}\;g\;f\;b\;x &=& \mathbf{case}\;g\;b\;\mathbf{of} \\ & & \quad \begin{array}[t]{@{}lcl@{}} \mathit{Just}\;(c,b') &\rightarrow& c : \mathit{stream}\;g\;f\;b'\;x \\ \mathit{Nothing} &\rightarrow& \mathbf{case}\;x\;\mathbf{of} \\ & & \quad \begin{array}[t]{@{}lcl@{}} a:x' &\rightarrow& \mathit{stream}\;g\;f\;(f\;b\;a)\;x' \end{array} \end{array} \end{array}$

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

$\displaystyle \mathit{stream}\;\mathit{next}\;\mathit{step}\;(0,1)$

where

$\displaystyle \begin{array}{@{}lcl@{}} \mathit{next}\;(u,v) &=& \begin{array}[t]{@{}l} \mathbf{let}\;y = \lfloor{7 \times u \times v}\rfloor\;\mathbf{in} \\ \mathbf{if}\;\lfloor{y}\rfloor = \lfloor{7 \times (u+1) \times v}\rfloor\;\begin{array}[t]{@{}l@{\;}l}\mathbf{then}&\mathit{Just}\;(y,(u - y/(v \times 7), v \times 7))\\\mathbf{else}&\mathit{Nothing} \\ \end{array} \end{array} \\ \mathit{stepl}\;(u,v)\;d &=& (u \times 3 + d, v / 3) \end{array}$

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

Series for pi

Leibniz showed that

$\displaystyle \displaystyle \frac{\pi}{4} = \sum_{i=0}^{\infty} \frac{(-1)^i}{2i+1}$

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

$\displaystyle \pi = \sum_{i=0}^{\infty} \frac{(i!)^2\,2^{i+1}}{(2i+1)!}$

or equivalently

$\displaystyle \pi = 2 + \frac{1}{3} \times \biggl(2 + \frac{2}{5}\times \biggl(2 + \frac{3}{7}\times \biggl( \cdots \biggl(2 + \frac{i}{2i+1}\times \biggl(\cdots\biggr)\biggr)\biggr)\biggr)\biggr)$

This can be seen as the number ${(2;2,2,2...)}$ in a funny mixed-radix base ${(\frac{1}{3}, \frac{2}{5}, \frac{3}{7}...)}$ , just as the usual decimal expansion

$\displaystyle \pi = 3 + \frac{1}{10} \times \biggl(1 + \frac{1}{10}\times \biggl(4 + \frac{1}{10}\times \biggl( \cdots\biggr)\biggr)\biggr)$

is represented by the number ${(3;1,4,1...)}$ in the fixed-radix base ${(\frac{1}{10},\frac{1}{10},\frac{1}{10}...)}$ . Computing the decimal digits of ${\pi}$ 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

$\displaystyle x = \frac{1}{m} \times \biggl(x_0 + \frac{1}{m}\times \biggl(x_1 + \frac{1}{m}\times \biggl( \cdots\biggr)\biggr)\biggr)$

in base ${m}$ , in which ${0 \le x_i < m}$ for each digit ${x_i}$ . We are to convert it to a similar representation

$\displaystyle x = \frac{1}{n} \times \biggl(y_0 + \frac{1}{n}\times \biggl(y_1 + \frac{1}{n}\times \biggl( \cdots\biggr)\biggr)\biggr)$

in base ${n}$ , in which ${0 \le y_j < n}$ for each output digit ${y_j}$ . The streaming process maintains a state ${(u,v)}$ , a pair of rationals; the invariant is that after consuming ${i}$ input digits and producing ${j}$ output digits, we have

$\displaystyle x = \frac{1}{n} \times \biggl(y_0 + \cdots + \frac{1}{n}\times \biggl(y_{j-1} + v \times (u + \frac{1}{m} \times \biggl( x_i + \frac{1}{m} \times \biggl(x_{i+1} + \cdots \biggr)\biggr)\biggr)\biggr)\biggr)$

so that ${(u,v)}$ represents a linear function ${(v\times) \cdot (u+)}$ that should be applied to the value represented by the remaining input.

We can initialize the process with ${i=0, j=0, u=0, v=1}$ . At each step, we first try to produce another output digit. The remaining input digits ${x_i, x_{i+1},...}$ represent a value in the unit interval; so if ${n \times v \times (u+0)}$ and ${n \times v \times (u+1)}$ have the same integer part, then that must be the next output digit, whatever the remaining input digits are. Let ${y_j = \lfloor n \times v \times u \rfloor}$ be that integer. Now we need to find ${(u',v')}$ such that

$\displaystyle \frac{1}{n} \times \biggl(y_j + v' \times (u' + r)\biggr) = v \times (u + r)$

for any remainder ${r}$ ; then we can increment ${j}$ and set ${(u,v)}$ to ${(u',v')}$ and the invariant is maintained. A little algebra shows that we should take ${v' = n \times v}$ and ${u' = u - y_j/v'}$ .

If ${v \times u}$ and ${v \times (u+1)}$ 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 ${(u',v')}$ such that

$\displaystyle v \times \biggl(u + \frac{1}{m} \times \biggl(x_i + r\biggr)\biggr) = v' \times (u' + r)$

for any remainder ${r}$ ; then we can increment ${i}$ and set ${(u,v)}$ to ${(u',v')}$ and the invariant is again maintained. Again, algebraic manipulation leads us to ${v' = v/m}$ and ${u' = m \times u + x_i}$ .

For example, ${\frac{1}{4} = 0.020202..._3 = 0.151515..._7}$ , and the conversion starts as follows:

$\displaystyle \begin{array}{c|c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}cc} x_i & & 0 && 2 && 0 && && 2 && && 0 && 2 \\ \hline (u,v) & \bigl(\frac{0}{1},\frac{1}{1}\bigr) && \bigl(\frac{0}{1},\frac{1}{3}\bigr) && \bigl(\frac{2}{1},\frac{1}{9}\bigr) && \bigl(\frac{6}{1},\frac{1}{27}\bigr) && \bigl(\frac{15}{7},\frac{7}{27}\bigr) && \bigl(\frac{59}{7},\frac{7}{81}\bigr) && \bigl(\frac{8}{49},\frac{49}{81}\bigr) && \bigl(\frac{24}{49},\frac{49}{243}\bigr) && \bigl(\frac{170}{49},\frac{49}{729}\bigr) & \cdots \vrule height 2.5ex depth 1.5ex width 0pt \\ \hline y_j & & && && && 1 && && 5 \end{array}$

That is, the initial state is ${u_0=\frac{0}{1}, v_0=\frac{1}{1}}$ . This state does not yet determine the first output digit, so we consume the first input digit 0 to yield the next state ${u_1 = \frac{0}{1}, v_1 = \frac{1}{3}}$ . 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 ${u_3 = \frac{6}{1}, v_3 = \frac{1}{27}}$ . This state does determine the next digit: ${v_3 \times u_3 = 0.020_3 = 0.136..._7}$ and ${v_3 \times (u_3+1) = 0.021_3 = 0.154..._7}$ both start with a 1 in base 7. So we can produce a 1 as the first output digit, yielding state ${u_4 = \frac{15}{7}, v_4 = \frac{7}{27}}$ . And so on.

The process tends to converge. Each production step widens the non-empty window ${[n \times v \times u, n \times v \times (u+1))}$ by a factor of ${n}$ , so it will eventually contain multiple integers; therefore we cannot produce indefinitely. Each consumption step narrows the window by a factor of ${m}$ , so it will tend towards eventually producing the next output digit. However, this doesn’t always work. For example, consider converting ${0.333..._{10}}$ to base 3:

$\displaystyle \begin{array}{c|c@{}c@{}c@{}c@{}c@{}c@{}cc} x_i & & 3 && 3 && 3 & \\ \hline (u,v) & \bigl(\frac{0}{1},\frac{1}{1}\bigr) && \bigl(\frac{3}{1},\frac{1}{10}\bigr) && \bigl(\frac{33}{1},\frac{1}{100}\bigr) && \bigl(\frac{333}{1},\frac{1}{1000}\bigr) & \cdots \vrule height 2.5ex depth 1.5ex width 0pt \\ \hline y_j & & && \end{array}$

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 ${\pi}$ . We have the input

$\displaystyle \pi = 2 + \frac{1}{3} \times \biggl(2 + \frac{2}{5}\times \biggl(2 + \frac{3}{7}\times \biggl( \cdots \biggl(2 + \frac{i}{2i+1}\times \biggl(\cdots\biggr)\biggr)\biggr)\biggr)\biggr)$

which we want to convert to decimal. The streaming process maintains a pair ${(u,v)}$ of rationals—but this time representing the linear function ${(u+) \cdot (v\times)}$ , since this time our expression starts with a sum rather than a product. The invariant is similar: after consuming ${i-1}$ input digits and producing ${j}$ output digits, we have

$\displaystyle \pi = y_0 + \frac{1}{10} \times \biggl(\cdots y_{j-1} + \frac{1}{10} \times \biggl(u + v \times \biggl(x_i + \frac{i}{2i+1} \times \biggl(x_{i+1} + \frac{i+1}{2i+3} \times \cdots\biggr)\biggr)\biggr)\biggr)$

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

We can initialize the process with ${i=1, j=0, u=0, v=1}$ . At each step, we first try to produce an output digit. What value might the remaining input

$\displaystyle r = 2 + \frac{i}{2i+1} \times \biggl(2 + \frac{i+1}{2i+3} \times \cdots \biggr)$

represent? Each of the bases is at least ${\frac{1}{3}}$ , so it is clear that ${r_{\mathrm{min}} \le r}$ , where

$\displaystyle r_{\mathrm{min}} = 2 + \frac{1}{3} \times r_{\mathrm{min}}$

which has unique solution ${r_{\mathrm{min}} = 3}$ . Similarly, each of the bases is less than ${\frac{1}{2}}$ , so it is clear that ${r < r_{\mathrm{max}}}$ , where

$\displaystyle r_{\mathrm{max}} = 2 + \frac{1}{2} \times r_{\mathrm{max}}$

which has unique solution ${r_{\mathrm{max}} = 4}$ . So we consider the bounds ${\lfloor u + v \times 3 \rfloor}$ and ${\lfloor u + v \times 4 \rfloor}$ ; if these have the same integer part ${y_j}$ , then that is the next output digit. Now we need to find ${(u',v')}$ such that

$\displaystyle y_j + \frac{1}{10} \times (u' + v' \times r) = u + v \times r$

for any remainder ${r}$ , so we pick ${u' = 10 \times (u - y_j)}$ and ${v' = 10 \times v}$ . Then we can increment ${j}$ and set ${(u,v)}$ to ${(u',v')}$ , 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 ${(u',v')}$ such that

$\displaystyle u' + v' \times r = u + v \times \biggl(x_i + \frac{i}{2i+1} \times r\biggr)$

for all ${r}$ , so we pick ${u' = u + v \times x_i}$ and ${v' = v \times i / (2i+1)}$ . Then we can increment ${i}$ and set ${(u,v)}$ to ${(u',v')}$ , and again the invariant is maintained.

The conversion starts as follows:

$\displaystyle \begin{array}{c|c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}c@{}cc} x_i & & 2 && && 2 && 2 && && 2 && 2 \\ \hline (u,v) & \bigl(\frac{0}{1},\frac{1}{1}\bigr) && \bigl(\frac{2}{1},\frac{1}{3}\bigr) && \bigl(\frac{-10}{1},\frac{10}{3}\bigr) && \bigl(\frac{10}{3},\frac{4}{3}\bigr) && \bigl(\frac{-2}{3},\frac{4}{7}\bigr) && \bigl(\frac{-50}{3},\frac{40}{7}\bigr) && \bigl(\frac{-110}{21},\frac{160}{63}\bigr) && \bigl(\frac{-10}{63},\frac{800}{693}\bigr) & \cdots \vrule height 2.5ex depth 1.5ex width 0pt \\ \hline y_j & & && 3 && && && 1 && && \end{array}$

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 ${\mathit{stream}}$ function above:

$\displaystyle \mathit{piDigits} = \mathit{stream}\;g\;f\;(1,0,0\%1,1\%1)\;(\mathit{repeat}\;2)$

where

$\displaystyle g\;(i,j,u,v) = \begin{array}[t]{@{}l} \mathbf{if}\;y == \mathit{floor}\;(u + v \times 4) \\ \mathbf{then}\;\mathit{Just}\;(y, (i,j+1, 10 \times (u - \mathit{fromIntegral}\;y), 10 \times v)) \\ \mathbf{else}\;\mathit{Nothing} \\ \mathbf{where}\;y = \mathit{floor}\;(u + v \times 3) \end{array}$

and

$\displaystyle f\;(i,j,u,v)\;x = \begin{array}[t]{@{}l} (i+1,j,u + v \times \mathit{fromIntegral}\;x, v \times i' / (2 \times i' + 1)) \\ \mathbf{where}\;i' = \mathit{fromIntegral}\;i \end{array}$

(The ${\%}$ s make rational numbers in Haskell, and force the ambiguous fractional type to be ${\mathit{Rational}}$ rather than ${\mathit{Double}}$ .)

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 ${j}$ 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 ${i}$ component as a rational in the first place:

$\displaystyle \mathit{piDigits} = \begin{array}[t]{@{}l} \mathit{go}\;((1,0,1) :: (\mathit{Rational},\mathit{Rational},\mathit{Rational}))\;\mathbf{where} \\ \qquad \mathit{go}\;(i,u,v) = \begin{array}[t]{@{}ll} \mathbf{if} & y == \mathit{floor}\;(u+v \times 4) \\ \mathbf{then} & y : \mathit{go}\;(i,10 \times (u-\mathit{fromIntegral}\;y),10 \times v) \\ \mathbf{else} & \mathit{go}\;(i+1,u+2 \times v, (v \times i) / (2 \times i+1)) \\ \multicolumn{2}{@{}l}{\qquad \mathbf{where}\; y = \mathit{floor}\;(u+v \times 3)} \end{array} \end{array}$

Then we have

$\displaystyle \mathit{piDigits} = [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3...$

4 Responses to The Digits of Pi

Peter Milley says:

Thursday, November 9th, 2017 at 7:44 pm

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.

- jeremygibbons says:
  
  Friday, November 10th, 2017 at 1:29 pm
  
  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.
  
zenzike says:

Friday, November 10th, 2017 at 9:19 am

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

The Digits of Pi

Series for pi

Conversion from a fixed base

Conversion from a mixed base

Code

Related

About jeremygibbons

4 Responses to The Digits of Pi

Leave a Reply Cancel reply

Where to Start

Recent Posts

Archives

Meta

Where to go next?