We can approximate π using nothing more than random numbers and some simple geometry: we randomly throw darts at a square board of side r; within the square we inscribe a quadrant of a circle of radius r with its centre at (0, 0). We count all of the 'throws'; if a dart lands within the quadrant, we also count a 'hit'.
For a large number of throws, we see that:
Some half-remembered geometry tells us that:
I first solved this problem as an undergraduate sometime in 1994 as part of a Computational Physics module. Using FORTRAN.
#!/usr/bin/env python import random import math # class representing a single throw of a dart class Throw: def __init__(self): # generate two random coordinates and work out how far away we are from # the origin self.x = random.random() self.y = random.random() self.distance = self.distance() def distance(self): # the distance from the origin is the hypotenuse of a right-angled # triangle with sides of length and x and y. Pythagoras told us that: # distance = sqrt((x^2) + (y^2)) # which looks like this in python return math.sqrt(self.x**2 + self.y**2) # did we land inside the quadrant? def is_a_hit(self): return self.distance <= 1.0 # main class class MonteCarlo: def __init__(self): # initialise everything self.hits = 0 self.throws = 0 self.pi = 0 # this method is called on every throw def increment(self, throw): # we always clock up another throw self.throws += 1 # and accumulate a hit if we scored if throw.is_a_hit(): self.hits += 1 # then get a new value of pi self.calculate_pi() # explanation can be found here: http://icanhaz.com/montecarlo def calculate_pi(self): self.pi = 4 * (float(self.hits) / float(self.throws)) # we use this in determining whether to print our status def divides_by(self, number): return float(self.throws) % number == 0 # represent the current state as a string def __repr__(self): return "Throws: %10d, Hits: %10d, Pi: %10f, Actual Pi: %10f" % (self.throws, self.hits, self.pi, math.pi) # if we're called on the command line if __name__ == '__main__': import sys try: step = int(sys.argv) except IndexError: step = 1000 # construct a new instance m = MonteCarlo() # loop forever while 1: # keep throwing darts m.increment(Throw()) # only print on every nth iteration if m.divides_by(step): print m
There's now a ruby version, too.
Published under a Creative Commons Attribution-NonCommercial licence