Pi Calculating Programs

Quick links: none at the moment — see below for some code snippets!

For a long time, one of my side interests has been writing computer programs to calculate the value of the mathematical constant $\pi$. In middle school, I wrote Java implementations of the well-known Chudnovsky algorithm to calculate lots of digits of $\pi$:

\[\frac{1}{\pi} = 12 \sum^\infty_{q=0} \frac{(-1)^q (6q)! (545140134q + 13591409)}{(3q)! (q!)^3 (640320)^{3q+\frac{3}{2}}}\]

and of my favorite such algorithm, Gauss-Legendre:

\[\begin{align} a_0 &= 1 & a_{n+1} &= \frac{a_n+b_n}{2}\\ b_0 &= \frac{1}{\sqrt{2}} & b_{n+1} &= \sqrt{a_n b_n}\\ t_0 &= \frac{1}{4} & t_{n+1} &= t_n-p_n(a_n-a_{n+1})^2\\ p_0 &= 1 & p_{n+1} &= 2p_n\\ & & \pi_{n+1} &= \frac{(a_{n+1}+b_{n+1})^2}{4t_{n+1}}\text{.} \end{align}\]

Later, I became interested in goals more exotic than calculating as many digits as possible. The Monte Carlo method imagines a circle inscribed in a square; it has one generate random points within the square and test using the distance formula whether they fall within the circle. From this, the ratio of the area of the circle to the area of the square can be approximated, and anyone who understands the formula for the area of a circle can do a little algebra to calculate $\pi$ from that. Thus, while nowhere near as efficient as the above formulas, the Monte Carlo method is much more explainable, a virtue in its own right. While experimenting with a Java implementation of this algorithm, I took on the secondary goal of making my code as short as possible — an objective known as “code golf” — and ended up with this:

for(long t=0,n=0,y=1;;n-=y*y+(y=(int)(t*y+n>>7))*y>>62)System.out.println(n++*4d/t++);

It may look nothing like an implementation of the Monte Carlo method to calculate $\pi$, but it is one (see the video at the top of the page). Every step taken from the description of the method to this code can be fully explained to someone without a deep math or programming background (hints: the t*y+n>>7 part functions as a custom pseudorandom number generator and the >>62 uses long integer overflow to test whether the point is in the circle).

I was interested in how fast I could make this algorithm run, and in learning to write code in assembly, so I rewrote it in C and then directly in x86_64. My assembly version (blog post, GitHub repo), a fairly direct translation, performs better than GCC’s most optimized version, and I’m not quite sure why.

Keeping the code golf goal but switching back to digits over explainability, I manipulated a neat little formula of Simon Plouffe’s:

\[\pi+3 = \sum^\infty_{n=1} \frac{2^n (n!)^2 n}{(2n)!}\]

to produce this, in Python:

b=i=j=k=1;a=-3
while 1:
 j*=i;k*=2*i-1;s=a*k+i*j*b;t=b*k;i+=1
 if a/b==s/t:print(s//t,end='',flush=1);j*=10;s=s%t*10
 while b:a,b=b,a%b
 b=t//a;a=s//a

which takes advantage of the fact that Python integers can get arbitrarily large to endlessly produce digits of $\pi$ as they become available (see the video at the top of the page). In fact, the last two lines merely reduce the giant fraction at the heart of the expression to speed up the computation; the algorithm still works without them:

t=i=j=k=1;s=-3
while 1:
 d=s/t;j*=i;k*=2*i-1;s=s*k+i*j*t;t*=k;i+=1
 if d==s/t:print(s//t,end='',flush=1);j*=10;s=s%t*10

In the future, I hope to produce more thorough explanations of how these algorithms work and to experiment more with so-called “digit extraction algorithms” that let one calculate a particular digit of $\pi$ without calculating all the digits before it.