Fórmula de Bailey-Borwein-Plouffe

La fórmula de Bailey-Borwein-Plouffe (o fórmula BBP) permite calcular el enésimo dígito de π en base 2 (o 16) sin necesidad de hallar los precedentes, de una manera rápida y utilizando muy poco espacio de memoria en la computadora. Simon Plouffe junto con David Bailey y Peter Borwein hallaron esta fórmula el 19 de septiembre de 1995 usando un programa informático llamado PSLQ que busca relaciones entre números enteros.^[1]​

${displaystyle pi =sum _{k=0}^{infty }{frac {1}{16^{k}}}left({frac {4}{8k+1}}-{frac {2}{8k+4}}-{frac {1}{8k+5}}-{frac {1}{8k+6}} ight)}$

La demostración de esta fórmula se encuentra más abajo.

A continuación se muestra el cálculo del enésimo dígito hexadecimal de π.

Primero se debe observar que el dígito ubicado en la posición N+1 de π en base 16 es el mismo que el primer dígito hexadecimal de 16^Nπ. En efecto, como en la base 10, multiplicar un número en base 16 por 16 equivale a desplazar la coma decimal un lugar hacia la derecha. De esta manera, multiplicando por 16^N, la coma se desplaza N lugares hacia la derecha. El problema original se reduce al cálculo del primer dígito de 16^Nπ. Usando la fórmula BBP:

El cálculo de los primeros dígitos hexadecimales a la derecha de la coma de este número no es sencillo por dos razones: el número es muy grande y la suma es infinita.

Supongamos que ${displaystyle S_{N}(a)=sum _{k=0}^{infty }{frac {16^{N-k}}{8k+a}}}$. El cálculo de los primeros dígitos hexadecimales de S_N(a) permitirá obtener los de 16^Nπ, a través de la relación:

Descomponiendo la suma S_N(a) en dos:

se pueden calcular A_N(a) y B_N(a) en forma independiente.

${displaystyle B_{N}(a)=sum _{k=N}^{infty }{frac {16^{N-k}}{8k+a}}}$

Aunque se trata de una suma infinita, es muy fácil de calcular, porque sus términos son pequeños y decrecen rápidamente.

Finalmente, la suma B_N(a) es de la forma (en el peor caso):

${displaystyle B_{N} = 0,**********. . . . }$

Por lo tanto, para obtener B_N(a) con una precisión de P cifras detrás de la coma, es suficiente calcular los P primeros términos de la suma, agregándose algunos términos más para evitar errores que aparecen al realizar cálculos con valores aproximados.

