bailey borwein plouffe Algorithm
The Bailey-Borwein-Plouffe (BBP) algorithm is a remarkable mathematical formula and algorithm that allows for the computation of individual hexadecimal or binary digits of mathematical constants, such as π (pi), log(2), and other irrational numbers, without having to compute any of the preceding digits. It was discovered by David H. Bailey, Peter Borwein, and Simon Plouffe in 1995, and it represented a significant breakthrough in the field of computational mathematics. The BBP algorithm relies on a spigot algorithm, which is a type of algorithm that produces a stream of digits as output, rather than calculating all the digits at once. This method is especially useful for computing digits of irrational numbers that have an infinite decimal expansion, as it allows for the extraction of specific digits without the need to calculate all the preceding ones.
At the core of the BBP algorithm is a formula that expresses constants like π as an infinite series involving simple mathematical operations and rational numbers. The algorithm uses this formula to calculate the nth digit of the constant by summing up a small number of terms in the series, where n is the desired digit position. By exploiting the properties of modular arithmetic and the radix representation of numbers, the algorithm can efficiently compute individual digits with minimal computational effort. This method has been used to calculate billions of digits of π and other mathematical constants, and it has furthered our understanding of their properties and distribution of digits. Moreover, the discovery of the BBP algorithm has inspired the development of similar digit-extraction algorithms for other constants, expanding the field of computational mathematics and opening up new possibilities for research.
def bailey_borwein_plouffe(digit_position: int, precision: int = 1000) -> str:
"""
Implement a popular pi-digit-extraction algorithm known as the
Bailey-Borwein-Plouffe (BBP) formula to calculate the nth hex digit of pi.
Wikipedia page:
https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
@param digit_position: a positive integer representing the position of the digit to extract.
The digit immediately after the decimal point is located at position 1.
@param precision: number of terms in the second summation to calculate.
A higher number reduces the chance of an error but increases the runtime.
@return: a hexadecimal digit representing the digit at the nth position
in pi's decimal expansion.
>>> "".join(bailey_borwein_plouffe(i) for i in range(1, 11))
'243f6a8885'
>>> bailey_borwein_plouffe(5, 10000)
'6'
>>> bailey_borwein_plouffe(-10)
Traceback (most recent call last):
...
ValueError: Digit position must be a positive integer
>>> bailey_borwein_plouffe(0)
Traceback (most recent call last):
...
ValueError: Digit position must be a positive integer
>>> bailey_borwein_plouffe(1.7)
Traceback (most recent call last):
...
ValueError: Digit position must be a positive integer
>>> bailey_borwein_plouffe(2, -10)
Traceback (most recent call last):
...
ValueError: Precision must be a nonnegative integer
>>> bailey_borwein_plouffe(2, 1.6)
Traceback (most recent call last):
...
ValueError: Precision must be a nonnegative integer
"""
if (not isinstance(digit_position, int)) or (digit_position <= 0):
raise ValueError("Digit position must be a positive integer")
elif (not isinstance(precision, int)) or (precision < 0):
raise ValueError("Precision must be a nonnegative integer")
# compute an approximation of (16 ** (n - 1)) * pi whose fractional part is mostly accurate
sum_result = (
4 * _subsum(digit_position, 1, precision)
- 2 * _subsum(digit_position, 4, precision)
- _subsum(digit_position, 5, precision)
- _subsum(digit_position, 6, precision)
)
# return the first hex digit of the fractional part of the result
return hex(int((sum_result % 1) * 16))[2:]
def _subsum(
digit_pos_to_extract: int, denominator_addend: int, precision: int
) -> float:
# only care about first digit of fractional part; don't need decimal
"""
Private helper function to implement the summation
functionality.
@param digit_pos_to_extract: digit position to extract
@param denominator_addend: added to denominator of fractions in the formula
@param precision: same as precision in main function
@return: floating-point number whose integer part is not important
"""
sum = 0.0
for sum_index in range(digit_pos_to_extract + precision):
denominator = 8 * sum_index + denominator_addend
exponential_term = 0.0
if sum_index < digit_pos_to_extract:
# if the exponential term is an integer and we mod it by the denominator before
# dividing, only the integer part of the sum will change; the fractional part will not
exponential_term = pow(
16, digit_pos_to_extract - 1 - sum_index, denominator
)
else:
exponential_term = pow(16, digit_pos_to_extract - 1 - sum_index)
sum += exponential_term / denominator
return sum
if __name__ == "__main__":
import doctest
doctest.testmod()
