Fórmula de Bailey-Borwein-Plouffe


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]

Contenido

La fórmula BBP

\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}\right)

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

Uso de la fórmula para calcular los decimales del número π

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 16Nπ. 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 16N, la coma se desplaza N lugares hacia la derecha. El problema original se reduce al cálculo del primer dígito de 16Nπ. Usando la fórmula BBP:

16^{N}\pi = \sum_{k=0}^\infty 16^{N-k}\left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)

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 S_N(a)=\sum_{k=0}^\infty \frac{16^{N-k}}{8k+a}. El cálculo de los primeros dígitos hexadecimales de SN(a) permitirá obtener los de 16Nπ, a través de la relación:

16^{N}\pi = 4\ S_N(1) - 2 S_N(4) - S_N(5) - S_N(6)

Descomponiendo la suma SN(a) en dos:

S_N(a)=\sum_{k=0}^\infty \frac{16^{N-k}}{8k+a}=\sum_{k=0}^{N-1} \frac{16^{N-k}}{8k+a}+\sum_{k=N}^ \infty \frac{16^{N-k}}{8k+a} = A_N(a) + B_N(a)

se pueden calcular AN(a) y BN(a) en forma independiente.

Cálculo de BN(a)

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.

  • En efecto, el primer término de la suma es: b_N=\frac{1}{8N+a}. Si se busca el enésimo dígito hexadecimal de π (N = 1 000 000 000 por ejemplo), el primer término es mucho menor que 1.
  • Además, cada término tiene un cero más a la derecha de la coma que el precedente, porque para k ≥ N, bk > 16 bk+1:
\frac{b_k}{b_{k+1}}=\frac{16^{N-k}}{16^{N-(k+1)}}\frac{8(k+1)+a}{8k+a}=16\left( 1+\frac{8}{8k+a}\right) \longrightarrow 16^{+}.

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

B_N\ =\ \ 0,**********.\ .\ .\ .\

+\ 0,0*********.\ .\ .\
+\ 0,00********.\ .\ .\
+\ 0,000********.\ .\ .\

Por lo tanto, para obtener BN(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: B_N'(a)=\sum_{k=N}^{N+P+10} \frac{16^{N-k}}{8k+a}

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

Cálculo de AN(a)

A_N(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}}{8k+a}

El problema para el cálculo de AN(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 \frac{16^{N-k}}{8k+a}. Para ello realizamos la división entera de 16N-k por 8k+a:

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

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

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

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

Así, se calcula: 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, 16Nk(mod 8k + a) se calcula rápidamente (con un tiempo de ejecución de O(log2(N-k)).

Conclusión

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

\pi_N=4\ S_N'(1)-2\ S_N'(4)-S_N'(5)-S_N'(6)

con 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}.

Complejidad de este método

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

Complejidad temporal

  • Bn'(a) se calcula en complejidad lineal (O(1))
  • An'(a) : utilizando el método de la exponenciación binaria, sus términos se calculan en O(log2(n)). Así la suma de n términos, An'(a), se calcula en O(n*log2(n)).

Así Sn'(a) se calcula en O(1)+O(n*log2(n))=O(n*log2(n)). Finalmente, πn se calcula en 4*O(n*log2(n))=O(n*log2(n)).

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

Complejidad espacial

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

Formulas derivadas

Simon Plouffe

Fórmula original : \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}\right)

\forall r \in \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}\right)

\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}\right)

\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}\right)

\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}\right)

\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}\right)

La última fórmula permite hallar dígitos aislados de π2 en base 3 ó 9.

Viktor Adamchick y Stan Wagon (1997)

\pi\ =\ \sum_{i=0}^\infty
\frac{(-1)^i}{4^{i}}\left(\frac{2}{4i+1}+\frac{2}{4i+2}+\frac{1}{4i+3}\right)

Fabrice Bellard

\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}\right)

Los récords

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).

  • 7 de octubre de 1996 (Fabrice Bellard): dígito número 400 mil millones en base 2
  • septiembre de 1997 (Fabrice Bellard) : billonésimo dígito en base 2
  • febrero de 1999 (Colin Percival) : dígito número 40 billones en base 2
  • 2001: dígito número 4000 billones en base 2
  • 2011: dígito número 10 billones en base 10

Cálculo del enésimo decimal

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(n3*log2(n)) pierde su utilidad en la práctica. Fabrice Bellard ha mejorado el algoritmo para alcanzar un nivel de complejidad en O(n2), pero no es suficiente para competir con los métodos convencionales que calculan todos los decimales.

Anexo : demostración de la fórmula BBP

Se demuestra primero el siguiente resultado general:

\forall n \in \mathbb{N}: \sum_{k=0}^\infty \frac{1}{16^k(8k+n)}=\sqrt{2}^{n}\sum_{k=0}^\infty \frac{\left(\frac{1}{\sqrt{2}}\right) ^{8k+n}}{8k+n}
=\sqrt{2}^{n}\sum_{k=0}^\infty \left[\frac{x^{8k+n}}{8k+n}\right]_0^{\frac{1}{\sqrt{2}}}
=\sqrt{2}^{n}\sum_{k=0}^\infty \left(\int_{0}^{\frac{1}{\sqrt{2}}} x^{n+8k-1}dx \right)
=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \left(\sum_{k=0}^\infty x^{n+8k-1}\right)dx
=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \bigg(x^{n-1}\underbrace{\sum_{k=0}^\infty \left(x^8\right)^k} \bigg)dx
=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \bigg(x^{n-1} \frac{1}{1-x^8}\bigg)dx


