"""Independent recomputation of first-escape counts for c_T = 1/4 + 1/(16 (T+1)^2), z0 = 0.
Written by the reviewer from the paper's definitions, without reusing the bundle's code.
Method A: outward-rounded dyadic interval arithmetic at several precisions (integers only).
Method B: mpmath interval arithmetic (mpmath.iv) at 300 bits.
Also: the largest critical iterate over 0..T (checks C1's bound 17/32 numerically), and
Klebanoff's asymptotic N ~ pi / sqrt(eps) for comparison."""
import json, math
from fractions import Fraction
import mpmath
from mpmath import iv

CUTOFFS = [1, 4, 16, 64, 256, 1024, 4096]

def c_of(T):
    return Fraction(1, 4) + Fraction(1, 16 * (T + 1) ** 2)

def dyadic(T, bits):
    c = c_of(T); D = 1 << bits
    clo = (c.numerator * D) // c.denominator
    chi = -((-c.numerator * D) // c.denominator)
    lo = hi = 0; n = 0; max_hi_to_T = 0
    while True:
        n += 1
        lo = (lo * lo) // D + clo
        hi = -((-hi * hi) // D) + chi
        if n <= T:
            max_hi_to_T = max(max_hi_to_T, hi)
        if lo > 2 * D:
            return n, Fraction(max_hi_to_T, D)
        if hi > 2 * D:
            return None, None  # undecided at this precision

def mp_iv(T, bits=300):
    iv.prec = bits
    c = iv.mpf(c_of(T).numerator) / c_of(T).denominator
    z = iv.mpf(0); n = 0
    while True:
        n += 1
        z = z * z + c
        if z.a > 2:
            return n
        if z.b > 2:
            return None

out = []
for T in CUTOFFS:
    a128, m128 = dyadic(T, 128)
    a256, m256 = dyadic(T, 256)
    b = mp_iv(T)
    eps = 1 / (16 * (T + 1) ** 2)
    out.append({"cutoff": T, "dyadic128": a128, "dyadic256": a256, "mpmath_iv300": b,
                "max_upper_iterate_through_T": float(m256), "bound_17_32": 17 / 32,
                "klebanoff_pi_over_sqrt_eps": round(math.pi / math.sqrt(eps), 1)})
    print(out[-1], flush=True)
json.dump(out, open("independent_check.json", "w"), indent=2)
