Оптимизация итерационного пи-алгоритма с порядком сходимости 2m+1 ⇐ Python

[Anonymous]

Я вывел эту рекуррентную формулу для π с порядком сходимости 2m+1:


Она контролирует, сколько дополнительных правильных цифр π вы получаете на каждом шаге, позволяя добавлять огромные порции цифр за итерацию, а не несколько за раз. время. Например,



Метод
Множитель цифр на итерацию




Фиксированная точка/простая итерация
× 1


Метод Ньютона
× 2


Метод Галлея
× 3


Метод порядка домохозяина (p)
× p


Брента–Саламина / Гаусса–Лежандра (AGM)
× 2


Алгоритм квадратичного π Борвейна
× 2


Квартичный π-алгоритм Борвейна
× 4


Нониковый π-алгоритм Борвейна
× 9


Эта формула (параметр) m)
× 2m+1



Моя цель — реализовать и оптимизировать формулу на Python. Вот что я пробовал на данный момент:

Код: Выделить всё

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")

На моем компьютере выводятся правильные цифры числа «пи».

Код: Выделить всё

Matching digits: 100001
Time Taken: 0.093684 seconds

Однако я не думаю, что он еще полностью оптимизирован — математически, алгоритмически или синтаксически.
Есть ли какие-либо советы по оптимизации для повышения скорости итерационного вычисления? На данный момент я хотел бы вывести 100001 цифру за более короткий период времени, если это возможно.

Подробнее здесь: https://stackoverflow.com/questions/798 ... -order-2m1