From f0dfc4f46d47102bf34c09a3920138caf9cd9ae5 Mon Sep 17 00:00:00 2001 From: TheSuperNoob Date: Wed, 12 Feb 2020 15:04:59 +0100 Subject: [PATCH] Add Chudnovskys algorithm for calculating many digits of pi (#1752) * Add Chudnovskys algorithm for calculating many digits of pi * Update return value type hint * Initialize partial sum to be of type Decimal * Update chudnovsky_algorithm.py Co-authored-by: Christian Clauss --- maths/chudnovsky_algorithm.py | 61 +++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 maths/chudnovsky_algorithm.py diff --git a/maths/chudnovsky_algorithm.py b/maths/chudnovsky_algorithm.py new file mode 100644 index 000000000..fb188cd6a --- /dev/null +++ b/maths/chudnovsky_algorithm.py @@ -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)}")