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

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[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:

    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

    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)
    
    $ ./piest
    2.70801280154532
    2.805375174576454
    3.122498999199199
    3.154036191466773
    3.1458988695068744
    3.1415919037069258
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: