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 which the formula was first published, David H. Bailey, Peter Borwein, and Plouffe. [cite journal
author = Bailey, David H., Borwein, Peter B., and Plouffe, Simon
year =1997 | month = April
title = On the Rapid Computation of Various Polylogarithmic Constants
journal = Mathematics of Computation
volume = 66 | issue = 218 | pages = 903–913
url = http://crd.lbl.gov/~dhbailey/dhbpapers/digits.pdf
doi = 10.1090/S0025-5718-97-00856-9 ] [ [http://groups.google.com/group/sci.math.symbolic/msg/5b7e62ce42ae0cb6 Controversy surrounding who among the three actually invented this algorithm] ]

The discovery of this formula came as a surprise. For centuries it had been assumed that there was no way to compute the "n"th digit of π without calculating all of the preceding "n" − 1 digits.

BBP-type formulae

A great many formulae that sum to other fundamental irrational constants have since been discovered. They are of the form

: $alpha\; =\; sum\_\{k\; =\; 0\}^\{infty\}\; frac\{p(k)\}\{b^kq(k)\}$

where α is a constant and "p" and "q" are polynomials in integer coefficients and b is an integer numerical base. Some formulae, such as the Leibniz formula for pi can be expressed in this form, but b is positive or negative one and thus do not display the spigot property, nor do they generally display fast convergence. In some cases, $alpha$ is an integer-scaled version of some well-known constant, such as 2π, rather than π itself.

Formulae in this form are known as BBP-type formulae. [MathWorld | urlname=BBPFormula | title=BBP Formula] Certain combinations of specific "p", "q" and "b" result in well-known constants, but there is no general algorithm for the mapping and the combinations are currently discovered via searches.

Some of the simplest formulae of this type that were well-known before BBP are

:$ln\; frac\{9\}\{10\}\; =\; -\; sum\_\{k\; =\; 1\}^\{infty\}\; frac\{1\}\{10^k\; k\}$:$ln\; 2\; =\; sum\_\{k=1\}^\{infty\}frac\{1\}\{2^k\; k\}$

Plouffe was also inspired by the related well-known formula

:$arctanfrac\{1\}\{b\}\; =\; frac\{1\}\{b\}\; -\; frac\{1\}\{b^3\; 3\}\; +\; frac\{1\}\{b^5\; 5\}\; -\; frac\{1\}\{b^7\; 7\}\; +\; frac\{1\}\{b^9\; 9\}\; +\; cdots$

A specialization of the general formula that has produced many results is

:$P(s,b,m,A)\; =\; sum\_\{k=0\}^\{infty\}left(frac\{1\}\{b^k\}sum\_\{j=1\}^\{m\}frac\{a\_j\}\{(mk+j)^s\}\; ight)$

where "s", "b" and "m" are integers and $A\; =\; (a\_1,\; a\_2,\; dots\; ,\; a\_m)$ is a vector of integers. This P function leads is a compact notation for some solutions. For the above two formulae, that would be

:$ln\; frac\{9\}\{10\}\; =\; -\; frac\{1\}\{10\}\; P(1,10,1,(1))$:$ln\; 2\; =\; frac\{1\}\{2\}\; P(1,2,1,(1))$

The search for new equalities

Using the P function mentioned above, the simplest known formula for π is for "s = 1" but "m > 1". Many now-discovered formulae are known for b as an exponent of 2 or 3 and m is an exponent of 2 or it is some other factor-rich value but where several of the terms of vector A are zero. The discovery of these fomulae involves a computer search for such linear combinations after computer the individual sums. The search procedure consists of choosing a range of parameter values for s, b and m, evaluating the sums out to many digits, and then using an integer relation finding algorithm (typically Helaman Ferguson's PSLQ algorithm) to find a vector "A" that adds up those intermediate sums to a well-known constant or perhaps to zero.

The BBP formula for π

The original BBP π summation formula was found in 1995 by Plouffe using PSLQ. It is

: $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),$

which reduces to this equivalent ratio of two polynomials:

: $pi\; =\; sum\_\{k\; =\; 0\}^\{infty\}\; frac\{1\}\{16^k\}left(frac\{120k^2\; +\; 151k\; +\; 47\}\{512k^4\; +\; 1024k^3\; +\; 712k^2\; +\; 194k\; +\; 15\}\; ight)$

and to the compact notation of π = "P"(1, 16, 8, (4, 0, 0, −2, −1,−1, 0, 0)). This formula has been shown through a rigorous and fairly simple proof to equal π. [http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf The Quest for Pi] ]

BBP digit-extraction algorithm for π

The formula yields an algorithm for extracting hexadecimal digits of π. In order to perform digit extraction first we must rewrite the formula as

: $pi\; =\; 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)\}.$

