Calculating $\pi$
$\pi$ is a mathematical constant defined as the ratio of the circumference of a circle to its diameter. There are many other definitions including the ratio of the area of a circle to the square of its radius. It is a constant that also appears in many formulae used in mathematics, physics and engineering.
$\pi$ is an irrational number. It cannot be expressed as a ratio of two whole numbers $a/b$. Throughout history ingenious ways have been devised to calculate this constant with increasing accuracy. In this article I'll describe three methods to determine a decimal representation of $\pi$ (3.14159...):
 The Monte Carlo method where we use a statistical approach to estimate the area of a circle
 The Leibniz formula for $\pi$ consisting of the evaluation of a simple series
 Machin's formula for $\pi$. A more advanced method using a trigonometric relationship
The decimal representation of $\pi$ to 11 decimal places is good enough to estimate the circumference of a circle that would fit inside the earth with an error of less than 1 mm (see reference (1)).
Monte Carlo method
One method of determining the value of $\pi$ is to calculate the area of a circle inscribed in a square. In the figure below, the sides of the square are of length d, equal to the diameter of the circle.
The area of a circle is defined as:
$$ A_{circle} = \pi r^2 \qquad _{...(1)} $$
$$ A_{circle} = \pi \left( \frac{d}{2} \right)^2 = \frac{\pi d^2}{4} \qquad _{...(2)}$$
The area of the square is $d^2$. Therefore $\pi$ can be calculate by the following expression
$$\pi = 4 . \frac{A_{circle}}{A_{square}} \qquad _{...(3)}$$
A statistical method can be used to determine the ratio of the area of the circle to the square. A number of random points are counted inside the square. The probability of random points found inside the circle is proportional to the area of the circle relative to the square.
The algorithm can be described as follows. Let N denote the total number of random points.
 Generate a random point inside the square with coordinates $(x_i, y_i)$. $x_i$ and $y_i$ are independent random numbers in the range r to r
 Calculate the distance of the point from the origin. $\sqrt{x_i^2+y_i^2}$

Count the number of points inside the circle. In the diagram below $(x_1, y_1)$ is inside the circle and $(x_2, y_2)$ is outside the circle
 Inside the circle: $\sqrt{x_1^2+y_1^2} \le r$
 Outside the circle: $\sqrt{x_2^2+y_2^2} \gt r$
 The estimate for $\pi$ is 4 x ( number of points inside the circle / N )
