Yichuan Shen

Blog

  1. «
  2. »

Monte Carlo Pi

Unit circle

Today we’re gonna calculate pi. You might say, it doesn’t seem too interesting; it’s an old hat. You may be true, but we’re going to calculate pi using a method, other than approximating a circle with a polygon. So here’s what we do:

  1. We draw the unit circle on our cartesian coordinate system.
  2. We randomly drop points on the square, which contains our unit circle.
  3. The probability of the point dropping on the circle is Pi/4

Explanation

Why? The explanation is easy. The area of the unit circle is pi and the area of the square is 2^2. As the points are randomly chosen, the probability of a point dropping on the circle is

p = area of the circle/area of the square = pi/4

This gives us an interesting computational method for “calculating” pi using random numbers… Such algorithms are called Monte Carlo methods, by the way. A highly inefficient way perhaps, but nevertheless interesting. I don’t want to implement this example, in fact I wanted to blog about another Monte Carlo method for calculating pi originally, which is easy on the outside, but more difficult to grasp in its entirety.

Another Monte Carlo

  1. Randomly pick two positive natural numbers.

  2. The probability of these two numbers being coprime is

    p = 6/pi^2 = 60.79%

This gives us a formula for pi:

pi = sqrt(6/p)

Wait, you don’t believe me? With a quick demonstration, you could easily verify it’s true. Here’s a quick implementation in C#:

static Random r = new Random();

static int getGcd(int a, int b) {
    // Euclid's algorithm
    return b == 0 ? a : getGcd(b, a % b);
}

static bool doExperiment() {
    // Pick two random numbers
    int a = r.Next();
    int b = r.Next();

    // Calculate greatest common divisor
    int gcd = getGcd(a, b);

    // two numbers are coprime iff their greatest common divisor is 1
    return gcd == 1;
}

static void Main(string[] args) {
    while (true) {
        Console.Write("Count: ");

        int count = int.Parse(Console.ReadLine());
        int successes = 0;

        if (count == 0) break; // End application

        for (int i = 0; i < count; i++) {
            if (doExperiment()) successes++;
        }

        // Calculate probability
        double p = (double)successes / count;

        // Calculate pi
        Console.WriteLine("Probability: {0}%", p * 100);
        Console.WriteLine("Calculated pi: {0}", Math.Sqrt(6 / p));
        Console.WriteLine();
    }
}

You can also download the source code.

Test run

Count: 100
Probability: 59%
Calculated pi: 3.1889640207164

Count: 1000
Probability: 59.8%
Calculated pi: 3.16756133579975

Count: 10000
Probability: 60.91%
Calculated pi: 3.1385664314759

Count: 100000
Probability: 60.819%
Calculated pi: 3.14091358277565

Count: 1000000
Probability: 60.8672%
Calculated pi: 3.13966971072351

Count: 10000000
Probability: 60.78924%
Calculated pi: 3.14168232203594

Count: 100000000
Probability: 60.796216%
Calculated pi: 3.14150207229303

Count: 1000000000
Probability: 60.7912067%
Calculated pi: 3.14163150221403

As you can see, the approximation of pi is very bad and highly inefficient. The last iteration took me more than a minute to calculate, but it shows an interesting result.

Also the probability is not an approximation of 6/pi^2, it is indeed exactly 6/pi^2. Is that logical? Yes, it’s been proven1. Is that intuitive? Hell fucking no!

  1. For the mathematically challenged: Here’s a proof sketch on Wikipedia, but beware! Infinite products and the Riemann zeta function awaits you. ↩

One comment

  1. Leave your opinion
  2. Fabian says:
    January 6, 2012 at 6:21 pm

    Hey,

    this is interesting:) I just learnt how to integrate numerically using the monte carlo method – but this is fascinating!

Leave your opinion