Open In Colab

Background

In this article, we will examine some interesting methods to calculate Pi value. This is a special article that I prepared to celebrate the Pi day, which is March 14th.

Since Python has a built-in pi value in the math module, we can always use it. Python uses double precision to store the values, and double precision numbers have 53 bits (16 digits) of precision. To understand this more, see the code below.

It looks like that there are 48 digits after the decimal point when we try to print 100 digits, but it may not be accurate after 16 digits! In fact, the value is wrong, because the 16th digit should be 2, yet it is shown as 1 above. To better understand this behavior, we can do a simple experiment below. You will see that the digits in Python representation of 0.1 can be random ones after 16 digits.

If we however know the true value of $\pi$, we can represent it using decimal module instead.

Below I describe several methods that can be used to calculate $\pi$ values. For each method, I first calculate it as floating point value, so the actual value calculted may not be accurate after 16th digit. After that, I will re-do the calculation using the decimal module to get higher precision.

Method 1

Ramanujan is an Indian mathematician. He left over 3000 equations to the mankind. His formula to calculate $\pi$ is below:

$\frac{1}{\pi} = \frac{\sqrt{8}}{9801} \displaystyle\sum_{n-0}^{\infty} \frac{(4n)! * (1103+26390n)}{(n!)^4 396^{4n}} $

This equation looks so alien to me. So let's first use the violent way to validate it before going any further.

What can you say? He is like a bug in human system. Nobody knows where this equation comes from or how to explain it and of course he did not give an explanation because he died so early in his age (32 years old).

But let's just be realistic here. If he was born in a rich European family rather than India, then he is probably a famous mathematician in textbooks (even if he still died at 32 years old). He had no formal training at all, and many of his mathematical inventions were already known oustdie of India in his time.

In his time, there is no computer. So he created all these formula in his head. Nowadays, with computers, people are trying to re-create these fancy equations by letting the computer find it by itself. There is a website for doing that.

Next, we try to do the same calculation using decimal module with a 100 digit accuray

Method 2

Chudnovsky formula is a variant of the Ramanujan formula.

$\frac{1}{\pi} = \frac{1}{53360\sqrt{640320}} \displaystyle\sum_{k=0}^{\infty} (-1)^k \frac{(6k)!}{(k!)^3 (3k)!} \times \frac{13591409+545140134k}{640320^{3k}} $

It feels as alien as Ramanujan, but let's test it in the violent way:

As you will see, the convergence is actually faster than Ramanujan formula. At cycle 7, the modification to $1/\pi$ is already in the 1e-100 range.

Basically, each cycle yield approximately 14 additional effective digits (see the change from 1 to 15 to 29 to 43 to 57 above). This formula was previously used to calculate world record of Pi precision. Yet in Ramanujan formula, only 8 digits are generated for each cycle. I know we are talking about $1/\pi$ rather than $\pi$ directly, but the principles of the precision of inverting a number should be the same.

Below is the more precise version using decimal module:

Method 3: Super Pi

People who play with CPU benchmarking and overclocking should be very familiar with the program. I use it a lot and I consider it the best proxy for computational power in bioinformatics research even though it is not an integer algorithm.

In previous Methods we see that infinite series typically increase the number of correct digits additively in successive terms. However, iterative algorithms generally multiply the number of correct digits at each step, so it is preferred over infinite series.

Based on its wiki page, SuperPi uses Gauss–Legendre algorithm and it calculates the decimal point up to a maximum of 32 million. Perhaps the best explanation that I can find online is given here, which describes many algorithms including the Gauss–Legendre algorithm. Basically it contains the following steps:

1, Initialize: let $a_0=1$, $b_0=1/\sqrt{2}$, $t_0=1/4$, $p_0=1$

2, Iterate: $a_{n+1}=\frac{a_n+b_n}{2}$, $b_{n+1}=\sqrt{a_nb_n}$, $t_{n+1}=t_n-p_n(a_n-a_{n+1})^2$, $p_{n+1}=2p_n$

3, Calculate $\pi=\frac{(a_n+b_n)^2}{4t_n}$

Again, let's use violent way to do test it:

It looks like the convergence already reached limit of double digit precision at cycle 4 above.

We next use the decimal module to calculate a more precise value for $\pi$ below. At cycle 7, the increase becomes zero, I think this is an indication that it already reached 100-digit precision at cycle 7. You can change the precision to 300 and see how results change.

Method 4: Borwein algorithm

This is a more modern algorithm that can achieve Nonic convergence, and because it is iterative algorithm, the correct digits are multiplied (rather than linear increase) with each iteration. In other words, each iteration approximately multiplies the number of correct digits by nine.

I copied the Nonic convergence formula from wiki.

1, Initiliaze: let $a_0=1/3$, $r_0=(\sqrt{3}-1)/2$, $s_0=(1-r_0^3)^{1/3}$

2, Iterate: let $t_{n+1}=1+2r_n$, $u_{n+1}=(9r_n(1+r_n+r_n^2))^{1/3}$, $v_{n+1}=t_{n+1}^2+t_{n+1}u_{n+1}+u_{n+1}^2$, $w_{n+1}=\frac{27(1+s_n+s_n^2)}{v_{n+1}}$, $a_{n+1}=w_{n+1}a_n+3^{2n-3}(1-w_{n+1})$, $s_{n+1}=\frac{(1-r_n)^3}{(t_{n+1}+2u_{n+1})v_{n+1}}$, $r_{n+1}=(1-s_{n+1}^3)^{1/3}$

3, Calculate $a_k$ which converges nonically to $1/\pi$

Again, we do the same calculation using decimal module and see what happens. You can change the prec value yourself to see how results change based on precision. The convergence is very fast as you will see below.

In conclusion, this Borwein algorithm is perhaps among the best ways to calculate $\pi$ values.

I hope you enjoy the Pi Day!