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.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
12 Comments »

12 Responses to “Calculating Pi”

  1. Programming Praxis – Calculating Pi « Bonsai Code said
    October 9, 2009 at 6:08 PM

    […] 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
    October 9, 2009 at 6:09 PM

    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
    October 10, 2009 at 3:58 AM

    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
    October 12, 2009 at 4:52 PM

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

  5. Rudi Angela said
    October 20, 2009 at 10:01 PM

    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. Wally Glutton said
    October 30, 2009 at 5:10 AM

    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. Ronald Kaiser said
    November 12, 2009 at 2:37 AM

    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. Calculating Sines « Programming Praxis said
    January 12, 2010 at 9:01 AM

    […] 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
    February 8, 2010 at 9:24 PM

    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. Berechnung von Pi | ComSubVies Creative Content said
    March 17, 2010 at 1:18 PM

    […] (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
    May 16, 2012 at 1:22 PM

    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 comment