Thursday, November 12, 2015

How to calculate Pi

$\pi$ is an irrational number: it cannot be expressed as the quotent of two integers. Furthermore the decimal representation of an irrational number does not terminate and does not consist of a repeating sequence. How can such a number ever be fully communicated? The informational content of an irrational number is infinite, therefore the number must be represented indirectly as with an expression.

A convenient way to represent $\pi$ is with an infinite series. When the terms of an infinite series are added together they can converge to a single value or they can diverge by increasing to infinity or by oscillating within a range of values. For the purposes of representing $\pi$ we need to use a convergent infinite sequence.

The Nilakantha Series is easy to understand and easy to calculate but converges to $\pi$ slowly. \[ \pi = 3 +\frac{4}{2\times3\times4} -\frac{4}{4\times5\times6} +\frac{4}{6\times7\times8} -\frac{4}{8\times9\times10}\cdots \] The Nilakantha Series can be expressed in sigma notation as: \[\pi = 3+\displaystyle\sum_{k=1}^{\infty} \frac{4(-1)^{k-1}}{2k\times (2k+1) \times (2k+2)} \] \[\pi = 3+4\displaystyle\sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{8k^{3}+12k^{2}+4k} \] The Chudnovsky Series is more complex but converges to $\pi$ quickly. \[\frac{1}{\pi} = 12\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}640320^{\frac{3}{2}}}\] This blog will demonstrate how to calculate the Chudnovsky series in C++ and will attempt to do so with sufficient detail such that any motivated novice can understand the process. The Chudnovsky series can be reformulated to allow for a more convenient calculation.
\[\frac{1}{\pi} = 12\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}640320^{\frac{3}{2}}}\] \[\frac{1}{\pi} = \frac{12}{640320^{\frac{3}{2}}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}}\] \[\frac{1}{\pi} = \frac{12}{640320^{1}\times640320^{\frac{1}{2}}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}}\] \[\frac{1}{\pi} = \frac{12}{640320\times\sqrt{640320}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}}\] \[\frac{1}{\pi} = \frac{12}{12\times53360\times\sqrt{64\times10005}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}}\] \[\frac{1}{\pi} = \frac{1}{53360\times8\times\sqrt{10005}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}}\] \[\frac{1}{\pi} = \frac{1}{426880\sqrt{10005}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!(13591409+545140134k)}{(3k)!(k!)^{3}640320^{3k}}\] \[\frac{1}{\pi} = \frac{1}{426880\sqrt{10005}}\displaystyle\sum_{k=1}^{\infty}\frac{(-1)^{k}(6k)!13591409}{(3k)!(k!)^{3}640320^{3k}}+\frac{(-1)^{k}(6k)!545140134k}{(3k)!(k!)^{3}640320^{3k}}\] \[a = \frac{(-1)^{k}(6k)!}{(3k)!(k!)^{3}640320^{3k}}\] \[b = \frac{(-1)^{k}(6k)!k}{(3k)!(k!)^{3}640320^{3k}}\] \[b = a\times k\] \[\frac{1}{\pi} = \frac{1}{426880\sqrt{10005}}\displaystyle\sum_{k=1}^{\infty}13591409a+545140134b\] \[\frac{1}{\pi} = \displaystyle\sum_{k=1}^{\infty}\frac{13591409a+545140134b}{426880\sqrt{10005}}\] \[\pi = \displaystyle\sum_{k=1}^{\infty}\frac{426880\sqrt{10005}}{13591409a+545140134b}\] At this point the series has a favorable form where knowing a allows b to be quickly calculated. This is further improved by finding a multiple such that ak can be calculated quickly from a(k-1). \[\frac{a_k}{a_{(k-1)}}=\frac{(-1)^{k}(6k)!}{(3k)!(k!)^{3}640320^{3k}}\div\frac{(-1)^{k-1}(6(k-1))!}{(3(k-1))!((k-1)!)^{3}640320^{3(k-1)}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{(-1)^{k}(6k)!}{(3k)!(k!)^{3}640320^{3k}}\times\frac{(3(k-1))!((k-1)!)^{3}640320^{3(k-1)}}{(-1)^{k-1}(6(k-1))!}\] \[\frac{a_k}{a_{(k-1)}}=\frac{(-1)^{k}(6k)!(3(k-1))!((k-1)!)^{3}640320^{3(k-1)}}{(-1)^{k-1}(6(k-1))!(3k)!(k!)^{3}640320^{3k}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{(-1)^{k}(6k)!(3k-3)!((k-1)!)^{3}640320^{3k-3}}{(-1)^{k}(-1)^{-1}(6k-6)!(3k)!(k!)^{3}640320^{3k}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{(6k)!(3k-3)!((k-1)!)^{3}640320^{3k}640320^{-3}}{(-1)^{-1}(6k-6)!(3k)!(k!)^{3}640320^{3k}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{-1(6k)!(3k-3)!((k-1)!)^{3}640320^{-3}}{(6k-6)!(3k)!(k\times(k-1)!)^{3}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{-1(6k)(6k-1)(6k-2)(6k-3)(6k-4)(6k-5)(6k-6)!(3k-3)!((k-1)!)^{3}}{(6k-6)!(3k)(3k-1)(3k-2)(3k-3)!k^{3}((k-1)!)^{3}640320^{3}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{-6k(6k-1)(6k-2)(6k-3)(6k-4)(6k-5)}{3k(3k-1)(3k-2)k^{3}640320^{3}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{-2(6k-1)(2)(3k-1)(3)(2k-1)(2)(3k-2)(6k-5)}{(3k-1)(3k-2)k^{3}640320^{3}}\] \[\frac{a_k}{a_{(k-1)}}=\frac{-24(6k-1)(2k-1)(6k-5)}{k^{3}640320^{3}}\] Here is a basic implementation of the algorithm in C++. The code is free to use on an as is basis. No attribution is required if you use it. No warranty of any kind is implied.

void chudnovsky (void)
{

long double const const_A = 13591409.0;
long double const const_B = 545140134.0;
long double const const_C = 640320.0;
long double const const_D = 426880.0;
long double const const_E = 10005.0;
long double squareRoot = 0;
long double c_Cubed = pow(const_C, 3); // 262537412640768000
long double numerator = 0;
long double denominator = 0;
long double piAprx = 0;
long double alpha = 1; // at k = 0
long double alpha_Sigma = 1;
long double beta_Sigma = 0; // at k = 0
int iterations = 1;

squareRoot = sqrt(const_E);
numerator = squareRoot*const_D;
denominator = (alpha_Sigma*const_A) + (beta_Sigma*const_B);
piAprx = numerator/denominator;
std::cout << "at k = 0, piAprx = " << std::setprecision(16) << piAprx << std::endl << std::endl;

for (long double k = 1; k <= iterations; k++)
{
// this starts calculating the series from k = 1

alpha *= -24*(6*k-1)*(2*k-1)*(6*k-5)/(pow(k,3)*c_Cubed);

alpha_Sigma += alpha;

beta_Sigma += alpha*k;

}
denominator = (alpha_Sigma*const_A) + (beta_Sigma*const_B);

piAprx = numerator/denominator;
std::cout << "at k = " << iterations << " piAprx = " << std::setprecision(16) << piAprx << std::endl << std::endl;

}

3 comments:

  1. This comment has been removed by the author.

    ReplyDelete
  2. Just came here to say... thank you immensely for this. I'm grateful that you actually took the time to help me understand this.

    ReplyDelete
  3. Just came here to say... thank you immensely for this. I'm grateful that you actually taking the time to help me understand this.

    ReplyDelete