In the 2nd figure (to the right) 78 random points located inside the circle compare to a total of 100. Using the equation (3) we get:
$$\pi_{(N=100)} \approx 4 \times \frac{78}{100} = 3.12 $$
The accuracy can be improved by increasing the number of random points (e.g. 1000, and 10000 points).
For N = 1000, number of points inside the circle = 783
$$\pi_{(N=1000)} \approx 4 \times \frac{783}{1000} = 3.132 $$
The experiment can be repeated a number of times to get a mean value with an error estimate. The table below consists of 5 measurements for a sample of N random points.
There are 10 times as many random points in the trial for each subsequent row.
Sample 
Trial 1  Trial 2  Trial 3  Trial 4  Trial 5  Mean ($\mu$)  $\sigma$ 
100  3.16000  3.04000  3.28000  3.28000  3.28000  3.208000  0.107331 
1000  3.07200  3.14400  3.22400  3.16800  3.23200  3.168000  0.065238 
10,000  3.14800  3.12600  3.17720  3.12400  3.14040  3.143120  0.021514 
100,000  3.13752  3.14116  3.13988  3.13392  3.14536  3.139568  0.004251 
1,000,000  3.14412  3.14187  3.14189  3.14191  3.13887  3.141734  0.001869 
10^{7}  3.14203  3.14099  3.14102  3.14221  3.14213  3.141677  0.000613 
10^{8}  3.14172  3.14181  3.14163  3.14189  3.14159  3.141728  0.000124 
10^{9}  3.14153  3.14165  3.14169  3.14167  3.14166  3.141640  0.000063 
10^{10}  3.14161  3.14158  3.14159  3.14158  3.14161  3.141593  0.000015 
10^{11}  3.14159  3.14159  3.14161  3.14159  3.14160  3.141594  0.000009 
Leibniz formula for $\pi$
One can determine the value of $\pi$ by summing an series of terms as shown in equation (4) below
$$ \frac{\pi}{4} = \sum_{n=0}^\infty \frac{(1)^n}{2n+1} = 1  \frac{1}{3} + \frac{1}{5}  \frac{1}{7} + \frac{1}{9}  \frac{1}{11} + \frac{1}{13}  \frac{1}{15} + ... \qquad _{...(4)} $$
The formula is easy to remember. The denominator of each term is a sequence of odd numbers and the sign alternates for each successive term. A major drawback of this method is that the sum converges very slowly. The table below indicates that a large number of terms are required to obtain an high precision estimate of $\pi$. After adding 10 million terms, the estimate is accurate to 7 decimal places.
Number of terms  Sum 
100  3.13159290 
1000  3.14059265 
10,000  3.14149265 
100,000  3.14158265 
1,000,000  3.14159165 
10,000,000  3.14159255 
Machin's formula for $\pi$
Despite the simplicity of the Leibniz formula it has been demonstrated in the previous section that the convergence rate is very slow. After summing 10 million terms, an accuracy of only 7 decimal places is achieved.
In 1706 John Machin found the following:
$$\frac{\pi}{4} = 4 \times \arctan(\frac{1}{5})  \arctan(\frac{1}{239}) \qquad_{...(5)}$$
To determine $\arctan(x)$ the following series can be used for $1 <= x <= 1$:
$$\arctan(x) = x  \frac{x^3}{3} + \frac{x^5}{5}  \frac{x^7}{7} + \frac{x^9}{9} ... \qquad_{...(6)}$$
For x = 1, the Leibniz formula for $\pi/4$ is obtained. However for small values of x the series will converge rapidly. A small number of terms are required to obtain a high precision estimate.
The table below shows the evaluation of the terms in equation (5)
Term  arctan(1/5)  sum  artan(1/239)  sum 
1^{st}  0.200000000000  0.200000000000  0.004184100418  0.004184100418 
2^{nd}  0.002666666667  0.197333333333  0.000000024417  0.004184076002 
3^{rd}  0.000064000000  0.197397333333     
4^{th}  0.000001828571  0.197395504762     
5^{th}  0.000000056889  0.197395561651     
6^{th}  0.000000001862  0.197395559789     
7^{th}  0.000000000063  0.197395559852     
8^{th}  0.000000000002  0.197395559850     
Summing 8 terms for $\arctan(\frac{1}{5})$ we get 0.197395559850. Only 2 terms for $\arctan(\frac{1}{239})$ are required to obtain 0.004184076002.
Substituting the evaluated sums into equation (5) we get:
4 x ( 4 x 0.197395559850  0.004184076002) = 3.14159265359
High Precision calculation
One can determine the value of $\pi$ using fixed point arithmetic (instead of floating point). Using this method the number one is represented as a very large number (i.e. 10^{100}). An arithmetic operation such as a multiplication or division performed on a very large number will result in a high precision calculation. Take the following example:
1 / 997 = 0.001003009027 (to 12 decimal places).
Representing 1 as 1 x 10^{100}:
1 x 10^{100} / 997 = 10030090270812437311935807422266800401203610832497492477432296890672016048144433299899699097291875
The first 10 digits match the floating point calculation. However an additional 90 digits have been determined not possible using conventional floating point arithmetic. Once the operation is complete we simply move the decimal place to the left into the correct position. This simple trick is applied to the evaluation of the terms in the series.
To avoid floating point arithmetic we can rewrite equation (5) as:
$$\frac{\pi}{4} = 4 \times arccot(5)  arccot(239) \qquad_{...(7)}$$
Substituting $u=1/x$ into equation (6) we get:
$$arccot(u) = \frac{1}{u}  \frac{1}{3u^3} + \frac{1}{5u^5}  \frac{1}{7u^7} + ... \qquad_{...(8)}$$
The problem has been reduced to the following steps:
 Evaluate arccot(5) and arccot(239) using equation (8) with fixed point arithmetic

