PI calculation

with random numbers

Random PI

Introduction

The idea was to approximate the value of PI by filling a squared area with randomly generated points and checking which points fell into the circle contained in the square, given that the circle and the square have a ratio of areas that is π/4,

Check the repo.

The maths

All you need to define a point in a 2-dimensional plane is a pair of coordinates:

In rectangular (or Cartesian) coordinates, you can uniquely define a point by the pair (X, Y).
In cylindrical coordinates, a point is defined instead by a different pair of coordinates, (ρ,φ), where ρ represents the distance from the point to the origin of coordinates, and φ represents the angle of the point with the X axis^[*].

You can easily transform from one system to the other if you write X and Y in terms of ρ and φ:

$X = ρ \cos (φ)$
$Y = ρ \sin (φ)$

We can also do the opposite and write the cylindrical coordinates in terms of the rectangular coordinates:

$ρ = \sqrt{X^{2} + Y^{2}}$
$\tan (φ) = \frac{Y}{X}$

An image is worth a million words, so check this diagram:

Diagram of the cylindrical coordinates

Image by Andeggs, Wikimedia Commons

In a circle, we can use the cylindrical coordinates to calculate things more easily. A circle has a fixed value of ρ that we call radius (r), and the area of a circle is defined by:

$A = π r^{2}$

The code

The total number of points that we launch to fill the squared area are represented by ntot, while nin are the points that fall into the circle.

A random number is generated for the X and the Y coordinates. The maximum value returned by rand() is defined as the constant RAND_MAX. This value is library-dependent, but is guaranteed to be at least 32767 on any standard library implementation.

If the radius of the circle is 1, the area will be π. So, first, we generate a random number between 0 and 1, for which we call the rand() function and divide by RAND_MAX + 1:

double x = (double)rand() / (RAND_MAX + 1.0);
double y = (double)rand() / (RAND_MAX + 1.0);

This will generate points that fall into the first quadrant only. Then, we transform these random numbers so that they are laid over a circle of radius 1 centered in the origin of coordinates:

xc = (2*x) - 1;
yc = (2*y) - 1;

Now, to check if this point (X, Y) is inside the circle of radius r = 1, we calculate the squared radius of the circle that has X and Y values between 0 and 1:

$r = \sqrt{X^{2} + Y^{2}}$

int nin = 0;
double r = sqrt( (xc*xc) + (yc*yc) );
	
if (r <= 1)
	nin++;

Since the ratio of areas is proportional to π/4 (i.e., nin/ntot = π/4), we can calculate π as:

double PI = 4*( (double)nin / (double)ntot );

Macros

The first macro calculates random numbers up to 1e8. The second calculates from 1000 to 1e10 in steps of 10. The third calculates from 1e4 to 2e5 in steps of 1e4 and writes the results to a file. The results can be seen in this plot:

Plot of PI

Check the repo.

Share on Twitter or follow me!

[*] This is so in order to avoid confusion with spherical coordinates, where there is also a coordinate that represents the distance from the origin (r), and the angle with the Z axis (θ). The cylindrical coordinates are then (r, θ, φ), where ρ = r senθ, so,
X = r sen(θ) cos(φ),
Y = r sen(θ) sen(φ),
Z = r cos(θ)
Spherical and cylindrical coordinates

Published: 2004-05-27 19:56