coderjack
diff --git a/‎maths/chudnovsky_algorithm.py‎
Lines changed: 61 additions & 0 deletions b/‎maths/chudnovsky_algorithm.py‎
Lines changed: 61 additions & 0 deletions
@@ -0,0 +1,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)}")