Bailey–Borwein–Plouffe formula

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:

${\displaystyle {\begin{aligned}\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{aligned}}}$

${\displaystyle {\begin{aligned}\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{aligned}}}$

(In fact, this identity holds true for a > 1:

${\displaystyle \ln {\frac {a}{a-1}}=\sum _{k=1}^{\infty }{\frac {1}{a^{k}\cdot k}}}$.)

Plouffe was also inspired by the arctan power series of the form (the P notation can be also generalized to the case where b is not an integer):

${\displaystyle {\begin{aligned}\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{aligned}}}$

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 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-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 original BBP π summation formula was found in 1995 by Plouffe using PSLQ. It is also representable using the P function:

${\displaystyle {\begin{aligned}\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{aligned}}}$

which also reduces to this equivalent ratio of two polynomials:

${\displaystyle \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].}$

This formula has been shown through a fairly simple proof to equal π.[7]

We would like to define a formula that returns the nth hexadecimal digit of π. A few manipulations are required to implement a spigot algorithm using this formula.

We must first rewrite the formula as

${\displaystyle \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)}}.}$

Now, for a particular value of n and taking the first sum, we split the sum to infinity across the nth term:

${\displaystyle \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)}}.}$

We now multiply by 16ⁿ, so that the hexadecimal point (the divide between fractional and integer parts of the number) is in the nth place:

${\displaystyle \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; conversely, 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 to reduce modulo  8k + 1. Our sum for the first fractional part then becomes

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

Notice how the modulus operator always guarantees that only the fractional sum will be kept. To calculate 16^n−k mod (8k + 1) quickly and efficiently, the modular exponentiation algorithm is used. When the running 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 π:

${\displaystyle 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 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.