Optimizing the Chudnovsky Algorithm for Implementations

~ bernard2518141184

The Chudnovsky Algorithm is a fast method for calculating the digits of \pi. Here is the definition:

\displaystyle\frac{1}{\pi}=12\sum_{k=0}^{\infty}\frac{(-1)^k(6k)!(545140134k+13591409)}{(3k)!(k!)^3640320^{3k+\frac{3}{2}}}

The complexity of the algorithm is very high. It involves multiple factorials, exponents and calculations with large numbers. You should not directly implement the it as its definition, since the program will be pretty slow. Here shows some optimization of the algorithm for computations.

Mathematical Optimizations

\displaystyle\frac{1}{640320^{\frac{3}{2}}} can be factorized from the infinite sum as a constant to be computed first before the summation. Such factorization reduces the computation of that constant to only one time. The equation then becomes:

\displaystyle\frac{640320^{\frac{3}{2}}}{\pi}=12\sum_{k=0}^{\infty}\frac{(-1)^k(6k)!(545140134k+13591409)}{(3k)!(k!)^3640320^{3k}}

It can be then changed into surd form.

\displaystyle640320^{\frac{3}{2}}\\=\sqrt{640320}^3\\=(8\sqrt{10005})^3\\=5122560\sqrt{10005}

We can further divide 12 by both sides. The equation then becomes:

\displaystyle\frac{426880\sqrt{10005}}{\pi}=\sum_{k=0}^{\infty}\frac{(-1)^k(6k)!(545140134k+13591409)}{(3k)!(k!)^3640320^{3k}}

According to the exponential rule \displaystyle a^{b\cdot c}\equiv(a^b)^c, we can compute \displaystyle640320^3 in the denominator first, which is equal to \displaystyle262537412640768000.

The \displaystyle(-1)^k can also be factorized with the part above. Since \displaystyle(-1)^k=(-1)^{-k}, we can move it to the denominator. Then, factorize the exponent \displaystyle k from \displaystyle(-1)^k and \displaystyle262537412640768000^k, which results in \displaystyle(-262537412640768000)^k. Now, the equation becomes:

\displaystyle\frac{426880\sqrt{10005}}{\pi}=\sum_{k=0}^{\infty}\frac{(6k)!(545140134k+13591409)}{(3k)!(k!)^3(-262537412640768000)^k}

\displaystyle\pi=\frac{426880\sqrt{10005}}{\displaystyle\sum_{k=0}^{\infty}\frac{(6k)!(545140134k+13591409)}{(3k)!(k!)^3(-262537412640768000)^k}}

Programmatic Optimizations

You won’t just hard implement the entire algorithm into a program. Every time a loop executes, it recalculates everything within the summation, which includes multiple factorials, exponential calculations, etc. There is a way to highly optimize this, is to accumulate the values according to the index of the sum.

The first one is \displaystyle545140134k+13591409. Let’s define it as the constant \displaystyle L_k, where \displaystyle k is the variable of the sum. Since \displaystyle13591409 is always here, we can define it as the base, or the initial value when the index of the sum is zero. Notice that when the index is zero, \displaystyle545140134 will be multiplied and become zero.

\displaystyle L_0=13591409

Now, as the index grows, \displaystyle545140134 will be added more and more times. Each time the index is increased by one, it will be added once. Hence, the constant is as the following:

\displaystyle L_{k+1}=L_k+545140134\text{ where }L_0=13591409

The second one is \displaystyle(-262537412640768000)^k. Let’s define it as the constant \displaystyle X_k. The number of times \displaystyle-262537412640768000 is multiplied increases by one as the index of sum increases by one. When the index is 0, the exponent is 0 and the result will be 1. Therefore, the base is 1. Now the constant is as the following:

\displaystyle X_{k+1}=X_k\cdot-262537412640768000\text{ where }X_0=1

