Squeezing $\pi$ from 122 Bits

Recently I had some excess of procrastination potential (exam phase), so I finally got around to a fun idea I had many months ago. The following binary term calculates Pi to infinite precision:

00010100011010000000010
11100111000011001011111
   01111       00000
   01110       01110
   01011       11110
   11010       00000
   10111       10000
  01010        000001110
 01110          0111010

…well, almost, let me explain :)

There’s this fractal called a Sierpiński Carpet (also known as its 3D variant, the Menger sponge). It works by dividing some area into $n$ (typically $n=9$) parts, and repeating this step on every subarea except the one in the center. Doing this infinite times will result in a fractal.

There’s also this concept called Lambda Screen that I’ve created. It can render fractals and other images using a recursive, quad-tree-like, purely functional product type. Specifically,

• a screen is a term $\lambda x.(x\ \textrm{tl}\ \textrm{tr}\ \textrm{bl}\ \textrm{br})$ where the terms tl, tr, bl, and br represent the top-left, top-right, bottom-left, and bottom-right quadrants of the image. Each of these terms can either reduce to another screen, or a pixel.
• a pixel is either $\textrm{w}=\lambda w\ b.w$ (on, white) or $\textrm{b}=\lambda w\ b.b$ (off, black).

Rendering of terms will stop at a certain depth by drawing the unreduced terms as grey in the current screen context. The outermost term also gets an empty screen as additional argument to allow a kind of monadic IO – we ignore this here using $\lambda\_$.

By strategic use of fixed-point recursion, you can write some beautifully elegant terms such as this 51-bit term (in Binary Lambda Calculus (BLC), currently the smallest in a codegolf challenge):

$$\lambda\_.(\textrm{y}\ \lambda f\ x.(x\ f\ \textrm{b}\ f\ f))$$

How does this relate to $\pi$? Well, there’s this researcher Ed Pegg who posted something fascinating in a blog post (archive) around 9 years ago. When you increase $n$ by two per “layer” of the Sierpiński Carpet, the total area of the fractal will eventually approach $\frac{\pi}4$! (no, not factorial)

This might seem counterintuitive at first, but it actually makes perfect sense. When splitting a screen and recursing it on itself in each of its $n$ sub-screens, then the total recursive area will represent a variation of the Wallis product

$$\frac{8}{9}\cdot\frac{24}{25}\cdot\frac{48}{49}\cdot\frac{120}{121}\cdot\dots=\frac{\pi}4.$$

So, if I can write a minimal lambda term that generates such fractal, we will only have to measure the recursive (grey) area that the Lambda Screen renderer draws and multiply it by $4$ in order to get an approximation of $\pi$. The only problem is that Lambda Screen does not support arbitrary area splits of $n\ne4$. This was a simple change in the renderer though that doesn’t even affect the behavior of existing quad-tree terms:

• A screen is now a term $\lambda x.(x\ S_{1,1}\ \dots\ S_{1,j}\ S_{2,1}\ \dots\ S_{2, j}\ \dots\ S_{i,1}\ \dots\ S_{i, j})$, where $i=j>0$, i.e. the amount of squares per screen $i\cdot j$ is always a perfect square.
• A pixel remains as before.

This allows us to draw the following $n=9$ Sierpiński Carpet:

$$\lambda\_.(\textrm{y}\ \lambda f\ x.(x\ f\ f\ f\ f\ \textrm{b}\ f\ f\ f\ f))$$

Now there are three remaining open problems:

1. How do we generate $n$ squares from a numeric encoding?
2. How do we exclude the square in the center from the recursion?
3. And how do we pack all of this into a teeny-tiny lambda term?!

The solution to the first problem is quite elegant: A Church numeral is a right-associative $n$-fold composition, while a (single-color) screen is a left-associative $n$-fold composition applied to some color. So, in order to convert between the two, you can left fold the Church numeral itself like so:

$$\lambda f.(\lambda x.(\overbrace{x\circ\dots\circ x}^{\mathclap{n}})\ \lambda x.(x\ f))\leadsto\lambda f\ x.(x\ \overbrace{f\ \dots\ f}^{\mathclap{n}})$$ and thus creating a screen of arbitrary size becomes trivial: $$\lambda f.(\lambda x.(\overbrace{x\circ\dots\circ x}^{\mathclap{n}})\ \lambda x.(x\ f)\ \textrm{b})\leadsto\lambda x.(x\ \overbrace{\textrm{b}\ \dots\ \textrm{b}}^{\mathclap{n}})$$ In our case, the screen will of course contain itself recursively instead of a single color, but the difference doesn’t matter here – it’s still the same argument.

