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
[…] Praxis – Calculating Pi By Remco Niemeijer In today’s Programming Praxis problem we have to implement two ways of calculating pi. Let’s get […]
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)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.
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…)
>perl -e “$m=pop;for(1..$m){$n+=rand()**2+rand()**2<1}print$n*4/$m" 5000000
3.1419256
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).
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)
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
[…] 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 […]
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!
[…] (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 […]
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