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 Aouter is r2. The area of the circle Acircle is πr2/4, and the area of the inner triangle Atriangle is r2/2. Therefore, we can make two equations for π based on the ratios of these areas: π = 4Acircle/Aouter and π = 2Acircle/Atriangle. 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 Aouter ≈ N, where N is the total number of darts thrown.
For a dart to land inside the circle, its coordinates must fulfill the equation x2 + y2 ≤ 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 Ncircle by 1, and if they satisfy the triangle equation, we increment the number of darts that hit the triangle Ntriangle by 1.
Because Acircle ≈ Ncircle and Atriangle ≈ Ntriangle, we can say that π ≈ 4Ncircle/N and π ≈ 2Ncircle/Ntriangle. 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) # 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.