For fixed point arithmetic operations:
 Represent 1 as a large number, depending on how many decimal places are required, e.g. for 1000 decimal places, 1 is represented as 1 x 10^{1000}
 Continue to evaluate terms until the term itself is equal to zero
 To avoid rounding errors, additional guard digits are used.^{(4)}
 Substitute arccot(5) and arccot(239) into equation (7) to get a final determination of $\pi$.
Python is a suitable high level scripting language that runs on most computing platforms and can perform high precision arithmetic using long integers. Combining the basic concepts from two python scripts (references (5) and (6)) equation (7) can be evaluated to a high precision.
The number of decimal places can be easily set in the main loop of the script.
Executing the python script we get $\pi$ to 1000 decimal places:
... calculating pi to 1000 decimal places 31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
Conclusion
Three relatively simple methods to determine the value of $\pi$ have been described and demonstrated in this article.
The Monte Carlo simulation uses random points to estimate the area of a square and an inscribed circle inside the square. Accurately determining the ratio of the areas of two geometric closed objects using random points is inefficient. Counting one million random points provides an estimate that is no more accurate than 2 decimal places. Taking the average of 5 trials with 100 billion random points results in a value accurate to 6 decimal places, allowing for the uncertainty in the average.
The Leibniz formula offers an easy way to compute $\pi$ but is also highly inefficient requiring the summation of 10 million terms to achieve an estimate accurate to 7 decimal places.
Machin's formula makes use of a trigonometric relationship containing $\arctan(x)$ which converges rapidly for small values of x. It forms the basis of more advanced methods for calculating $\pi$. One can overcome the limitations of floating point arithmetic and adopt fixed point arithmetic with long integers to estimate $\pi$ to many digits.
A summary of the results reported in this article is tabulated below:
Method  Estimate 
Monte Carlo simulation (5 trials of 100 billion points)  3.141594 +/ 0.000009 
Leibniz formula (sum of 10 million terms)  3.14159255 
Machin's formula (floating point)  3.14159265359 
Machin's formula (fixed point to 1000 decimal places)  3.141592653589793238...1989 
More sophisticated methods exist for calculating $\pi$. For further reading also check out an excellent article on wikipedia which covers a detailed cronology of the computation of pi
References
 Pi: http://en.wikipedia.org/wiki/Pi
 Monte Carlo Pi: http://math.fullerton.edu/mathews/n2003/montecarlopimod.html
 Calculating Pi with the Monte Carlo method: http://learntofish.wordpress.com/2010/10/13/calculatingpiwiththemont...
 Pi  Machin: http://www.craigwood.com/nick/articles/pimachin/
 Pi with Machin's formula: http://en.literateprograms.org/Pi_with_Machin's_formula_(Python)
 Maclaurin Series: http://www.intmath.com/seriesexpansion/2maclaurinseries.php
