Optimierung eines iterativen Pi-Algorithmus mit der Konvergenzordnung 2m+1 ⇐ Python

[Anonymous]

Ich habe diese Wiederholungsformel für π mit der Konvergenzordnung 2m+1 abgeleitet:


Sie steuert, wie viele zusätzliche korrekte π-Ziffern Sie bei jedem Schritt erhalten, sodass Sie pro Iteration große Ziffernblöcke hinzufügen können, anstatt nur einige auf einmal. Zum Beispiel:



Methode
Ziffernmultiplikator pro Iteration




Festkomma-/einfache Iteration
× 1


Newtons Methode
× 2


Halleys Methode
× 3


Haushaltsmethode der Ordnung (p)
× p


Brent–Salamin / Gauss–Legendre (AGM)
× 2


Borwein quadratischer π-Algorithmus
× 2


Quartischer π-Algorithmus von Borwein
× 4


nonischer π-Algorithmus von Borwein
× 9


Diese Formel (Parameter m)
× 2m+1



Mein Ziel ist es, die Formel in Python zu implementieren und zu optimieren. Folgendes habe ich bisher versucht:

Code: Select all

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[i] == s_pi[i]:
i += 1

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

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

Auf meinem PC werden diese korrekten Pi-Ziffern ausgegeben

Code: Select all

Matching digits: 100001
Time Taken: 0.093684 seconds

Ich glaube jedoch nicht, dass es noch vollständig optimiert ist – mathematisch, algorithmisch oder syntaktisch.
Gibt es Optimierungstipps zur Verbesserung der Iterations-Rechengeschwindigkeit? Im Moment möchte ich möglichst in kürzerer Zeit 100001 Ziffern ausgeben.