I was recently reading the book Fooled by Randomness, and the author made an interesting observation: You can compute Pi via Monte Carlo simulation.

The intuition is: imagine a circle inscribed tanged to the inside of a square. Now choose a point at random in the figure, so that the point is guaranteed to be inside the square, but may or may not be inside the circle. The probability that the point is inside the circle is the ratio of the area of the circle to the total area. And the ratio of the area of the total circle to the area of the square has Pi in it. So if we run this simulation a lot of times, and see how many points are inside the circle vs outside, we can get probabilistically close to Pi. (There are analytical ways of counting Pi, but this way is easier to code).

Concretely:

1. Let there be a unit circle inscribed on a Cartesian coordinate plane at the origin.
2. Let there be a square of side length 1 such that one corner is at the origin, one corner is a (0, 1), one is at (1,0), and one is at (1,1).
3. Let a be the area in the first quadrant inside the circle.
4. Let b be the area in the first quadrant inside the square and outside the circle
5. We can see that a + b =1, because the 1 by 1 square has the area of 1.
6. We can see that a is the area of the unit circle in the first quadrant, that is a = PI / 4
7. That means that PI / 4 + b = 1, which means that PI = 4 - 4b
8. b is approximated by the proportion of random points in the square but not in the circle.

That’s enough to make a good guess at Pi. Here’s some code!

use rand::prelude::*;

struct Point {
    x: f64,
    y: f64,
}

fn main() {
    let mut rng = rand::thread_rng();

    let iterations = 10_000_000;

    let b_count = (0..iterations)
        .map(|_| random_point(&mut rng))
        .filter(|point| !is_in_circle(&point))
        .count() as f64;

    let b = b_count / iterations as f64;

    let pi = 4.0 - (4.0 * b);

    println!("PI might be {}", pi);
}

fn random_point(rng: &mut rand::rngs::ThreadRng) -> Point {
    let x = rng.gen::<f64>();
    let y = rng.gen::<f64>();

    Point { x, y }
}

fn is_in_circle(point: &Point) -> bool {
    (point.x * point.x) + (point.y * point.y) <= 1.0
}

And running it, we see we get:

$ cargo build -- release
$ ./target/release/monte-carlo-pi
PI might be 3.1412868

At least the first several places are right. Not bad! In the future, I might try to figure out how to make this run better. For example, if the rng is biased, or there is floating point rounding, it would be thrown off. I also might try to use some parallelism to squeeze out more iterations.

Till next time, happy learning!
– Will

Discussion

Love this post? Hate it? Is someone wrong on the internet? Is it me? Please feel free to @me on mastodon

I used to host comments on the blog, but have recently stopped. Previously posted comments will still appear, but for new ones please use the mastodon link above.

Join the conversation!