Controlling pi digit generation

from mpmath import mp

def build_coeffs(m, prec):
    """
    Build coefficients c_r = r! / (2r+1)!! using the recurrence
        c_0 = 1
        c_r = c_{r-1} * r / (2r+1)
    at precision 'prec'.
    """
    old_dps = mp.dps
    mp.dps = prec

    coeffs = [mp.mpf('1')]
    for r in range(1, m):
        coeffs.append(coeffs[-1] * r / (2*r + 1))

    mp.dps = old_dps
    return coeffs

def step(x, coeffs):
    """
    One iteration of the recurrence:
        x_{n+1} = x_n + sin(x_n) * sum_{r=0}^{m-1} c_r * (cos x_n + 1)^r
    evaluated via Horner's rule in t = cos(x) + 1.
    """
    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

def matching_digits(x, dps=None):
    """
    Count how many decimal digits of x match pi (ignoring the '3.' prefix).
    """
    if dps is None:
        dps = mp.dps

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

    # Compare characters from the start
    L = min(len(s_x), len(s_pi))
    i = 0
    while i < L and s_x[i] == s_pi[i]:
        i += 1

    # Subtract 2 to ignore '3.' (the integer part and decimal point)
    return max(0, i - 2)

def run_pi_from_x0(m, x0_str, num_steps, safety=30):
    """
    Run the pi recurrence starting from x0 (as a string) and parameter m.
    The convergence order is 2m+1.

    We use the error model
        N_{n+1} ≈ (2m+1) N_n + log10((2m+1)! / (m!)^2)
    to:
      - estimate how many digits N_n we get after each step,
      - choose mp.dps for each step,
      - choose how many digits to use when building the coefficients.
    """
    # --- Initial estimate for digits in x0 ---
    mp.dps = 80
    x = mp.mpf(x0_str)
    N0 = -mp.log10(abs(x - mp.pi))      # N_0 ≈ -log10|x0 - pi|

    print(f"Initial x0 = {x0_str}")
    #print(f"Estimated N0 ≈ {float(N0):.3f} correct digits")

    # --- Digit-growth constants ---
    a = 2*m + 1
    B = mp.log10(mp.factorial(2*m + 1)) - 2*mp.log10(mp.factorial(m))
    #print(f"Convergence order = {a},  B ≈ {float(B):.6f}")
    print(f"Convergence order = {a}")

    # --- Predict max digits over all steps (to size coeff precision) ---
    N_tmp = N0
    max_N = N0
    for _ in range(num_steps):
        N_tmp = a * N_tmp + B
        if N_tmp > max_N:
            max_N = N_tmp

    coeff_prec = int(max_N + safety)
    if coeff_prec < 80:
        coeff_prec = 80

    print(f"Building coefficients at precision = {coeff_prec} dps")
    coeffs = build_coeffs(m, coeff_prec)

    # --- Main iteration loop ---
    N = N0
    for step_index in range(1, num_steps + 1):
        # Predicted digits after this step
        N = a * N + B

        # Working precision for this step
        dps_for_step = int(mp.ceil(N + safety))
        if dps_for_step < 50:
            dps_for_step = 50
        if dps_for_step > coeff_prec:
            dps_for_step = coeff_prec

        mp.dps = dps_for_step
        x = step(x, coeffs)

        # Actual digits that match pi at this precision
        actual = matching_digits(x, dps_for_step)

        print(
            f"Step {step_index}: "
            f"dps={dps_for_step}, "
            f"predicted ~{int(N)} digits, "
            f"actual {actual} matching digits"
        )

    return x

# Example usage:
m = 31          # convergence order = 2m+1 = 63
x0 = "3.0"      # starting approximation
num_steps = 4   # how many iterations to run
pi_approx = run_pi_from_x0(m, x0, num_steps)