Bailey-Borwein-Plouffe formula

From Wikipedia, the free encyclopedia

The Bailey-Borwein-Plouffe formula (BBP formula) provides a spigot algorithm for the computation of the nth 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.[1]

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

Contents

[edit] 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, α is a scaled version of some well-known constant, such as 2π, rather than π itself.

Formulae in this form are known as BBP-type formulae.[2] 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

\arctan\frac{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} ...

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

P(s,b,n,A) = \sum_{k=0}^{\infty}\left(\frac{1}{b^k}\sum_{j=1}^{n}\frac{a_j}{(nk+j)^s}\right)

where s, b and n are integers and A = (a1,a2,...,an) is a vector of integers. The simplest π formula is for s=1. Many now-discovered formulae are known for b as an exponent of 2 or 3 and n is an exponent of 2 or it is some other factor-rich value. For many now-known formulae, several of the terms of vector A are zero. 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 consists of choosing parameter values for s, b and n, evaluating the n terms of 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.

[edit] 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}\right),

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

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 proof to equal π.[3]

[edit] BBP algorithm

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

Now we would like to find hexadecimal digit n of π, so, taking the first sum we split the sum to infinity across the nth 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 16n so that the hexadecimal point (the divide between fractional and integer parts of the number) is in the nth 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 8k+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 8k+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.

[edit] Advantages of the BBP algorithm

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

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

[edit] Generalizations

D.J. Broadhurst provides a generalization of the BBP algorithm[4] 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, π3, log32, Apery's constant ζ(3) (where ζ(x) is the Riemann zeta function), π4, log42, log52, ζ(5), and various products of powers of π and log2. These results are obtained primarily by the use of polylogarithm ladders.

[edit] See also

[edit] References

  1. ^ Bailey, David H., Borwein, Peter B., and Plouffe, Simon (April 1997). "On the Rapid Computation of Various Polylogarithmic Constants". Mathematics of Computation 66 (218): 903–913. 
  2. ^ Eric W. Weisstein, BBP Formula at MathWorld.
  3. ^ The Quest for Pi
  4. ^ D.J. Broadhurst, "Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5)", (1998) arXiv math.CA/9803067

[edit] External links