View source for Bailey–Borwein–Plouffe formula

{{short description|Formula for computing the nth base-16 digit of π}}
The '''Bailey–Borwein–Plouffe formula''' ('''BBP formula''') is a formula for [[pi|{{pi}}]]. It was discovered in 1995 by [[Simon Plouffe]] and is named after the authors of the article in which it was published, [[David H. Bailey (mathematician)|David H. Bailey]], [[Peter Borwein]], and Plouffe.<ref name="BBP">{{cite journal
 | last1 = Bailey | first1 = David H. | last2 = Borwein | first2 = Peter B. | last3 = Plouffe | first3 = Simon
 | year = 1997
 | 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   PDF is free at doi link, and is better formatted. -->
 | doi = 10.1090/S0025-5718-97-00856-9
 | mr = 1415794
 | doi-access = free| hdl = 2060/19970009337| hdl-access = free}}</ref> The formula is:
:<math>  \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]</math>
The BBP formula gives rise to a [[spigot algorithm]] for computing the ''n''th [[Hexadecimal|base-16]] (hexadecimal) digit of {{pi}} (and therefore also the ''4n''th [[binary digit]] of {{pi}}) without computing the preceding digits. This does ''not'' compute the ''n''th decimal digit of {{pi}} (i.e., in base 10).<ref>{{cite web
 | last   = Gourdon
 | first  = Xavier
 | title  = N-th Digit Computation
 | date   = 12 February 2003
 | url    = http://www.plouffe.fr/simon/articles/nthdigit.pdf
 | accessdate = 4 November 2020}}</ref> But another formula discovered by Plouffe in 2022 allows extracting the ''n''th digit of {{pi}} in decimal.<ref>{{cite arXiv|eprint=2201.12601 |last1=Plouffe |first1=Simon |title=A formula for the $n^&#123;\rm th&#125;$ decimal digit or binary of $π$ and powers of $π$ |date=2022 |class=math.NT }}</ref> BBP and BBP-inspired algorithms have been used in projects such as [[PiHex]]<ref>{{cite web|title=PiHex Credits|url=http://wayback.cecm.sfu.ca/projects/pihex/credits.html|website=Centre for Experimental and Constructive Mathematics|publisher=Simon Fraser University|accessdate=30 March 2018|archiveurl=https://web.archive.org/web/20170610094408/http://wayback.cecm.sfu.ca/projects/pihex/credits.html|archivedate=2017-06-10|url-status=live|date=March 21, 1999}}</ref> for calculating many digits of {{pi}} using distributed computing. The existence of this formula came as a surprise. It had been widely believed that computing the ''n''th digit of {{pi}} is just as hard as computing the first ''n'' digits.<ref name="BBP" />

Since its discovery, formulas of the general form:
: <math>\alpha = \sum_{k = 0}^\infty \left[ \frac{1}{b^k} \frac{p(k)}{q(k)} \right]</math>
have been discovered for many other [[irrational number]]s <math>\alpha</math>, where <math>p(k)</math> and <math>q(k)</math> are [[polynomial]]s with integer coefficients and <math>b \geq 2</math> is an integer [[radix|base]].
Formulas of this form are known as '''BBP-type formulas'''.<ref name="Ref_a">{{MathWorld| title=BBP Formula| urlname=BBPFormula}}</ref> Given a number <math>\alpha</math>, there is no known systematic algorithm for finding appropriate <math>p(k)</math>, <math>q(k)</math>, and <math>b</math>; such formulas are discovered [[experimental mathematics|experimentally]].

== Specializations ==
A specialization of the general formula that has produced many results is:
: <math> P(s, b, m, A) = \sum_{k=0}^{\infty}\left[\frac{1}{b^k} \sum_{j=1}^m \frac{a_j}{(mk+j)^s} \right],</math>
where ''s'', ''b'', and ''m'' are integers, and <math>A = (a_1, a_2, \dots, a_m)</math> is a [[sequence]] of integers.
The ''P'' function leads to a compact notation for some solutions. For example, the original BBP formula: 
:<math> \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]</math>
can be written as:
: <math> \pi = P\bigl(1, 16, 8, (4, 0, 0, -2, -1, -1, 0, 0) \bigr).</math>

=== Previously known BBP-type formulae ===
Some of the simplest formulae of this type that were well known before BBP and for which the ''P'' function leads to a compact notation, are:
: <math>\begin{align}
 \ln\frac{10}{9} &= \frac{1}{10} + \frac{1}{200} + \frac{1}{3\ 000} + \frac{1}{40\,000} + \frac{1}{500\,000} + \cdots \\
  &= \sum_{k=1}^\infty \frac{1}{10^k \cdot k} = \frac{1}{10} \sum_{k=0}^\infty \left[ \frac{1}{10^k} \left( \frac{1}{k+1} \right) \right] \\
  &= \frac{1}{10} P\bigl(1, 10, 1, (1) \bigr),
\end{align}</math>
: <math>\begin{align}
 \ln 2 &= \frac{1}{2} + \frac{1}{2 \cdot 2^2} + \frac{1}{3 \cdot 2^3} + \frac{1}{4 \cdot 2^4} + \frac{1}{5 \cdot 2^5} + \cdots \\
  &= \sum_{k=1}^\infty \frac{1}{2^k \cdot k} = \frac{1}{2} \sum_{k=0}^\infty \left[ \frac{1}{2^k} \left( \frac{1}{k + 1} \right) \right] \\
  &= \frac{1}{2} P\bigl( 1, 2, 1, (1) \bigr).
\end{align}</math>

(In fact, this identity holds true for ''a'' > 1:
:<math>\ln\frac{a}{a-1} = \sum_{k=1}^\infty \frac{1}{a^k \cdot k}</math>.)

Plouffe was also inspired by the [[Inverse tangent#Infinite series|arctan power series]] of the form (the ''P'' notation can be also generalized to the case where ''b'' is not an integer):
: <math>\begin{align}
 \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} + \cdots \\
  &= \sum_{k=1}^\infty \left[ \frac{1}{b^{k}} \frac{\sin\frac{k\pi}{2}}{k} \right]
   = \frac{1}{b} \sum_{k=0}^\infty \left[ \frac{1}{b^{4k}} \left( \frac{1}{4k+1} + \frac{-b^{-2}}{4k+3} \right) \right] \\
  &= \frac{1}{b} P\left( 1, b^4, 4, \left(1, 0, -b^{-2}, 0\right) \right).
\end{align}</math>

=== The search for new equalities ===
Using the ''P'' function mentioned above, the simplest known formula for {{pi}} is for ''s''&nbsp;=&nbsp;1, but ''m''&nbsp;>&nbsp;1. Many now-discovered formulae are known for ''b'' as an exponent of 2 or 3 and ''m'' as an exponent of 2 or it some other factor-rich value, but where several of the terms of sequence ''A'' are zero. The discovery of these formulae involves a computer search for such linear combinations after computing 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 algorithm|integer relation-finding algorithm]] (typically [[Helaman Ferguson]]'s [[PSLQ algorithm]]) to find a sequence ''A'' that adds up those intermediate sums to a well-known constant or perhaps to zero.

=== The BBP formula for {{pi}} ===
The original BBP {{pi}} summation formula was found in 1995 by Plouffe using [[integer relation algorithm|PSLQ]]. It is also representable using the ''P'' function:
: <math>\begin{align}
 \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] \\
  &= P\bigl( 1, 16, 8, (4, 0, 0, -2, -1, -1, 0, 0) \bigr),