Appendix
Proofs
Leibniz formula for $\pi$
It can be shown from elementary trigonometry that $\tan (\pi / 4) = 1$. The inverse of this relationship is:
$$\frac{\pi}{4} = \tan^{1} \left( 1 \right) = \arctan ( 1 ) \qquad_{...(1)} $$
Any continuous function can written as a series of polynomials near some value $x=a$ using the Taylor series. A special case of the Taylor series is $a=0$ which is known as the Maclaurin series^{(5)}.
We can express $\arctan(x) = f(x)$ as an infinite sum of terms.
$$ f(x) = f(0) + f'(0)x + \frac{f''(0)}{2!}x^2 + \frac{f^{(3)}(0)}{3!}x^3 + ... + \frac{f^{n}(0)}{n!}x^n + ... \qquad_{...(2)} $$
Evaluating the higher order derivatives using the quotient rule in calculus we get the following:
Derivatives  Values 
$f(x) = \arctan(x)$  $f(0) = 0$ 
$f'(x) = \frac{1}{1+x^2}$  $f'(0) = 1$ 
$f''(x) = \frac{2x}{(1+x^2)^2} $  $f''(0) = 0$ 
$f^{(3)} (x) = \frac{2(3x^21)}{(1+x^2)^3} $  $f^{(3)}(0) = 2 $ 
$f^{(4)}(x) = \frac{24x(x^21)}{(1+x^2)^4}$  $f^{(4)}(0) = 0$ 
$f^{(5)}(x) = \frac{24(5x^410x^2+1)}{(1+x^2)^5} $  $f^{(5)}(0) = 24 $ 
If we now substitute the values into equation (2) we get the following:
$$\begin{align} \arctan(x) &= 0 + x + 0 + \frac{2}{3.2.1}x^3 + 0 + \frac{+24}{5.4.3.2.1}x^5 + ... \\ &= x  \frac{x^3}{3} + \frac{x^5}{5}  ... \end{align} \qquad_{...(3)}$$
Therefore,
$$\arctan(1) = 1  \frac{1}{3} + \frac{1}{5}  \frac{1}{7} + ... $$
Scripts
Monte Carlo simulation
The following script was used for the Monte Carlo simulation. The number of random points can be set with the variable N in the main section.
Notes:
 The number of random points, N, has been set to 1,000,000 in the main section of the script. This can be changed to whatever value is required. Obviously the larger the value of N, the longer the time required to reach a final value.
 In the simulation random points are selected inside a square of length 2 units. The inscribed circle has a radius of one unit.
 The simulation uses a whole circle even though a quadrant would have been more efficient.
#!/usr/bin/python # # Monte Carlo simulation to determine pi # # John Giovannis  Dec 2011 # import math import sys from random import random #  def hypot(x,y): a = math.sqrt(x*x + y*y) return a #  def calc_pi(max): n = 0 inside = 0 # loop until the value "max" is reached while (n < max): n += 1 # pick a random point inside a box x = 2*random()  1 y = 2*random()  1 # calculate the distance of the point from the centre h = hypot(x, y) # determine if the point is inside the circle if h <= 1.0: inside += 1 # calculate pi based on the number of points inside the circle pi = 4 * ( inside / float(n) ) print "# points inside circle = %10d " % inside print "# points inside square = %10d " % n return pi #  def main(): # number of random points to count N = 1000000 np = long(N) print "Monte Carlo simulation in progress" pi_estimate = calc_pi(np) print "pi (estimate) = %10.6f " % pi_estimate sys.stdout.flush() # =================================================== # Main if __name__ == '__main__': main()
High Precision calculation of $\pi$
The high precision python script has been slightly simplified from another source (see reference (5)) with additional comments for clarity. The main points to note are:
 The number one is represented as 10 ^{(number of decimal places + 10)}. 10 guard digits are used. In the final correction, pi is divided by 10^{10}
 In the arccot function the numerator and denominator are clearly defined
 The // is the python floor division operator. This will ensure that only long division operations are maintained in the calculation.
 The script has been set to calculate pi to 1000 decimal places. This can be set to any value.
#!/usr/bin/python # # Calculate pi to many decimal places using long arithmetic # .. based on python listing in # http://en.literateprograms.org/Pi_with_Machin's_formula_(Python) # # John Giovannis  Dec 2011 # #  def arccot(u): # calculate arccot of an integer using long arthimetic # initialise the sum and first term sum = num = one // u # precalculate u^2 u2 = u*u # initialise denominator den = 3 while 1: # evaluate the next numerator in the series num =  num // u2 # evaluate the term term = num // den # break out of the loop if the term is equal to zero if term == 0: break # update the denominator den += 2 # update the sum in the series sum += term return sum #  def calc_pi(ndp): global one # represent the number 1 (as a large number) one = 10**(ndp+10) print "... calculating pi to %d decimal places" % ndp # evaluate pi pi = 4 * (4 * arccot(5)  arccot(239)) # ... final division to allow for the guard digits pi = pi // (10**10) print "%s \n" % pi #  def main(): # calculate pi to 1000 decimal places calc_pi(1000000) # =================================================== # Main if __name__ == '__main__': main()
 Log in to post comments