pie.py

from mpmath import mp
import time

mp.dps = 2000  # working precision; adjust as you like

def compute_coeffs(m):
    """Compute c_{m,k} = (-1)^{m+k} * (1/k) * prod_{j!=k} j^2 / (k^2 - j^2)."""
    coeffs = []
    for k in range(1, m+1):
        prod = mp.mpf('1')
        for j in range(1, m+1):
            if j == k:
                continue
            prod *= (j*j) / (k*k - j*j)
        c = (-1)**(m+k) * (mp.mpf(1)/k) * prod
        coeffs.append(c)
    return coeffs

def step_direct(x, coeffs):
    """One iteration using direct sin(k*x) calls."""
    dx = mp.mpf('0')
    m = len(coeffs)
    for k in range(1, m+1):
        dx += coeffs[k-1] * mp.sin(k * x)
    return x + dx

def step_recur(x, coeffs):
    """One iteration using recurrence for sin(k*x)."""
    dx = mp.mpf('0')
    m = len(coeffs)

    # compute sin(x), cos(x) once
    s1 = mp.sin(x)      # sin(1*x)
    c  = mp.cos(x)      # cos(x)

    # k = 1 term
    s_prev = mp.mpf('0')   # sin(0*x)
    s_curr = s1            # sin(1*x)
    dx += coeffs[0] * s_curr

    # k = 2..m terms via recurrence: sin(kx) = 2 cos x sin((k-1)x) - sin((k-2)x)
    for k in range(2, m+1):
        s_next = 2 * c * s_curr - s_prev
        dx += coeffs[k-1] * s_next
        s_prev, s_curr = s_curr, s_next

    return x + dx

# --------- experiment parameters ---------

m = 30        # try 2, 3, 5, 10, 20, ...
N = 1000      # number of iterations in the timing loop

coeffs = compute_coeffs(m)

# sanity check: direct and recur match for a single step
x0 = mp.mpf(3)
y_direct = step_direct(x0, coeffs)
y_recur  = step_recur(x0, coeffs)
print("Sanity check diff:", y_direct - y_recur)

# --- time direct version ---
x = mp.mpf(3)
t0 = time.perf_counter()
for _ in range(N):
    x = step_direct(x, coeffs)
t1 = time.perf_counter()
print(f"direct time (m={m}, N={N}): {t1 - t0:.6f} s")

# --- time recurrence version ---
x = mp.mpf(3)
t2 = time.perf_counter()
for _ in range(N):
    x = step_recur(x, coeffs)
t3 = time.perf_counter()
print(f"recur  time (m={m}, N={N}): {t3 - t2:.6f} s")

# show final approximation (optional)
print("Final x (recur):", mp.nstr(x, 50))