Calculating Pi

October 9, 2009

The function that solves the first task uses a do loop accumulates the counts. The sand? function simulates casting a single grain of sand, and is #t if it lands inside the circle:

(define (pi n)
    (define (square n) (* n n))
    (define (sand?) (< (+ (square (rand)) (square (rand))) 1))
    (do ((i 0 (+ i 1)) (p 0 (+ p (if (sand?) 1 0))))
        ((= i n) (exact->inexact (* 4 p (/ n))))))

We use rand from the Standard Prelude, and convert from rational to real with the exact->inexact function. The process converges slowly; after casting a hundred thousand grains of sand, π is 3.14188, in which only the first four digits are correct.

> (pi 100000)
3.14188

The function that solves the second task is:

(define (bound-pi n)
  (let loop ((a (* 3 (sqrt 3)))
             (b (* 3 (sqrt 3) (/ 2)))
             (n n))
    (if (= n 1)
        (values b a)
        (let ((next-a (/ 2 (+ (/ a) (/ b)))))
          (loop next-a (sqrt (* next-a b)) (- n 1))))))

Archimedes’ calculation is:

> (bound-pi 6)
3.14103195089051
3.1427145996453683

The calculation converges quickly; (bound-pi 27) = 3.1415926535897936 is the limit of double-precision arithmetic, using inscribed and circumscribed polygons of 201326592 sides to calculate the bounds.

You can run the code at http://programmingpraxis.codepad.org/3LMdEHBZ.

About these ads

Pages: 1 2

