## 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.

Pages: 1 2

### 8 Responses to “Estimating Pi”

1. Zack said

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)

``````if iseven(x) && iseven(y); return true; end

for d = 3:2:m
if (x / d == round(x / d)) && (y / d == round(y / d)); return true; end
end

return false
``````

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

2. Quick perl coding – a golfish gcd calculator….

```use strict;
use warnings;
use feature 'say';
use Const::Fast qw(const);

const my\$N=>100000;
const my\$C=>100000;

say sqrt\$C*6/(\$C-grep{g(r(),r())-1}1..\$C);

sub r{1+int\$N*rand}
sub g{my(\$a,\$b)=@_;\$b>\$a?g(\$b,\$a):\$a%\$b?g(\$b,\$a%\$b):\$b}
```
3. Daniel said

Here’s a solution in C.

```#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define DEFAULT_N_TRIALS 1000000

unsigned gcd(unsigned a, unsigned b) {
if (a > b) {
unsigned t = a;
a = b;
b = t;
}
while (true) {
if (a == 0) return b;
unsigned t = a;
a = b % a;
b = t;
}
}

int main(int argc, char* argv[]) {
if (argc > 2) {
fprintf(stderr, "Usage: %s [N_TRIALS]\n", argv);
return EXIT_FAILURE;
}
int n_trials = DEFAULT_N_TRIALS;
if (argc == 2) n_trials = atoi(argv);
struct timespec t;
clock_gettime(CLOCK_MONOTONIC, &t);
srand(t.tv_nsec);
unsigned count = 0;
for (int i = 0; i < n_trials; ++i) {
int a = 0;
int b = 0;
while (!(a = rand()));
while (!(b = rand()));
count += gcd(a, b) == 1;
}
double pi = sqrt(6.0 * n_trials / count);
printf("%f\n", pi);
return EXIT_SUCCESS;
}
```

Output:

```3.142567
```
4. Milbrae said

Practically what Daniel did, but in Python:

```import random
import math

def gcd(a, b):
''
while b:
a, b = b, a % b

return abs(a)

if __name__ == "__main__":
''
N = 10000000
C = 0

for k in range(N):
a = random.randrange(1, N * 100)
b = random.randrange(1, N * 100)
if gcd(a, b) == 1:
C += 1

print (math.sqrt(6.0 * N / C))
```

3.14155312612

5. Kevin said

Here’s a solution in Racket.

```(define (est-pi sz)
(let loop ((i sz) (tot 0))
(if (zero? i)
(let* ((avg (/ tot sz)) (prob (- 1 avg)) (pie (sqrt (/ 6 prob))))
(display (format "n = ~a\np = ~a\n" pie (exact->inexact prob))))
(let* ((sd 4294967087) (hit? (> (gcd (random sd) (random sd)) 1)))
(loop (sub1 i) (+ (if hit? 1 0) tot))))))
```

Testiing.

(est-pi 100000)
n = 3.1474149880266986
p = 0.60568

6. Paul said

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

```from math import gcd, pi, sqrt

A, M = 2, 10000
score = sum(1 for i in range(A, A+M) for j in range(A, A+M) if gcd(i,j) == 1)
calc_pi = sqrt((6*M**2) / score)
print("calculated: {}  vs. pi: {}".format(calc_pi, pi))
# -> calculated: 3.1415450131521814  vs. pi: 3.141592653589793
```
7. axiac said

There’s a nice piece of math (it was either on NumberPhile or 3brown1blue, but I can’t find it again)

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

8. Globules said

```import qualified Control.Foldl as F
import Data.Bool (bool)
import Data.List.Ordered (isectBy)
import Data.Word (Word64)
import System.Random (RandomGen, randoms)
import System.Random.TF.Init (initTFGen)

type W = Word64

-- A sequence of (hopefully) converging estimates to π.
piEstimates :: RandomGen r => r -> [Double]
piEstimates r = let rs = randoms r :: [W]
ps = zip rs (tail rs)
in F.postscan (approx <\$> sumLen) \$ map coprime ps

-- 1 if m and n are coprime, else 0.
coprime :: Integral a => (a, a) -> Int
coprime (m, n) = bool 0 1 (gcd m n == 1)

-- An approximation to π, given the fraction of random numbers that are coprime.
approx :: Integral a => (a, a) -> Double
approx (x, y) = sqrt (6.0 / (fromIntegral x / fromIntegral y))

main :: IO ()
main = do
rand <- initTFGen
mapM_ print \$ pick (piEstimates rand) [10, 100, 1000, 10000, 100000, 1000000]

-- --------------------------------- Utilities ---------------------------------

-- Fold over integers, producing their sum and length.
sumLen :: F.Fold Int (Int, Int)
sumLen = (,) <\$> F.sum <*> F.length

-- xs `pick` ns is the sublist of xs having the indices ns.  It's assumed that
-- ns is strictly increasing.
pick :: [a] -> [Int] -> [a]
pick xs = map snd . isectBy (\(i, _) j -> i `compare` j) (zip [0..] xs)
```
```\$ ./piest
2.70801280154532
2.805375174576454
3.122498999199199
3.154036191466773
3.1458988695068744
3.1415919037069258
```