\end{align}</math>
which also reduces to this equivalent ratio of two polynomials:
: <math>\pi = \sum_{k = 0}^\infty \left[ \frac{1}{16^k} \left( \frac{120k^2 + 151k + 47}{512k^4 + 1024k^3 + 712k^2 + 194k + 15} \right) \right].</math>

This formula has been shown through a fairly simple proof to equal {{pi}}.<ref name=QuestForPi>{{cite journal | last1 = Bailey | first1 = David H. | last2 = Borwein | first2 = Jonathan M. | last3 = Borwein | first3 = Peter B. | last4 = Plouffe | first4 = Simon | year = 1997 | journal = [[Mathematical Intelligencer]] | volume = 19 | issue = 1 | pages = 50–57 | doi=10.1007/BF03024340 | mr=1439159 | title=The quest for pi| s2cid = 14318695 }}</ref>

==== BBP digit-extraction algorithm for {{pi}} ====
We would like to define a formula that returns the (<math>n + 1</math>)-th (with <math>n \geq 0</math>) [[hexadecimal]] digit of {{pi}}.  A few manipulations are required to implement a [[spigot algorithm]] using this formula.

We must first rewrite the formula as:

: <math>\pi = 4 \sum_{k=0}^\infty \frac{1}{\left(16^k\right)(8k+1)} - 2 \sum_{k=0}^\infty \frac{1}{\left(16^k\right)(8k+4)} - \sum_{k=0}^\infty \frac{1}{\left(16^k\right)(8k+5)} - \sum_{k=0}^\infty \frac{1}{\left(16^k\right)(8k+6)}.</math>

Now, for a particular value of ''n'' and taking the first sum, we split the [[series (mathematics)|sum]] to [[infinity]] across the ''n''th term:

: <math>\sum_{k=0}^\infty \frac{1}{\left(16^k\right)(8k+1)} = \sum_{k=0}^n \frac{1}{\left(16^k\right)(8k+1)} + \sum_{k=n+1}^\infty \frac{1}{\left(16^k\right)(8k+1)}.</math>

