Programming Praxis

8 Responses to “Estimating Pi”

1. Zack said

   June 5, 2018 at 12:57 PM

   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. James Curtis-Smith said

   June 5, 2018 at 4:09 PM

   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

   June 6, 2018 at 2:17 AM

   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

   June 7, 2018 at 2:29 PM

   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

   June 8, 2018 at 5:57 PM

   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

   June 11, 2018 at 7:31 AM

   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

   June 14, 2018 at 9:19 AM

   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

   July 14, 2018 at 6:01 AM

   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