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

### 12 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

```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);

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
```