11 Responses to “Calculating Pi”

  1. [...] Praxis – Calculating Pi By Remco Niemeijer In today’s Programming Praxis problem we have to implement two ways of calculating pi. Let’s get [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2009/10/09/programming-praxis-calculating-pi/ for a version with comments):

    import Control.Monad
    import System.Random
    
    monteCarloPi :: Int -> IO ()
    monteCarloPi n = do
        hits <- fmap sum $ liftM2 (zipWith checkHit) rs rs
        print $ fromIntegral hits / fromIntegral n
        where rs = replicateM n $ randomRIO (0,1 :: Double)
              checkHit x y = if x*x + y*y < 1 then 4 else 0
    
    boundPi :: Int -> (Double, Double)
    boundPi n = iterate f (3/2 * sqrt 3, 3 * sqrt 3) !! (n - 1)
                where f (b, a) = let x = 2 / (1 / a + 1 / b)
                                 in (sqrt $ b * x, x)
    
  3. kbob said

    Instead of the Monte Carlo method, I started by generating a uniform grid of n x n points on the first quadrant (0 <= x, y <= 1). pi/4 of them fell inside the unit circle. That ran rather slowly and wasn't converging very well.

    Then I realized I could use the Bresenham circle algorithm and just calculate one point in each row. Since the basic BCA only finds points through a 45 degree angle, I ran it for 1/8th of the circle to find pi/8.

    Here's how quickly it converged, using up to a billion iterations.

    1: -8
    10: 2.8
    100: 3.1
    1000: 3.137576
    10000: 3.1411928
    100000: 3.1415526136
    1000000: 3.141588652488
    10000000: 3.14159225361224
    100000000: 3.14159261359012
    1000000000: 3.14159264958976
    

    The code is deceptively simple. In C, based on the sample code in Wikipedia.

    #include 
    
    double pi(int n)
    {
        int f = 1 - n;
        int ddF_x = 1;
        int ddF_y = -2 * n;
        int x = 0;
        int y = n;
        long long int in;
        while (x < y) {
    	if (f >= 0) {
    	    --y;
    	    ddF_y += 2;
    	    f += ddF_y;
    	}
    	x++;
    	ddF_x += 2;
    	f += ddF_x;
    	in += y - x;
        }
        return 8.0 * in / ((double)n * n);
    }
    
    main()
    {
        int i;
    
        for (i = 1; i <= 1000000000; i *= 10)
    	printf("%d: %.15g\n", i, pi(i));
    }
    
    

    (I wonder if the formatting will survive…)

  4. Maurits said

    >perl -e “$m=pop;for(1..$m){$n+=rand()**2+rand()**2<1}print$n*4/$m" 5000000
    3.1419256

  5. Rudi Angela said

    Here’s an implementation in Erlang.
    Approach:
    Function calc() takes parameter Tally that is the number of times to produce a random coordinate (X and Y between 0 and 1).
    Calculate the distance to the origin and if it is less than 1 then count it as a match (inside circle with radius 1).
    Keep a count on the total number of tries (Done) executed.
    At the end (when Tally = 0) just calculate Pi as 4 * Matches / Done.

    calc(Tally) ->
    calc(Tally, 0, 0).

    calc(0, Matches, Done) -> 4 * Matches / Done;

    calc(Tally, Matches, Done) ->
    X = random:uniform(), Y = random:uniform(),
    calc (Tally – 1, if X*X + Y*Y Matches + 1; true -> Matches end, Done + 1).

  6. Just for fun, a visual approach to the Monte Carlo simulation using Processing (http://processing.org). Unfortunately this program will never truly converge towards Pi as the whitespace is only a pixel-based approximation of a quarter-circle.

    You can view the results as an applet here: http://stungeye.com/processing/1009/

    final int R = 400;
    final color BLUE  = color(0, 0, 255);
    final color WHITE = color(255, 255, 255);
    final color RED   = color(255, 0, 0);
    final color BLACK = color(0, 0, 0);
    
    double n = 0;
    double count = 0;
    
    void setup () {
      size(R,R);
      background(0);
      fill(255);
      stroke(255);
      ellipse(0, R, R*2, R*2);
    }
    
    void draw() {
      count++;  
      
      int random_x = (int)random(R);
      int random_y = (int)random(R);
      
      loadPixels();
      color current_pixel_colour = pixels[random_y * R + random_x];
      
      if ((current_pixel_colour == WHITE) || (current_pixel_colour == BLUE)) {
        n++;
        set(random_x, random_y, BLUE);
        println(4 * n / count);
      } else {
        set(random_x, random_y, RED);
      }  
    }
    

    (The applet version of the code is slightly more complex: http://stungeye.com/processing/1009/monte_carlo_pi.pde)

  7. Monte Carlo method in python:

    from random import random
    def pi(n):
    return 4*float(sum(1 if (random()**2 + random()**2) <= 1 else 0 for i in range(n)))/n

  8. [...] calculated the value of pi, and logarithms to seven digits, in two previous exercises. We continue that thread in the current exercise with a function that calculates [...]

  9. Brian said

    Hi-

    What is the advantage of the monte carlo method – randomly picking n points where the x,y are independent and uniformly distributed – and the grid method, where you set up points on some defined lattice and then just keep making the lattice finer and finer?

    Thanks!

  10. [...] (diesmal von Programming Praxis). Eine sehr beliebte Aufgabe ist die Berechnung von Pi, eine Übung die ich auch meinen Schülern gerne gebe – oft mit der Zusatzaufgabe zu [...]

  11. David said

    I thought it would be fun to do first exercise using integer arithmetic only (fixed point in other words) random-ui is a FORTH word which returns a random 32 bit integer, which can be interpreted as a 32 bit binary “floating point” number between 0 and 1 (“ui” suffix stands for “unit interval”…)

    \ multiply with carry RNG
    \ quickie MWC-0  (no lag, 2^63 period is good enough)
    
    create seed  counter ,
    create carry 48313 ,
    
    : random-ui  ( -- r )
       seed @ 4164903690 um*  carry @ 0 d+
       carry ! dup seed ! ;
    
    : random  ( n -- r )
       random-ui swap um* nip ;
    
    
    \ Calculating pi in FORTH!
    \ Monte-Carlo method, fixed point
    
    : sq-ui  ( n -- n )   \ square unit interval fixed point
        dup um* nip ;
    
    : sand?  ( -- ? )       \ TRUE - in circle, FALSE - in square
        random-ui sq-ui 0   \ square and convert to double
        random-ui sq-ui 0
        d+ nip 0= ;
    
    : monte-pi  ( n -- )
        locals| n |
        0
        n 0 DO  sand? 1 and  +  LOOP
        4000000 n */
        0 <# # # # # # # [char] . hold #s #> type space ;
    

    A few runs, last one with max 32 bit signed integer value

    include pi  ok
    100000000 monte-pi 3.141614  ok
    100000000 monte-pi 3.141501  ok
    $7FFFFFFF monte-pi 3.141553  ok
    

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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 612 other followers

%d bloggers like this: