forked from TheAlgorithms/Python
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathchudnovsky_algorithm.py
More file actions
60 lines (48 loc) · 2 KB
/
Copy pathchudnovsky_algorithm.py
File metadata and controls
60 lines (48 loc) · 2 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
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()
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)}")