|
| 1 | +from decimal import Decimal, getcontext |
| 2 | +from math import ceil, factorial |
| 3 | + |
| 4 | + |
| 5 | +def pi(precision: int) -> str: |
| 6 | + """ |
| 7 | + The Chudnovsky algorithm is a fast method for calculating the digits of PI, |
| 8 | + based on Ramanujan’s PI formulae. |
| 9 | +
|
| 10 | + https://en.wikipedia.org/wiki/Chudnovsky_algorithm |
| 11 | +
|
| 12 | + PI = constant_term / ((multinomial_term * linear_term) / exponential_term) |
| 13 | + where constant_term = 426880 * sqrt(10005) |
| 14 | +
|
| 15 | + The linear_term and the exponential_term can be defined iteratively as follows: |
| 16 | + L_k+1 = L_k + 545140134 where L_0 = 13591409 |
| 17 | + X_k+1 = X_k * -262537412640768000 where X_0 = 1 |
| 18 | +
|
| 19 | + The multinomial_term is defined as follows: |
| 20 | + 6k! / ((3k)! * (k!) ^ 3) |
| 21 | + where k is the k_th iteration. |
| 22 | +
|
| 23 | + This algorithm correctly calculates around 14 digits of PI per iteration |
| 24 | +
|
| 25 | + >>> pi(10) |
| 26 | + '3.14159265' |
| 27 | + >>> pi(100) |
| 28 | + '3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706' |
| 29 | + >>> pi('hello') |
| 30 | + Traceback (most recent call last): |
| 31 | + ... |
| 32 | + TypeError: Undefined for non-integers |
| 33 | + >>> pi(-1) |
| 34 | + Traceback (most recent call last): |
| 35 | + ... |
| 36 | + ValueError: Undefined for non-natural numbers |
| 37 | + """ |
| 38 | + |
| 39 | + if not isinstance(precision, int): |
| 40 | + raise TypeError("Undefined for non-integers") |
| 41 | + elif precision < 1: |
| 42 | + raise ValueError("Undefined for non-natural numbers") |
| 43 | + |
| 44 | + getcontext().prec = precision |
| 45 | + num_iterations = ceil(precision / 14) |
| 46 | + constant_term = 426880 * Decimal(10005).sqrt() |
| 47 | + multinomial_term = 1 |
| 48 | + exponential_term = 1 |
| 49 | + linear_term = 13591409 |
| 50 | + partial_sum = Decimal(linear_term) |
| 51 | + for k in range(1, num_iterations): |
| 52 | + multinomial_term = factorial(6 * k) // (factorial(3 * k) * factorial(k) ** 3) |
| 53 | + linear_term += 545140134 |
| 54 | + exponential_term *= -262537412640768000 |
| 55 | + partial_sum += Decimal(multinomial_term * linear_term) / exponential_term |
| 56 | + return str(constant_term / partial_sum)[:-1] |
| 57 | + |
| 58 | + |
| 59 | +if __name__ == "__main__": |
| 60 | + n = 50 |
| 61 | + print(f"The first {n} digits of pi is: {pi(n)}") |
0 commit comments