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.
[…] 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
[…] https://programmingpraxis.com/2009/10/09/calculating-pi/ […]