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.
Programming Praxis 
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….
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}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[0]); return EXIT_FAILURE; } int n_trials = DEFAULT_N_TRIALS; if (argc == 2) n_trials = atoi(argv[1]); 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:
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
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
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.141592653589793Check this one on the “standupmaths” YouTube channel: https://youtu.be/RZBhSi_PwHU
A Haskell version.
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)