# Calculating pi

I’ve been reading a wonderful little book called ‘A history of pi‘ by Petr Beckman. A few highlights:

• Archimedes’ approach to estimating pi, which was based on approximating a circle with multi-sided polygons, was the best method for 19 centuries until techniques based on infinite series were developed in the 17th century
• Dabbling in mathematics in the Middle Ages could get you into big trouble.  The Spanish mathematician Valmes was burned at the stake because he claimed to have found the solution of the quartic equation (although this story is disputed). The solution to the quartic equation is now on Wikipedia.
• As in many areas of mathematics, Newton and Euler made truly exceptional contributions.

I wondered about estimating pi using R.

R knows about pi and can correctly produce the first 16 digits which is the limit of the IEEE 754 double-precision floating-point numbers that R uses.

pi
# [1] 3.141593

op <- options(digits = 22)
pi
# [1] 3.141592653589793115998

# The correct value is
# 3.141592653589793238462

op # restore default number of digits

The  Rmpfr library provides the functionality for many more digits.

# Using Rmpfr
# http://cran.r-project.org/web/packages/Rmpfr/vignettes/Rmpfr-pkg.pdf
library(Rmpfr)
Const("pi",500) # 500 bits
[1] 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081283


The first infinite series for pi was developed by James Gregory.  He showed that $\arctan(x)$ could be expressed as follows.
$\arctan(x) = x-\frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \ldots$

If we set $x = 1$, $\arctan(x) = \frac{\pi}{4}$, that is:
$\pi = 4 \left(1 - 1/3 + 1/5 - 1/7 + \ldots \right)$

The problem is that this series converges very slowly.

In 1706, John Machin used the difference of two arctan series to develop a practical approach to estimating many digits for $\pi$.

$\frac{pi}{4} = 4 \left( \frac{1}{5} - \frac{1}{3 \cdot 5^3} + \frac{1}{5 \cdot 5^5} - \ldots \right) - \left( \frac{1}{239} - \frac{1}{3 \cdot 239^3} + \frac{1}{5 \cdot 239^5} - \ldots \right)$

We can implement this in R as follows:

i <- 1:100 # first 100 terms in the series
# pre calculate the series for arctan(1/5) and arctan(1/239)
x.5 <- mpfr(5,300)^ -(2*i-1) # using 300 bits precision
x.239 <- mpfr(239,300)^ -(2*i-1)

my.pi <- 4*(4*(sum((-1)^(i-1)*x.5/(2*i-1))) - sum((-1)^(i-1)*x.239/(2*i-1) ))
my.pi
# 1 'mpfr' number of precision  300   bits
# [1] 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348327
# which is correct up to the last 3 decimal places

Faster approaches have been developed since 1706.

Srinivasa Ramanujan’s formula for $\pi$ converges rapidly.

$\frac{1}{\pi} = \frac{2 \sqrt{2}}{9801} \sum_{n=0}^{\infty} \frac{(4n)!(1103 + 26390n)}{(n!)^4396^{4n}}$

After 1 term the approximation is correct to six decimal places, after 2 terms it’s correct to 15.

op <- options(digits = 22)

piRamanujan <- function(n){
n <- 0:n
1/(2*sqrt(2)/(9801) * sum(factorial(4*n)  * (1103 + 26390 * n) / (factorial(n)^4 *396^(4*n))))
}

piRamanujan(0)
[1] 3.141592 730013305523329

piRamanujan(1)
[1] 3.141592653589793 560087
# We are already at the limit of double-precision floating-point arithmatic

A faster converging formula was developed by Gregory and David Chudnovsky (as discussed in Alex’s Adventures in Numberland) and Wikipedia.

The formula listed in Alex’s Adventures in Numberland is:

$\frac{1}{\pi} = \sum_{n=0}^{\infty} (-1)^n \times \frac{(6n)!}{(3n)!n!^3} \times \frac{163096908 + 6541681608n}{(262537412640768000)^{n + \frac{1}{2}}}$

piChudnovsky <- function(n){

There is a good article of the history of calculating pi in Wikipedia.   $\pi$has been calculated to 1.21 trillion digits (as of 28 Dec 2013).