Así, se calcula: ${displaystyle B_{N}'(a)=sum _{k=N}^{N+P+10}{frac {16^{N-k}}{8k+a}}}$

Como esta suma solo posee una pequeña cantidad de términos, el tiempo que insume esta operación es insignificante para una computadora.

${displaystyle A_{N}(a)=sum _{k=0}^{N-1}{frac {16^{N-k}}{8k+a}}}$

El problema para el cálculo de A_N(a) es que los primeros términos son muy grandes (N cifras de base 16 antes de la coma). Sin embargo, al igual que las primeras cifras detrás de la coma, no importa la parte entera, que también es grande. Por lo tanto, puede "eliminarse" usando aritmética modular.

Toda la dificultad se reduce a hallar la parte fraccionaria de ${displaystyle {frac {16^{N-k}}{8k+a}}}$. Para ello realizamos la división entera de 16^N-k por 8k+a:

${displaystyle exists qin mathbb {Z} ,exists r<8k+a, 16^{N-k}=q(8k+a)+r}$

Así ${displaystyle {frac {16^{N-k}}{8k+a}}=q+{frac {r}{8k+a}}}$

${displaystyle {frac {r}{8k+a}}}$ es menor que 1, por lo tanto, es la parte fraccional de ${displaystyle {frac {16^{N-k}}{8k+a}}}$.

Y ${displaystyle {frac {r}{8k+a}}={frac {16^{N-k}{pmod {8k+a}}}{8k+a}}}$

Así, se calcula: ${displaystyle A_{N}'(a)=sum _{k=0}^{N-1}{frac {16^{N-k}{pmod {8k+a}}}{8k+a}}}$.

Utilizando el método de la exponenciación binaria, ${displaystyle 16^{N-k}{pmod {8k+a}}}$ se calcula rápidamente (con un tiempo de ejecución de O(log₂(N-k)).

Al fin y al cabo, para obtener los primeros dígitos de π base 16 (o 2), se deben calcular los primeros dígitos de:

${displaystyle pi _{N}=4 S_{N}'(1)-2 S_{N}'(4)-S_{N}'(5)-S_{N}'(6)}$

con ${displaystyle S_{N}'(a)=sum _{k=0}^{N-1}{frac {16^{N-k}{pmod {8k+a}}}{8k+a}}+sum _{k=N}^{N+P+10}{frac {16^{N-k}}{8k+a}}}$.

Para calcular el enésimo dígito de π en base 16 (o el 4n-ésimo dígito en base 2):

Así S_n'(a) se calcula en O(1)+O(n*log₂(n))=O(n*log₂(n)). Finalmente, π_n se calcula en 4*O(n*log₂(n))=O(n*log₂(n)).

Por lo tanto el tiempo de cálculo es proporcional a n*log₂(n), es decir, casi lineal.

La complejidad en el uso de memoria es constante, ya que solo se realizan sumas sucesivas de pequeños números (con una precisión de unos diez decimales es suficiente).

Fórmula original : ${displaystyle pi = sum _{i=0}^{infty }{frac {1}{16^{i}}}left({frac {4}{8i+1}}-{frac {2}{8i+4}}-{frac {1}{8i+5}}-{frac {1}{8i+6}} ight)}$

${displaystyle forall rin mathbb {C} pi = sum _{i=0}^{infty }{frac {1}{16^{i}}}left({frac {4+8r}{8i+1}}-{frac {8r}{8i+2}}-{frac {4r}{8i+3}}-{frac {2+8i}{8i+4}}-{frac {1+2r}{8i+5}}-{frac {1+2r}{8i+6}}+{frac {r}{8i+7}} ight)}$

${displaystyle pi {sqrt {2}} = sum _{i=0}^{infty }{frac {(-1)^{i}}{8^{i}}}left({frac {4}{6i+1}}+{frac {1}{6i+2}}+{frac {1}{6i+3}} ight)}$

${displaystyle pi ^{2} = sum _{i=0}^{infty }{frac {1}{16^{i}}}left({frac {16}{(8i+1)^{2}}}-{frac {16}{(8i+2)^{2}}}-{frac {8}{(8i+3)^{2}}}-{frac {16}{(8i+4)^{2}}}-{frac {4}{(8i+5)^{2}}}-{frac {4}{(8i+6)^{2}}}-{frac {2}{(8i+7))^{2}}} ight)}$

${displaystyle pi ^{2} = {frac {9}{8}} sum _{i=0}^{infty }{frac {1}{64^{i}}}left({frac {16}{(6i+1)^{2}}}-{frac {24}{(6i+2)^{2}}}-{frac {8}{(6i+3)^{2}}}-{frac {6}{(6i+4)^{2}}}-{frac {1}{(6i+5)^{2}}} ight)}$

${displaystyle pi ^{2} = {frac {2}{27}} sum _{i=0}^{infty }{frac {1}{729^{i}}}left({frac {243}{(12i+1)^{2}}}-{frac {405}{(12i+2)^{2}}}-{frac {81}{(12i+4)^{2}}}-{frac {27}{(12i+5)^{2}}}-{frac {72}{(12i+6)^{2}}}-{frac {9}{(12i+7)^{2}}}-{frac {9}{(12i+8)^{2}}}-{frac {5}{(12i+10)^{2}}}+{frac {1}{(12i+11)^{2}}} ight)}$

La última fórmula permite hallar dígitos aislados de ${displaystyle pi ^{2}}$ en base 3 o 9.

${displaystyle pi = sum _{i=0}^{infty }{frac {(-1)^{i}}{4^{i}}}left({frac {2}{4i+1}}+{frac {2}{4i+2}}+{frac {1}{4i+3}} ight)}$

${displaystyle pi = {frac {1}{64}} sum _{i=0}^{infty }{frac {(-1)^{i}}{2^{10i}}}left(-{frac {32}{4i+1}}-{frac {1}{4i+3}}+{frac {256}{10i+1}}-{frac {64}{10i+3}}-{frac {4}{10i+5}}-{frac {4}{10i+7}}+{frac {1}{10i+9}} ight)}$

Para poder comparar, hasta el año 2008 se habían obtenido los primeros 1,241 billones de decimales de π (aproximadamente 4,123 billones de bits).

Actualmente, no existe ninguna fórmula eficaz para hallar el enésimo decimal de π en base 10. Simon Plouffe ha desarrollado en diciembre de 1996, a partir de una serie muy antigua que calcula π basado en los coeficientes de binomio de Newton, un método para calcular cifras aisladas base 10, pero debido a su complejidad O(n³*log₂(n)) pierde su utilidad en la práctica. Fabrice Bellard ha mejorado el algoritmo para alcanzar un nivel de complejidad en O(n²), pero no es suficiente para competir con los métodos convencionales que calculan todos los decimales.

Se demuestra primero el siguiente resultado general:

Aplicando este resultado a la fórmula BBP:

Sustituyendo ${displaystyle y={sqrt {2}}x}$:

Descomponiendo en fracciones simples:

---