# Reviewer's independent check at p=24 (binary32) and p=11 (binary16) with numpy native types,
# plus exact first escape at p=24 by outward-rounded dyadic intervals (256 bits).
import numpy as np
from fractions import Fraction
def run(dt, p):
    c = dt(0.25) + dt(2.0**-(p+1)); assert float(c) == 0.25 + 2.0**-(p+1)
    out = {}
    for mode in ("separate", "fused"):
        x = dt(0); n = 0
        while True:
            n += 1
            if mode == "separate":
                nx = dt(dt(x*x) + c)
            else:
                nx = dt(float(x)*float(x) + float(c))  # exact in binary64 here, single rounding
            assert 0 <= nx <= 0.5
            if nx == x: out[mode] = (n, float(nx)); break
            x = nx
    return out
print("binary32", run(np.float32, 24))
print("binary16", run(np.float16, 11))
def esc(p, bits=256):
    D = 1 << bits; c = Fraction(1,4) + Fraction(1, 2**(p+1))
    clo = c.numerator*D // c.denominator; chi = -(-c.numerator*D // c.denominator)
    lo = hi = 0; n = 0
    while True:
        n += 1; lo = lo*lo//D + clo; hi = -(-hi*hi//D) + chi
        if lo > 2*D: return n
        assert hi <= 2*D
print("exact first escape p=24:", esc(24), " p=11:", esc(11))