$$* * *$$

Since taking the center element of such a functional structure is difficult and ugly, the second problem requires a little trick: We don’t actually care that the non-recursive term in each layer is precisely at the center – as long as it’s always exactly a single term. So, instead we map the first element of the structure we previously generated to a fixed color while keeping the rest.

We first eta-expand the screen term:

$$\lambda f\ x.(x\ f\ \dots\ f)\htmlClass{katex-flip-x}{\leadsto}\lambda f'\ x'.(\lambda f\ x.(x\ f\ \dots\ f)\ f'\ x')$$

Then the first element applied to $x$ can be mapped to $\textrm{b}$ by injecting a rewrite via $x'$:

$$\lambda f'\ x'.(\lambda f\ x.(x\ f\ \dots\ f)\ f'\ \lambda h\ t.(x'\ \textrm{b}\ t))$$

where we ignore the first argument $h$ and replace it while we keep the rest of the left-associative application (the other screens) as-is (argument $t$).

We also need to implement the actual Church numeral logic, such that we always generate $(n+2)^2$ sub-screens when splitting a layer with $n$ sub-screens. We can do this by giving the initial term an argument $n$, incrementing it two times at every layer, and – before we convert the Church numeral to a screen with $n$ sub-screens using the technique shown above – multiply $n$ by itself.

All the $n$ $f$ above will need to be mapped to the term itself – only that it will now take the incremented Church numeral as argument. All of this comes together in a single term:

$$\lambda\_.(\textrm{y}\ \lambda f\ n\ x.(\lambda f'.(\underbrace{\textrm{2}\ n}_{n^2}\ \lambda x.(x\ f'))\ (f\ \underbrace{(\textrm{inc}\ (\textrm{inc}\ n))}_{\mathrm{next}\ n})\ \lambda h\ t.(x\ \textrm{b}\ t))\ \underbrace{\textrm{3}}_{\mathclap{\mathrm{initial}\ n}})$$

$$* * *$$

Now onto the final problem – making the term tiny! This is done using incredibly complex and highly confidential techniques (i.e. carefully beta-reducing the non-recursive parts with some off-the-shelf strong reducer) and results in the following final term (in bruijn notation):

π [[0 0] [[[1 (1 [0 (3 3 [[1 (1 (4 1 0))]])]) [[2 [[0]] 0]]]]] (+3u)]

Converting the term to BLC will result in the 122-bit program from the beginning. Unfortunately I was not able to golf the Church three into the inner term, so the initialization number alone takes $17.5\%$ of the size!

Here’s the rendering of the term (this may crash your tab if you’re using a phone, chrome, safari, or something else that I don’t use)

While this doesn’t look like an incremented Sierpiński Carpet anymore, it’s actually fully equivalent to a quadrant of Ed Pegg’s design – just with the stable, non-recursive area at the top left.

$$* * *$$

So you probably also want to know how to get a concrete approximation of $\pi$ from this and how quickly it converges.

The simplest way is to count the area that you’re drawing black within the renderer, calculate $(1-\mathrm{area})$ to get the recursive area (grey), and multiply this number by four.¹ This way you will get the following approximations, based on how many layers deep you go:

As you can see, this converges quite slowly. Notably, I’ve never gotten deeper than layer 5 using my JavaScript implementation of Lambda Screen, so you probably won’t get some new world record approximation using this technique. You can see the measurement of the rendering above in your browser’s console, it will probably be around 3.302 (layer 4).

I’ve also implemented a Haskell higher-order reducer for this specific type of problem which counts the total area drawn for any term that uses a perfect square amount of sub-screens. I did some experiments with parallel reduction, but it also does not get much further than the JavaScript renderer.

One problem is that at every level we have to “lower” down the HOAS encoding in order to count the area and continue on all sub-screens – I think this might impact sharing of reductions.

This also leads me to a challenge I have for you! Please let me know if you come up with a reducer that measures the area of such fractals more efficiently. Due to its inherent parallelism and sharing of redexes I believe that optimal reduction could be the optimal approach to this problem. My experiments with the HVM optimal runtime failed though, since the term constantly duplicates non-affine terms.

Also let me know if you find smaller encodings of the fractal (and π in general!), I’m fairly certain that it could be golfed below 100 bit.

You can find the Haskell program in a GitHub Gist – have fun!

$$❦$$

Contact me via email. Support on Ko-fi. Subscribe on RSS. Follow on Mastodon.