## Estimating Pi

### June 5, 2018

It’s easy enough to write the program:

(define (rand50) (inexact->exact (random #e1e50)))

(define (pi n) (let loop ((m n) (k 0)) (if (zero? m) (sqrt (* 6 n (/ k))) (let ((a (rand50)) (b (rand50))) (loop (- m 1) (if (= 1 (gcd a b)) (+ k 1) k))))))

We use the native Chez Scheme random number generator. The result is not particularly quick to converge, and I’m really surprised at the result for *n* = 100:

> (do ((e 1 (+ e 1))) ((= e 9)) (display e) (display #\tab) (display (pi (expt 10 e))) (newline)) 1 3.1622776601683795 2 2.9925280083228984 3 3.1517891481565017 4 3.1523110971785706 5 3.142283033927125 6 3.1420141701765227 7 3.1416836140755615 8 3.141576664431344

You can run the program at https://ideone.com/0jqRSq.

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.