Sunbelt Computer Software

PL/B Language Development and Support

TheAlgorithmsInPython/maths/bailey_borwein_plouffe.py at master · unconnect/TheAlgorithmsInPython · GitHub

90 lines (82 loc) · 3.54 KB

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

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()