Random Bites of Pi(e)

Wed 14 March 2018 by Moshe Zadka

In today's edition of Pi day post, we will imagine we have a pie. (If you lack imagination, go out and get a pie.) (Even if you do not lack imagination, go out and get a pie.)

As is traditional, we got a round pie. Since pies are important, we will base our unit of measure on this pie -- the diameter of the pie will be 1.

Since we had to carry it home, we put it in a box. We are all ecologically minded, of course, so we put it in a box which is square -- with its length size 1.

We note something interesting -- the box's bottom's area is 1x1 -- or 1. The radius of the pie is 1/2, so its area is pi * 0.25.

As we are driving home, the pie on our passenger seat, we start wondering how we can estimate Pi. Luckily, we got some sugar. What if we sprinkled some sugar, and took notes for each grain, whether it was on the pie, or not?

Let's use Python to simulate:

import random
import attr

First, we need some randomness generator. Then, we also want to use the attrs library, because it makes everything more fun.

We make a Point class. Other than the basics, we give it a class method -- a named constructor which will generate a random point on the unit square -- this is where our sugar grain falls.

We also give it a way to calculate distance from another point, using the Pythagorean theorem.

@attr.s(frozen=True)
class Point:
    x = attr.ib()
    y = attr.ib()

    def distance(self, pt):
        return ((self.x - pt.x) ** 2 + (self.y - pt.y) ** 2) ** 0.5

    @classmethod
    def unit_square_random(cls):
        return cls(x=random.random(), y=random.random())

The center of the pie is at 0.5 by 0.5.

center = Point(0.5, 0.5)

A point is inside the pie if it is less than 0.5 away from the center.

def is_in_circle(pt):
    return center.distance(pt) < 0.5

Now we are ready. Even with just 100,000 grains of sugar, we get 2-digit accuracy.

inside = total = 0
for _ in range(10 ** 5):
    total += 1
    inside += int(is_in_circle(Point.unit_square_random()))
print((inside / total) * 4)
3.14052