Desde la antigüedad la fascinación por calcular las cifras decimales de $\pi$ motivó a los matemáticos. Arquímedes, Ptolomeo, Liu Hui, Zu Chongzhi, al-Kashi, Viète, Leibniz, Newton, … y un largo etcétera pusieron su empeño en encontrar procedimientos que abreviasen el cómputo de $\pi$, o los obtuviesen con mayor precisión. La fórmula más precisa que se consiguió antes de la eclosión de los computadores, la descubrió el genio indio Srinivasa Ramanujan(Ramanujan y el número pi):
$${\frac {1}{\pi }}={\frac {2{\sqrt {2}}}{9801}}\sum _{k=0}^{\infty }{\frac {(4k)!(1103+26390k)}{(k!)^{4}396^{4k}}}$$
que obtiene 8 decimales exactos de $\pi$ en cada iteración. Pero estamos en la era de los computadores, ¿cómo lo hacemos ahora?
En 1984, los hermanos Jonathan y Peter Borwein publican una iteración de convergencia de orden 4: dados
$${\begin{aligned}a_{0}&=2{\big (}{\sqrt {2}}-1{\big )}^{2}\\y_{0}&={\sqrt {2}}-1\end{aligned}}$$
entonces,
$${\begin{aligned}y_{k+1}&={\frac {1-(1-y_{k}^{4})^{1/4}}{1+(1-y_{k}^{4})^{1/4}}}\\a_{k+1}&=a_{k}(1+y_{k+1})^{4}-2^{2k+3}y_{k+1}(1+y_{k+1}+y_{k+1}^{2})\end{aligned}}$$
La sucesión $a_k$ converge a $1/\pi$. Los dos hermanos habían encontrado alguna más, con orden de convergencia mayores. En 1986, David Harold Bailey utilizó una de ellas para desarrollar un programa de cálculo de $\pi$, escrito para la NASA, como test para detectar problemas en la Cray-2[1].
En 1995 Peter Borwein y Simon Plouffe encuentran una forma de calcular dígitos en binario del logaritmo de 2, empezando en una cifra arbitraria. Plouffe se da cuenta que tiene un filón y produce la fórmula:
$$\pi =\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)\right]$$
El estudio de esta fórmula y su aplicación para calcular $\pi$ lo publican David H. Bailey, Peter Borwein, and Simon Plouffe en 1997 «The BBP Algorithm for Pi» (PDF).
Con esta fórmula, uno puede derivar un algoritmo que compute dígitos de $\pi$ en una posición arbitraria y, además, en base hexadecimal. En 1997, Frabrice Bellard de INRIA, calculó 152 dígitos binarios de $\pi$ empezando en la billonésima posición digital binaria. El cálculo tardó 12 días en 20 estaciones de trabajo que trabajan en paralelo a través de Internet. Repetiría la proeza en último día del año 2009, con 9 Desktop PCs, Core i7 CPU a 2.93 GHz, durante 131 días, para obtener más de $2\times 10^{12}$ decimales.
En los años 80 del pasado siglo, Yasumasa Kanada también trabajaba en el cálculo de decimales de $\pi$, utilizando HITAC S-820/80, con procesador vectorial, llegando a obtener 1.073.740.799 decimales. Pero los hermanos Chudnovsky no le iban a la zaga, utilizando un IBM 3090, y la fórmula:
$${\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320^{3})^{k+1/2}}}.\!
$$
El algoritmo dado por los hermanos Chudnovsky se usó en diciembre de 2013 para alcanzar los 12.1 billiones de dígitos[2].
Así que Kanada, viendo como estaban los competidores, se escurrió el cerebelo y desarrolló Super PI, un programa (ejecutable en Windows) que permitía calcular un número determinado de decimales a partir de uno dado; es decir, el mismo sistema que se le ocurrió a Plouffe. Super PI utiliza el algoritmo de Gauss-Legendre:
Primero inicializamos: $a_{0}=1\qquad b_{0}={\frac {1}{{\sqrt {2}}}}\qquad t_{0}={\frac {1}{4}}\qquad p_{0}=1.$
Segundo, iteramos hasta que la diferencia de $a_n$ y $b_n$ sea de la precisión deseada:
$${\begin{aligned}a_{{n+1}}&={\frac {a_{n}+b_{n}}{2}},\\b_{{n+1}}&={\sqrt {a_{n}b_{n}}},\\t_{{n+1}}&=t_{n}-p_{n}(a_{{n}}-a_{{n+1}})^{2},\\p_{{n+1}}&=2p_{n}.\end{aligned}}$$
Por último aproximamos
$$\pi \approx {\frac {(a_{{n+1}}+b_{{n+1}})^{2}}{4t_{{n+1}}}}.$$
Con este algoritmo Kanada batiría el récord en 2002, 1.241.100.000.000 dígitos[3]
Super Pi se vio pronto superado por el programa y-cruncher de Alexander Yee, con el que Shigeru Kondo comenzó a desbordar récord, desde agosto de 2010, hasta alcanzar los 12.1 billones de digitos – Diciembre de 2013. En la actualidad y-cruncher mantiene el récord, obtenido en octubre de 2014, con 13.3 billones de dígitos. El programa utiliza la fórmula de los hermanos Chudnovsky. El último récord se consiguió con 2 x Xeon E5-4650L @ 2.6 GHz, 192 GB DDR3 @ 1333 MHz y 24 x 4 TB + 30 x 3 TB.
Como veis, utilizando una fórmula y paralelizando sus desarrollos podemos encontrar las cifras de $\pi$ que deseemos, sólo es cuestión de tiempo. Por ejemplo, en A very simple simulation program: Parallel computation of PI, encontráis un programa sencillo que se basa en una aproximación por Simpson de la integral
$$\int_0^1\frac{4}{1+x^2}dx. $$