Calculating Pi
October 9, 2009
Pi, or π, is a mathematical constant with a value that is the ratio of a circle’s circumference to its diameter. It has been known since antiquity that the ratio of the circumference of a circle to its diameter is the same for all circles; the ancient Egyptians, Indians, and Babylonians had all calculated approximations for π about two millenia before the birth of Christ. π, which is approximately equal to 3.14159, is one of the most important constants in math, science and engineering; it pops up regularly in studies of geometry, trigonometry and calculus, Einstein used π in his field equation of general relativity, and it appears in many statistical distributions. In a previous exercise we used a spigot algorithm to calculate the digits of π; our exercise today will use two different methods to calculate the value of π.
An interesting method for calculating π uses Monte Carlo simulation. If a circle of radius r is inscribed in a square with sides of length 2r, the area of the circle will be πr2 and the area of the square will be (2r)2, so the ratio of the area of the circle to the area of the square will be π/4. Another way of looking at this, as on the diagram, is to consider just the first quadrant of the circle; the square has an area of r 2, and the portion of the circle within the square has an area of πr2/4.
By taking a large number of points randomly distributed throughout the square and counting how many are within the inscribed circle, we can estimate the value of π. We could do that by building a model, scattering sand over it, and counting the individual grains of sand, but since we are programmers, it is easier to write a program to do the counting for us.
Your first task is to implement a program to calculate the value of π using the Monte Carlo method described above.
The second method is due to Archimedes (287–212 BC), a Greek mathematician who lived in Syracuse, who famously bounded the value of π within a small range by measuring the perimeters of inscribed and circumscribed regular polygons with ninety-six sides: 223/71 < π < 22/7.
Consider a circle with radius 1 and circumference 2π in which regular polygons of 3 × 2n-1 sides are inscribed and circumscribed; the diagram for n = 2 is shown at right. If bn is the semiperimeter of the inscribed polygon, and an is the semiperimeter of the circumscribed polygon, then as n increases, b1, b2, b3, … defines an increasing sequence, and a1, a2, a3, … defines a decreasing sequence, each with limit π.
Given K = 3 × 2n-1, the semiperimeters are an = K tan(π/K) and bn = K sin(π/K) by the definitions of sine and tangent. Likewise, an+1 = 2K tan(π/2K) and a n+1 = 2K sin(π/2K).
Then, simple trigonometry allows us to calculate (1/an + 1/bn) = 2/an+1 and an+1bn = (bn+1)2. Archimedes started with a1 = 3 tan(π/3) = 3√3 and b1 = 3 sin(π/3) = 3√3/2 and calculated b6 < π < a6.
Archimedes, of course, didn’t have trigonometry available to him, as it hadn’t been invented yet; he had to work out the geometry directly, as well as making all the calculations by hand!
Your second task is to write a function that calculates the bounds of π using Archimedes’ algorithm. You can test your function for n = 6, as Archimedes did. When you are finished, you are welcome to read or run a suggested solution, or to post your solution or discuss the exercise in the comments below.
Programming Praxis 
[…] 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
[…] https://programmingpraxis.com/2009/10/09/calculating-pi/ […]