A few manipulations are required to implement a spigot algorithm using this formula. We would like to find hexadecimal digit "n" of π, so, taking the first sum we split the sum to infinity across the "n"th term

: $sum\_\{k\; =\; 0\}^\{infty\}\; frac\{1\}\{(16^k)(8k+1)\}\; =\; sum\_\{k\; =\; 0\}^\{n\}\; frac\{1\}\{(16^k)(8k+1)\}\; +\; sum\_\{k\; =\; n\; +\; 1\}^\{infty\}\; frac\{1\}\{(16^k)(8k+1)\}.$

We now multiply by 16 ^"n" so that the hexadecimal point (the divide between fractional and integer parts of the number) is in the "n"th place.

: $sum\_\{k\; =\; 0\}^\{infty\}\; frac\{16^\{n-k\{8k+1\}\; =\; sum\_\{k\; =\; 0\}^\{n\}\; frac\{16^\{n-k\{8k+1\}\; +\; sum\_\{k\; =\; n\; +\; 1\}^\{infty\}\; frac\{16^\{n-k\{8k+1\}.$

Since we only care about the fractional part of the sum, we look at our two terms and realise that only the first sum is able to produce whole numbers whereas the second sum cannot produce whole numbers since the numerator can never be larger than the denominator for "k" > "n". Therefore we need a trick to remove the whole numbers for the first sum. That trick is mod 8"k" + 1. Our sum for the first fractional part then becomes:

: $sum\_\{k\; =\; 0\}^\{n\}\; frac\{16^\{n-k\}\; mod\; 8k+1\}\{8k+1\}\; +\; sum\_\{k\; =\; n\; +\; 1\}^\{infty\}\; frac\{16^\{n-k\{8k+1\}.$

Notice how the modulo operator always guarantees that only the fractional sum will be kept. To calculate 16 ^{"n" − "k"} mod 8"k" + 1 quickly and efficiently, use the modular exponentiation algorithm. When the running product becomes greater than one, take the modulo just as you do for the running total in each sum.

Now to complete the calculation you must apply this to each of the four sums in turn. Once this is done, take the four summations and put them back into the sum to π.

: $4\; Sigma\_1\; -\; 2\; Sigma\_2\; -\; Sigma\_3\; -\; Sigma\_4.\; ,!$

Since only the fractional part is accurate, extracting the wanted digit requires that one removes the integer part of the final sum and multiplies by 16 to "skim off" the hexadecimal digit at this position (in theory the next few digits up to the accuracy of the calculations used would also be accurate).

This process is similar to performing long multiplication, but only having to perform the summation of some middle columns. While there are some carries that are not counted, computers usually perform arithmetic for many bits (32 or 64) and they round and we are only interested in the most significant digit(s). There is a vanishingly small possibility that a particular computation will be akin to failing to add a small number (e.g. 1) to the number 999999999999999 and that the error will propagate to the most significant digit, but being near this situation is obvious in the final value produced.

This algorithm has proven popular and has many software implementations. [ [http://www.nedprod.com/programs/portable/Pi/index.html A C++ implementation of the BBP algorithm for π (portable, SSE2 and OpenMP versions)] ] [ [http://en.literateprograms.org/Pi_with_the_BBP_formula_(Python) A Python implementation of the BBP algorithm for π] ] [ [http://www.314159265358979323846264338327950288419716939937510.net A Ruby implementation of the BBP algorithm for π] ]

Advantages of the BBP algorithm

This algorithm computes π without requiring custom data types having thousands or even millions of digits. The method calculates the "n"th digit "without" calculating the first "n" − 1 digits, and can use small, efficient data types.

The algorithm is the fastest way to compute the "n"th digit (or a few digits in a neighborhood of the "n"th), but π-computing algorithms using large data types remain faster when the goal is to compute all the digits from 1 to "n".

Generalizations

D.J. Broadhurst provides a generalization of the BBP algorithm [D.J. Broadhurst, " [http://arxiv.org/abs/math.CA/9803067 Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5)] ", (1998) "arXiv" math.CA/9803067] that may be used to compute a number of other constants in nearly linear time and logarithmic space. Explicit results are given for Catalan's constant, $pi^3$, $log^32$, Apery's constant $zeta(3)$ (where $zeta(x)$ is the Riemann zeta function), $pi^4$, $log^42$, $log^52$, $zeta(5)$, and various products of powers of $pi$ and $log2$. These results are obtained primarily by the use of polylogarithm ladders.

ee also

*Computing π
*Experimental mathematics
*Bellard's formula

References