## Estimating Pi

### June 5, 2018

There’s a nice piece of math (it was either on NumberPhile or 3brown1blue, but I can’t find it again, and in any event it is a well-known equality) that says if you pick two positive integers at random, the odds of them having no common divisors are 6 ÷ π² ≈ 0.61. Let’s test that.

Your task is to estimate π using the formula shown above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Advertisements

Here is a quick and dirty way of working out this approximation, using Julia (without any dependencies whatsoever):

function HaveCommonDivisors(x::Int64, y::Int64)

m = min(x, y)

end

n = 100000 # number of runs in simulation

M = 1000000 # maximum number involved

c = 0 # counter variable

for i = 1:n

x = rand(1:M)

y = rand(1:M)

if HaveCommonDivisors(x, y); c += 1; end

end

p = 1 – c / n

PI = sqrt(6 / p)

println(p)

println(PI)

Since there is randomness involved, the results are bound to differ every time this is run. Here are some representative results:

0.60643 # probability

3.145468109123545 # π approximation

Quick perl coding – a golfish gcd calculator….

Here’s a solution in C.

Output:

Practically what Daniel did, but in Python:

3.14155312612

Here’s a solution in Racket.

Testiing.

(est-pi 100000)

n = 3.1474149880266986

p = 0.60568

There is no need to randomly sample. It is much quicker to sample a large enough range.

Check this one on the “standupmaths” YouTube channel: https://youtu.be/RZBhSi_PwHU

A Haskell version.