Consider a square with an inscribed circle of radius r which itself has an inscribed square rotated 45 degrees. If we look at just the first quadrant, then the area of the outermost square A_outer is r². The area of the circle A_circle is πr²/4, and the area of the inner triangle A_triangle is r²/2. Therefore, we can make two equations for π based on the ratios of these areas: π = 4A_circle/A_outer and π = 2A_circle/A_triangle. Notice how neither equation depends on the radius; so, if we can figure out an alternate way to calculate these areas, we can determine the value of π. Python to the rescue!

Imagine that you’re throwing darts at this collection of shapes, and you’re guaranteed that they’ll all land in the upper right section (a stretch depending on the quality of your dart game, but bear with me). Some of the darts will land within the inscribed circle, and some of those darts will also land in the inner triangle. By adding up the numbers of darts that land in each section, we can approximate the areas of each of these sections.

For simplicity, let r = 1. Since we’re restricting ourselves to the first quadrant, the (x, y) coordinates of each dart will be positive real numbers on the interval [0, 1]. We’ve constrained our system so that every dart thrown lands in the outer square, thus A_outer ≈ N, where N is the total number of darts thrown.

For a dart to land inside the circle, its coordinates must fulfill the equation x² + y² ≤ 1 (remembering that r = 1). Darts that land inside in the triangle have coordinates such that x + y ≤ 1 (where we’ve constrained x and y to be positive). To simulate a dart throw, we can randomly generate two floating-point numbers between 0 and 1 inclusive for the x- and y-coordinates respectively. Then, if they satisfy the circle equation, we increment the number of darts that hit the circle N_circle by 1, and if they satisfy the triangle equation, we increment the number of darts that hit the triangle N_triangle by 1.

Because A_circle ≈ N_circle and A_triangle ≈ N_triangle, we can say that π ≈ 4N_circle/N and π ≈ 2N_circle/N_triangle. Now let’s code up this simulation!

import random
import sys

def point_in_circle(x, y):
	return True if x**2 + y**2 <= 1 else False

def point_in_triangle(x, y):
	return True if x + y <= 1 else False # assume x and y are positive

def get_point():
	return random.random(), random.random()

def get_point_in_shapes():
	x, y = get_point()
	in_circle = point_in_circle(x, y)
	in_triangle = False
	if in_circle: # if the point isn't in the circle, we know that it can't be in the triangle
		in_triangle = point_in_triangle(x, y)
	return in_circle, in_triangle

def get_points_in_shapes(N_darts):
	N_circle, N_triangle = 0, 0
	i = 0 # keep track of how many darts we've thrown so far
	while i < N_darts:
		in_circle, in_triangle = get_point_in_shapes()
		if in_circle:
			N_circle += 1
		if in_triangle:
			N_triangle += 1
		i += 1
	return N_circle, N_triangle

def pi_formula(N_darts, N_circle):
	return 4.0 * N_circle / N_darts

def pi_formula2(N_darts, N_circle, N_triangle):
	if 0 == N_triangle: # if no darts landed in the triangle, avoid dividing by zero by falling
			    # back on the equation that doesn't use the area of the triangle
		return pi_formula(N_darts, N_circle)
	# otherwise, average the two approximations
	return N_circle * (2.0 / N_darts + 1.0 / N_triangle)

def calculate_pi(N_darts):
	N_circle, N_triangle = get_points_in_shapes(N_darts)
	return pi_formula2(N_darts, N_circle, N_triangle)

# driver
N_darts = int(sys.argv[1]) # the total number of darts to throw (given by a command-line argument)
print calculate_pi(N_darts)

(Note how each function is responsible for doing one specific task. This will simplify unit testing, a topic I’ll discuss in a later post.)

While fun to think about, this method takes a while to produce decent results. For example, I can consistently lock in the first three digits of π by “throwing” 10,000,000 darts on my computer, but it takes about 12.3 seconds to run. Additionally, the ability of this method to give a good approximation of π is limited by the precision of Python’s floating-point numbers. The darts’ positions are effectively quantized because Python only uses 53 bits to represent floats. In my next post, I’ll present a variation on this technique that allows us to take advantage of Python’s arbitrarily large integers to improve the approximation.