The most difficult part to optimize is \displaystyle\frac{(6k)!}{(3k)!(k!)^3}. Since there is no addition or subtraction in this term, there will be only accumulation of multiplication. Let’s define it as the constant \displaystyle M_k. Instead of thinking how we can accumulate, let’s directly simplify \displaystyle\frac{M_{k+1}}{M_k} to find out what we need to multiply each time the index of sum is increased by one.

\displaystyle\frac{M_{k+1}}{M_k}
\displaystyle=\frac{\displaystyle\frac{(6(k+1))!}{(3(k+1))!((k+1)!)^3}}{\displaystyle\frac{(6k)!}{(3k)!(k!)^3}}
\displaystyle=\frac{(6(k+1)!(3k)!(k!)^3}{(3(k+1))!((k+1)!)^3(6k)!}
\displaystyle=\frac{(6k+6)!(3k)!(k!)^3}{(3k+3)!((k+1)!)^3(6k)!}
\displaystyle=\frac{(6k)!(6k+1)(6k+2)(6k+3)(6k+4)(6k+5)(6k+6)(3k)!(k!)^3}{(3k)!(3k+1)(3k+2)(3k+3)(k!)^3(k+1)^3(6k)!}
\displaystyle=\frac{(6k+1)(6k+2)(6k+3)(6k+4)(6k+5)(6k+6)}{(3k+1)(3k+2)(3k+3)(k+1)^3}
\displaystyle=\frac{2^3(6k+1)(3k+1)(6k+3)(3k+2)(6k+5)(3k+3)}{(3k+1)(3k+2)(3k+3)(k+1)^3}
\displaystyle=\frac{2^3(6k+1)(6k+3)(6k+5)}{(k+1)^3}
\displaystyle=\frac{(12k+2)(12k+6)(12k+10)}{(k+1)^3}

It is safe to cancel the terms, because \displaystyle k must be negative and fractional to result terms in 0, which is impossible. As a result, the constant is as the following:

\displaystyle M_{k+1}=M_k\cdot\frac{(12k+2)(12k+6)(12k+10)}{(k+1)^3}\text{ where }M_0=1

Notice that multiplication is still going on with the constant. Therefore, further optimizations can be done. Let’s define a new constant \displaystyle K. It will be used to accumulate the value of the second term in \displaystyle M_k. Although you can accumulate the value of all three terms, it will be more expensive than accumulating one only, as you will have to do more additions, and use more memory.

Storing one only is more optimized, since we can produce the value for other terms by comparing and manipulating the base or initial value of the term, provided that the addition and subtraction are cheap.

There is a pattern going on with those terms, is that the first term has the constant four smaller than the constant in the second term, and the third term has the constant four larger than the constant in the second term. Using this pattern to find out the value for the first and third term is the most optimal, for less bits are manipulated, though it won’t matter much.

\displaystyle K_{k+1}=K_k+12\text{ where }K_0=6

Now, we can change the definition of \displaystyle M_k.

\displaystyle(12k+2)(12k+6)(12k+10)
\displaystyle=(K_k-4)(K_k)(K_k+4)
\displaystyle=(K_k^2-4^2)K_k
\displaystyle=K_k^3-16K_k

Hence, the new definition is:

\displaystyle M_{k+1}=M_k\cdot\frac{K_k^3-16K_k}{(k+1)^3}\text{ where }M_0=1

Finally, the Chudnovsky Algorithm is optimized both mathematically and programmatically. For higher readability, let’s define \displaystyle C=426880\sqrt{10005}. Here is the summary of the optimized algorithm.

\displaystyle\pi=C\Bigg(\sum_{k=0}^\infty\frac{M_k\cdot L_k}{X_k}\Bigg)^{-1}

\displaystyle C=426880\sqrt{10005}
\displaystyle K_{k+1}=K_k+12\text{ where }K_0=6
\displaystyle M_{k+1}=M_k\cdot\frac{K_k^3-16K_k}{(k+1)^3}\text{ where }M_0=1
\displaystyle L_{k+1}=L_k+545140134\text{ where }L_0=13591409
\displaystyle X_{k+1}=X_k\cdot-262537412640768000\text{ where }X_0=1

Published by bernard2518141184

Leave a comment