[Python] Optimizing an iterative pi algorithm with convergence order 2m+1

I've derived this recurrence formula for π with convergence order 2m+1:



It controls how many extra correct π digits you gain each step, letting you add huge chunks of digits per iteration instead of just a few at a time. For example,

Method Digits multiplier per iteration
Fixed-point / simple iteration × 1
Newton’s method × 2
Halley’s method × 3
Householder method of order (p) × p
Brent–Salamin / Gauss–Legendre (AGM) × 2
Borwein quadratic π algorithm × 2
Borwein quartic π algorithm × 4
Borwein nonic π algorithm × 9
This formula (parameter m) × 2m+1

My goal is to implement and optimize the formula in Python. Here's what I've tried so far:

from mpmath import mp
import time

# High precision
mp.dps = 316000

m = 31 # method has convergence order 2*m + 1 (= 63)


def make_coeffs(m: int):
"""
Precompute c_r = r! / (2r+1)!! for r = 0..m-1 via the recurrence
c_0 = 1
c_r = c_{r-1} * r / (2r+1)
so we never call factorial or double-factorial at runtime.
"""
coeffs = [mp.mpf('0')] * m
c = mp.mpf('1') # c_0 = 0! / 1!! = 1
coeffs[0] = c
for r in range(1, m):
c = c * r / (2 * r + 1)
coeffs[r] = c
return coeffs


COEFFS = make_coeffs(m)


def step(x, coeffs=COEFFS):
"""
One iteration of:
x_{n+1} = x_n + sin(x_n) * sum_{r=0}^{m-1} c_r * (cos x_n + 1)^r
where c_r = r! / (2r+1)!!, evaluated by Horner's rule.
"""
c, s = mp.cos_sin(x)
t = c + 1
total = coeffs[-1]
for a in reversed(coeffs[:-1]):
total = total * t + a
return x + s * total


x = mp.mpf('3.0')
start = time.perf_counter()
with mp.workdps(10): x = step(x)
with mp.workdps(100): x = step(x)
with mp.workdps(5000): x = step(x)
with mp.workdps(100000): x = step(x)
end = time.perf_counter()

s_x = mp.nstr(x, 100000)
s_pi = mp.nstr(mp.pi, 100000)

i = 0
while i < len(s_x) and s_x == s_pi:
i += 1

#print("x =", s_x)
#print("pi =", s_pi)

print("Matching digits:", i)
print(f"Time Taken: {end - start:.6f} seconds")


On my PC, it outputs these correct pi digits

Matching digits: 100001
Time Taken: 0.093684 seconds


However, I don't think it's fully optimized yet - mathematically, algorithmically, or syntactically.

Are there any optimization tips to improve the iteration computational speed? For now, I'd like to output 100001 digits in a shorter period of time, if possible.

Continue reading...