## 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):

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.

(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/

(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”…)

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