We now multiply by 16<sup>''n''</sup>, so that the hexadecimal point (the divide between fractional and integer parts of the number) shifts (or remains, if ''n = 0'') to the left of the ''(n+1)''-th fractional digit:

: <math>\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}.</math>

Since we only care about the fractional part of the sum, we look at our two terms and realise that only the first sum contains terms with an integer part; conversely, the second sum doesn't contain terms with an integer part, since the numerator can never be larger than the denominator for ''k''&nbsp;>&nbsp;''n''. Therefore, we need a trick to remove the integer parts, that we don't need, from the terms of the first sum, in order to speed up and increase the precision of the calculations. That trick is to reduce modulo &nbsp;8''k''&nbsp;+&nbsp;1. Our first sum (out of four) to compute the fractional part then becomes:
: <math>\sum_{k=0}^n \frac{16^{n-k} \bmod (8k+1)}{8k+1} + \sum_{k=n+1}^\infty \frac{16^{n-k}}{8k+1}.</math>

Notice how the [[Modular arithmetic|modulus]] operator always guarantees that only the fractional parts of the terms of the first sum will be kept.  To calculate 16<sup>''n''−''k''</sup>&nbsp;mod&nbsp;(8''k''&nbsp;+&nbsp;1) quickly and efficiently, the [[modular exponentiation]] algorithm is done at the same loop level, not [[Nested loop|nested]]. When its running 16''x'' product becomes greater than one, the modulus is taken, just as for the running total in each sum.

Now to complete the calculation, this must be applied to each of the four sums in turn. Once this is done, the four summations are put back into the sum to {{pi}}:
: <math>4 \Sigma_1 - 2 \Sigma_2 - \Sigma_3 - \Sigma_4.</math>

Since only the fractional part is accurate, extracting the wanted digit requires that one removes the integer part of the final sum, multiplies it by 16 and keeps the integer part to "skim off" the hexadecimal digit at the desired 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 [[carry (arithmetic)|carries]] that are not counted, computers usually perform arithmetic for many bits (32 or 64) and round, and we are only interested in the most significant digit(s). There is a 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.

== BBP compared to other methods of computing {{pi}} ==
This algorithm computes {{pi}} without requiring custom data types having thousands or even millions of digits. The method calculates the ''n''th digit ''without'' calculating the first ''n''&nbsp;−&nbsp;1 digits and can use small, efficient data types. [[Fabrice Bellard]] found a variant of BBP, [[Bellard's formula]], which is faster.

Though the BBP formula can directly calculate the value of any given digit of {{pi}} with less computational effort than formulas that must calculate all intervening digits, BBP remains [[Time complexity#Linearithmic_time|linearithmic]] (<math>O(n \log n)</math>), whereby successively larger values of ''n'' require increasingly more time to calculate; that is, the "further out" a digit is, the longer it takes BBP to calculate it, just like the standard {{pi}}-computing algorithms.<ref>{{cite web
 | last   = Bailey
 | first  = David H.
 | title  = The BBP Algorithm for Pi
 | date   = 8 September 2006
 | quote  = Run times for the BBP algorithm ... increase roughly linearly with the position ''d''.
 | url    = http://www.experimentalmath.info/bbp-codes/bbp-alg.pdf
 | accessdate = 17 January 2013}}</ref>

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

== See also ==
* [[Approximations of π|Approximations of {{pi}}]]
* [[Experimental mathematics]]
* [[Bellard's formula]]

== References ==
{{reflist}}

== Further reading ==
* D. J. Broadhurst, [https://arxiv.org/abs/math.CA/9803067 "Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5)"], (1998) ''[[arXiv]]'' math.CA/9803067

== External links ==
* [[Richard J. Lipton]], "[http://rjlipton.wordpress.com/2010/07/14/making-an-algorithm-an-algorithm-bbp/ Making An Algorithm An Algorithm — BBP]", weblog post, July 14, 2010.
* [[Richard J. Lipton]], "[http://rjlipton.wordpress.com/2009/03/15/cooks-class-contains-pi/ Cook’s Class Contains Pi]", weblog post, March 15, 2009.
* {{cite web
|url=https://www.davidhbailey.com/dhbpapers/bbp-formulas.pdf
|first1=David H.
|last1=Bailey
|title=A compendium of BBP-type formulas for mathematical constants, updated 15 Aug 2017
|accessdate=2019-03-31}}
* [[David H. Bailey (mathematician)|David H. Bailey]], "[http://www.experimentalmath.info/bbp-codes/ BBP Code Directory]", web page with links to Bailey's code implementing the BBP algorithm, September 8, 2006.

{{DEFAULTSORT:Bailey-Borwein-Plouffe formula}}
[[Category:Pi algorithms]]
[[Category:Experimental mathematics]]