FilesExpand file tree

master

/

chudnovsky_algorithm.py

Copy path

Blame

More file actions

Blame

More file actions

Latest commit

 

History

History

History

61 lines (49 loc) · 2.03 KB

master

/

chudnovsky_algorithm.py

Copy path

Top

File metadata and controls

• Code
• Blame

61 lines (49 loc) · 2.03 KB

Raw

Copy raw file

Download raw file

Open symbols panel

Edit and raw actions

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

from decimal import Decimal, getcontext

from math import ceil, factorial

def pi(precision: int) -> str:

"""

The Chudnovsky algorithm is a fast method for calculating the digits of PI,

based on Ramanujan’s PI formulae.

https://en.wikipedia.org/wiki/Chudnovsky_algorithm

PI = constant_term / ((multinomial_term * linear_term) / exponential_term)

where constant_term = 426880 * sqrt(10005)

The linear_term and the exponential_term can be defined iteratively as follows:

L_k+1 = L_k + 545140134 where L_0 = 13591409

X_k+1 = X_k * -262537412640768000 where X_0 = 1

The multinomial_term is defined as follows:

6k! / ((3k)! * (k!) ^ 3)

where k is the k_th iteration.

This algorithm correctly calculates around 14 digits of PI per iteration

>>> pi(10)

'3.14159265'

>>> pi(100)

'3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706'

>>> pi('hello')

Traceback (most recent call last):

...

TypeError: Undefined for non-integers

>>> pi(-1)

Traceback (most recent call last):

...

ValueError: Undefined for non-natural numbers

"""

if not isinstance(precision, int):

raise TypeError("Undefined for non-integers")

elif precision < 1:

raise ValueError("Undefined for non-natural numbers")

getcontext().prec = precision

num_iterations = ceil(precision / 14)

constant_term = 426880 * Decimal(10005).sqrt()

multinomial_term = 1

exponential_term = 1

linear_term = 13591409

partial_sum = Decimal(linear_term)

for k in range(1, num_iterations):

multinomial_term = factorial(6 * k) // (factorial(3 * k) * factorial(k) ** 3)

linear_term += 545140134

exponential_term *= -262537412640768000

partial_sum += Decimal(multinomial_term * linear_term) / exponential_term

return str(constant_term / partial_sum)[:-1]

if __name__ == "__main__":

n = 50

print(f"The first {n} digits of pi is: {pi(n)}")