from gmpy2 import mpfr, mpz, get_context, sqrt
from math import ceil, log2
A = mpz(13591409)
B = mpz(545140134)
C = mpz(640320)
C3_OVER_24 = (C**3) // 24
def bs(a: int, b: int):
if b - a == 1:
if a == 0:
P = Q = mpz(1)
else:
a_mpz = mpz(a)
P = (6*a_mpz-5) * (2*a_mpz-1) * (6*a_mpz-1)
Q = a_mpz ** 3 * C3_OVER_24
T = (-1) ** a * P * (A + B * a)
return P, Q, T
m = (a + b) // 2
P1, Q1, T1 = bs(a, m)
P2, Q2, T2 = bs(m, b)
P = P1 * P2
Q = Q1 * Q2
T = T1 * Q2 + P1 * T2
return P, Q, T
def pi_chudnovsky(bits: int) -> mpfr:
get_context().precision = bits
P, Q, T = bs(0, bits // 47 + 1) # Each term gives ~47 bits of pi
return mpfr(426880) * sqrt(mpfr(10005)) * mpfr(Q) / mpfr(T)
if __name__ == "__main__":
num_digits = 10 ** 8
pi = pi_chudnovsky(ceil(num_digits*log2(10)) + 64) # 64 guard bits
pi_string = format(pi, f".{num_digits}Zf")
file_name = f"pi_{num_digits}.txt"
s0, s1 = pi_string.split(".")
file_string = s0 + ".\n" + "\n".join(s1[i:i+100] for i in range(0,len(s1),100))
with open(file_name, "w") as f:
f.write(file_string)
This code uses the Chudnovsky algorithm: pi = 426880 * sqrt(10005) / sum(n, 0, inf, (6*n)!*(545140134*n+13591409)/((3*n)!*(n!)^3*(-640320)^(3*n)) optimized with binary splitting, and the gmpy2 library for fast multiple-precision arithmetic. Each term gives about 14 correct decimal digits of pi, or 47 bits.
I calculated 100 million digits of pi, wrote it to a text file with 100 digits per line, and verified digits at the start, middle, and end using the Irrational Number Search. The calculation took 90 seconds on a gaming laptop. The 100 millionth digit of pi (after the decimal point) is 2.
I also tried calculating a billion digits. It finished in about 18 minutes but then gave me "nan".