por lo tanto: \sum_{k=0}^\infty \frac{1}{16^k(8k+n)}=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{n-1}}{1-x^8}dx


Aplicando este resultado a la fórmula BBP:

\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}\right)
= 4\sum_{k=0}^\infty \frac{1}{16^k(8k+1)}-2\sum_{k=0}^\infty \frac{1}{16^k(8k+4)}-\sum_{k=0}^\infty \frac{1}{16^k(8k+5)}-\sum_{k=0}^\infty \frac{1}{16^k(8k+6)}
= 4\left(\sqrt{2}^{1}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{1-1}}{1-x^8}dx \right)-2\left(\sqrt{2}^{4}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{4-1}}{1-x^8}dx \right)-\left(\sqrt{2}^{5}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{5-1}}{1-x^8}dx \right)-\left(\sqrt{2}^{6}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{6-1}}{1-x^8}dx \right)
= \int_{0}^{\frac{1}{\sqrt{2}}} \frac{4\sqrt{2}-8x^3-4\sqrt{2}x^4-8x^5}{1-x^8}dx

Sustituyendo y=\sqrt{2}x:

= \int_{0}^{1} \left(\frac{4\sqrt{2}-\frac{8}{\sqrt{2}^3}y^3-\frac{4\sqrt{2}}{\sqrt{2}^4}y^4-\frac{8}{\sqrt{2}^5}y^5}{1-\frac{1}{\sqrt{2}^8}y^8}\right)\frac{dy}{\sqrt{2}}
= 16 \int_{0}^{1} \frac{4-2y^3-y^4-y^5}{16-y^8}dy
= 16 \int_{0}^{1} \frac{y-1}{(y^2-2y+2)(y^2-2)}dy

Descomponiendo en fracciones simples:

= \int_{0}^{1}\left( \frac{8-4y}{y^2-2y+2}+\frac{4y}{y^2-2}\right) dy
= -2\int_{0}^{1}\frac{2y-2}{y^2-2y+2}dy+4\int_{0}^{1}\frac{1}{1+(y-1)^2}dy+2\int_{0}^{1}\frac{2y}{y^2-2}dy
= -2\left[ln(y^2-2y+2)\right]_0^1+4\left[arctan(y-1)\right]_0^1+2\left[ln(2-y^2)\right]_0^1
= -2\ ln(1)+2\ ln(2)+4\ arctan(0)-4\ arctan(-1)+2\ ln(1)-2\ ln(2)
= 4\ arctan(1)
= \ \pi

Véase también

Referencias

  1. The Quest for Pi (en inglés)

Wikimedia foundation. 2010.

Mira otros diccionarios:

  • Fórmula de Bailey-Borwein-Plouffe — La fórmula de Bailey Borwein Plouffe (BBP formula) está relacionada con el cálculo del número p. Fue descubierta en 1995 por Simon Plouffe. Junto con esta fórmula se han ido descubriendo muchas otras que poseen la siguiente forma: y que se… …   Enciclopedia Universal

  • Bailey–Borwein–Plouffe formula — The Bailey–Borwein–Plouffe formula (BBP formula) provides a spigot algorithm for the computation of the n th binary digit of π. This summation formula was discovered in 1995 by Simon Plouffe. The formula is named after the authors of the paper in …   Wikipedia

  • Borwein's algorithm — In mathematics, Borwein s algorithm is an algorithm devised by Jonathan and Peter Borwein to calculate the value of 1/ pi;. The most prominent and oft used one is explained under the first section.Borwein s algorithmStart out by setting: a 0 = 6… …   Wikipedia

  • David H. Bailey — For other people named David Bailey, see David Bailey (disambiguation). David Bailey (2010) David Harold Bailey (born 1948) is a mathematician and computer scientist. He received his B.S. in mathematics from Brigham Young University in 1972 and… …   Wikipedia

  • Simon Plouffe — is a Quebec mathematician born on June 11 1956 in Saint Jovite, Quebec. He discovered the formula for the BBP algorithm (the Bailey–Borwein–Plouffe formula ) which permits the computation of the n th binary digit of pi;, in 1995. Plouffe is also… …   Wikipedia

  • Peter Borwein — Peter Benjamin Borwein (St. Andrews, Scotland, 1953) is a Canadian mathematician, co developer of an algorithm for calculating π to the n th digit, PiHex, co discoverer of the trillionth, four trillionth, 40th trillionth, and quadrillionth digits …   Wikipedia

  • Simon Plouffe — (11 juin 1956 à Saint Jovite, Québec, Canada) est un mathématicien canadien. Sommaire 1 Travaux 2 Récompenses 3 Voir aussi …   Wikipédia en Français

  • Bellard's formula — Bellard s formula, as used by PiHex, the now completed distributed computing project, is used to calculate the n th digit of π in base 2. It is a faster version (about 43% faster [ [http://oldweb.cecm.sfu.ca/projects/pihex/credits.html PiHex… …   Wikipedia

  • Pi — This article is about the number. For the Greek letter, see Pi (letter). For other uses, see Pi (disambiguation). The circumference of a ci …   Wikipedia

  • Approximations of π — Timeline of approximations for pi …   Wikipedia