Spigot algorithm

A spigot algorithm is a type of algorithm used to compute the value of a mathematical constant such as $π$ or e. Spigot algorithms are unique because they do not require the total number of digits to be fixed beforehand, and do not require the computation of several intermediate results which are combined to produce the final result.^[1] There are two kinds of spigot algorithms: (1) those that can produce a single, arbitrary digit (also called digit extraction algorithm); and (2) those that produce a sequence of digits, one after the other. The Bailey-Borwein-Plouffe formula is a digit extraction algorithm for $π$ which produces hexadecimal digits. A sequential spigot algorithm for π was produced by Stanley Rabinowitz and Stanley Wagon (this algorithm is sometimes referred to as "the spigot algorithm for $π$ ").^[2] Spigot algorithms that produce a sequence of digits begin producing digits and produce them continuously, rather than waiting for the entire algorithm to finish. Spigot algorithms are known for $π$ and e. Spigot algorithms typically work in a particular radix, such as hexadecimal or binary number.

Example

This example illustrates the working of a spigot algorithm by calculating the binary digits of the natural logarithm of 2 (sequence A068426 in the OEIS) using the identity

\ln(2)=\sum _{k=1}^{\infty }{\frac {1}{k2^{k}}}\,.

To start calculating binary digits from, say, the 8th place we multiply this identity by 2⁷(since 7 = 8 - 1):

2^{7}\ln(2)=2^{7}\sum _{k=1}^{\infty }{\frac {1}{k2^{k}}}\,.

We then divide the infinite sum into a "head", in which the exponents of 2 are greater than or equal to zero, and a "tail", in which the exponents of 2 are negative:

2^{7}\ln(2)=\sum _{k=1}^{7}{\frac {2^{7-k}}{k}}+\sum _{k=8}^{\infty }{\frac {1}{k2^{k-7}}}\,.

We are only interested in the fractional part of this value, so we can replace each of the summands in the "head" by

{\frac {2^{7-k}\mod k}{k}}\,.

Calculating each of these terms and adding them to a running total where we again only keep the fractional part, we have:

k	A = 2^7-k	B = A mod k	C = B / k	Sum of C mod 1
1	64	0	0	0
2	32	0	0	0
3	16	1	1/3	1/3
4	8	0	0	1/3
5	4	4	4/5	2/15
6	2	2	1/3	7/15
7	1	1	1/7	64/105

We add a few terms in the "tail", noting that the error introduced by truncating the sum is less than the final term:

k	D = 1/k2^k-7	Sum of D	Maximum error
8	1/16	1/16	1/16
9	1/36	13/144	1/36
10	1/80	37/360	1/80

Adding the "head" and the first few terms of the "tail" together we get:

2^{7}\ln(2)\mod {1}\approx {\frac {64}{105}}+{\frac {37}{360}}=0.10011100\cdots _{2}+0.00011010\cdots _{2}=0.1011\cdots _{2}\,,

so the 8th to 11th binary digits in the binary expansion of ln(2) are 1, 0, 1, 1. Note that we have not calculated the values of the first seven binary digits - indeed, all information about them has been intentionally discarded by using modular arithmetic in the "head" sum.

The same approach can be used to calculate digits of the binary expansion of ln(2) starting from an arbitrary n^th position. The number of terms in the "head" sum increases linearly with n, but the complexity of each term only increases with the logarithm of n if an efficient method of modular exponentiation is used. The precision of calculations and intermediate results and the number of terms taken from the "tail" sum are all independent of n, and only depend on the number of binary digits that are being calculated - single precision arithmetic can be used to calculate around 12 binary digits, regardless of the starting position.

Footnotes

^ Arndt, pp 37, 77–85;
^ Arndt, pp 37, 77–85;

References

Arndt, Jorg; Haenel, Christoph, $π$ unleashed, Springer Verlag, 2000.

External links

[1] Arndt, pp 37, 77–85;

[2] Arndt, pp 37, 77–85;

[1]

[2]