Write a code for estimating Pi up to 15 digits in under 1 second!
Introduction
Currently, the record for the best estimate of Pi approximates it with a precision of 100 trillion decimals, which is about 99 trillion 999 billion 999 million 999 thousand 960 more decimals than what we need for any physical or practical computation in the universe.
This question asks for a lot less! While the record holder had 157 days to run the algorithm, you only have one second.
Another limiting factor for this question is the language you use to write the simulation. C++ would be the obvious choice, with its native speed. But this would be no fun or challenging! So let us try this in Python 🐍 !
Preparation
Let’s start by creating a testing environment; it will consist of two key items: a way of checking the time, and a way of verifying the correctness of our estimation.
For the former, we will rely on the signal‘s alarm function to interrupt our process at the 1-second mark.
For the latter, for simplicity, we will simply hardcode our pi digits. A potentially better way would be to use a slower algorithm, but this is not in our current scope.
Putting those together, we arrive at:
import signal
import math
from decimal import *
from IPython.display import Markdown
# Create a timeout handler & incorporate it in the signal
def timeout_handler(num, stack):
raise Exception("TIMEOUT")
signal.signal(signal.SIGALRM, timeout_handler)
# Hardcoded first 50 pi digits; could be derived via slower algorithm
pi_str = "3.14159265358979323846264338327950288419716939937510"
def check_result(pi_est, it):
# We pad our estimation with zeroes
pi_est = str(pi_est).ljust(100, '0')
wrong_digits = []
result_formatted = ""
# We look through the first 23 digits
# Storing bad digits + preparing text
for i in range(23):
digit_ok = True if pi_est[i] == pi_str[i] else False
text_color = "green" if digit_ok else "red"
result_formatted+=f'{pi_est[i]}'
if digit_ok==False:
wrong_digits.append(i-2)
# Preparing final output
output_md = f"Time is up! Ran {it} iterations."
output_md += f"Result: {result_formatted}"
output_md += f"Correct Digits: {min(wrong_digits)}"
# Returning final output as markdown
return Markdown(output_md)
Leveraging the above, we can use the following structure:
def pi_estimation():
it = 0
try:
while True:
it+=1
[...]
except Exception as ex:
if "TIMEOUT" in ex.args:
pi_est = [...]
display(check_result(pi_est, it))
signal.alarm(1)
pi_estimation()Simulations
Monte Carlo Simulation
When we think of pi, we almost automatically think of circles.
Similarly, computational estimates are closely linked to Monte Carlo methods. As such, the first possible solution should come by combining them. We use Monte Carlo Methods when we have a problem with a probabilistic interpretation.
What probability can we use that involves the value of pi?
Draw the unit circle, and circumscribe it in the square determined by the points (1, 1), (-1, 1), (-1,-1), and (1, -1). The area of the circle, given the radius of one, is pi; the size of the square is 4. For simulated uniform points inside the square, the probability that they are also in the circle is pi/4.
We have the probability needed to run a Monte Carlo estimate, so we generate pairs of coordinates between -1 and 1 until the time is up. Then, compute the ratio of points inside the circle, multiply it by 4, and Voila! We have an estimate for pi.
import numpy as np
# Monte carlo classic
def pi_estimation():
it = 0
in_circ = 0
try:
while True:
it+=1
x = np.random.uniform(-1,1)
y = np.random.uniform(-1,1)
in_circ += int(x**2+y**2<1)
except Exception as ex:
if "TIMEOUT" in ex.args:
pi_est = 4*in_circ/it
display(check_result(pi_est, it))
signal.alarm(1)
pi_estimation()Time is up! Ran 291763 iterations.
Result: 3.142578051363606600000
Correct Digits: 2
In our case, while we were able to run a bit less than 300 thousand simulations, the estimate is not better than what anyone on the street that has heard of geometry could tell you.
But, there is a significant bottleneck in this function. We are generating one value at each step when we could be able to do a lot more. We reduce the number of iterations through the while loop, which improves the overall speed. For example, we could now generate 1000 values at once, increasing the number of simulated points more than 100 times.
import numpy as np
# Monte carlo, faster
def pi_estimation():
it = 0
n = 1000
in_circ = 0
try:
while True:
it+=n
x = np.random.uniform(-1,1,n)
y = np.random.uniform(-1,1,n)
in_circ += np.sum((x**2 + y**2)<1)
except Exception as ex:
if "TIMEOUT" in ex.args:
pi_est = 4*in_circ/it
display(check_result(pi_est, it))
signal.alarm(1)
pi_estimation()Time is up! Ran 33117000 iterations.
Result: 3.141646767521212500000
Correct Digits: 3
Still, this does not produce a meaningful improvement in the estimation, increasing the accurate digits only by one.
We could generate more values at a time, but looking at the results so far, it is improbable that we will get to 15 digits in under 1 second if we start on this road.
Infinite series:
Monte Carlo might seem like a perfect algorithm for computing this estimate, but we can go even 500 years in the past and get a formula, that coupled with the power of modern computers, can give us a better estimate.
Some examples include:
Depending on which one of the above or another you might remember (for example, the more modern Ramanujan formula), you can estimate the value of pi. We are going to look now at the last two, but most of them will give you similar results in terms of speed and accuracy.
Let’s start with Viete.
#Viete
def pi_estimation():
it = 0
s=0
p=1
try:
while True:
it+=1
s = math.sqrt(2+s)
p = p*s/2
except Exception as ex:
if "TIMEOUT" in ex.args:
pi_est = 2/p
display(check_result(pi_est, it))
signal.alarm(1)
pi_estimation()Time is up! Ran 7822541 iterations.
Result: 3.141592653589794000000
Correct Digits: 14
The basic implementation of Viete manages to compute more than 5 million values in the series, which results in 14 accurate digits of pi.
One issue at this point is Python’s accuracy of floats. A possible improvement, bringing both speeds and correctly rounded decimal floating point arithmetic, is to use Decimals. The new algorithm can run 50% more iterations and adds one correct digit. Under the scope of this question, this is enough!
# Viete + decimal
def pi_estimation():
it = 0
s=Decimal(0)
p=Decimal(1)
try:
while True:
it+=1
s = Decimal(math.sqrt(2+s))
p = p*s/Decimal(2)
except Exception as ex:
if "TIMEOUT" in ex.args:
pi_est = str(2/p)
display(check_result(pi_est, it))
signal.alarm(1)
pi_estimation()Time is up! Ran 1126299 iterations.
Result: 3.141592653589793435774
Correct Digits: 15
Given that the use of decimal provided a low-effort, no-cost improvement to our computation, test it with another formula, Nilakantha.
# Nilakantha + Decimal
def pi_estimation():
it = 0
s = Decimal(0)
m = Decimal(1)
try:
while True:
it+=1
x=Decimal(2)*it
s += Decimal(m)/(x*(x+1)*(x+2))
m=m*-1
except Exception as ex:
if "TIMEOUT" in ex.args:
pi_est = (s*4+3)
display(check_result(pi_est, it))
signal.alarm(1)
pi_estimation()Time is up! Ran 1124401 iterations.
Result: 3.141592653589793238638
Correct Digits: 18
Running a similar number of iterations, we get a slightly better estimation, nailing the missing sixteenth decimal in the Viete computation.
Conclusion
With the imposed limit of one second, we’ve succeeded in estimating pi with the required precision. One possible improvement is to use a convergence accelerator, such as a Shanks, that can transform the slowly converging Leibniz series into a much faster accurate approximation.
One thing to remember is that the limit was regarding run time, not complexity. Unfortunately, the dominant factor is the machine that you are running it on. One slightly better comparison is between the number of iterations and the accuracy you obtain, but even that is not able to fully present the differences between algorithms.