"""
series_newton.py -- FAST, STABLE Laurent coefficients b_m of the exterior map of the
Mandelbrot set, and the Gronwall area-theorem upper bound A_n = pi(1 - sum_{m<=n} m b_m^2).

WHY A NEW METHOD (the open edge series.py/series_fast.py named).
  The forward map Phi_M(c) = c + sum a_j c^{-j} has Laurent coefficients that GROW like ~2^j
  (M reaches c=-2, so the series in 1/c has radius 1/2).  Any method that forms Phi's
  coefficients and reverses the series must recover the O(1) inverse coefficients b_m from
  ~2^N-sized numbers -- it needs ~0.31 N digits of precision, so float64 dies by N~=50
  (series_fast.py, the "recorded dead end") and mpmath dps=60 stalls near N~=160 (series.py).

  This file sidesteps the wall.  Write t=1/w, u=t, x=1/c.  With z_k = w^{2^{k-1}} Z_k the
  Boettcher product  w = c prod_{k>=1}(1 + c/z_k^2)^{2^-k}  becomes a FIXED POINT for
  V0 := c/w = c*u, whose coefficients ARE the b_m (V0 = 1 + b_0 u + b_1 u^2 + ...):

        1/V0 = prod_{k>=1} (1 + V0 u^{2^k-1} / Z_k^2)^{2^-k},
        Z_1 = V0,   Z_{k+1} = Z_k^2 + V0 u^{2^k-1}.

  Every Z_k has Z_k(0)=1, so 1/Z_k^2 is a BOUNDED series -- no 2^N growth anywhere.  We solve
  E(V0) := log V0 + sum_k 2^-k log(1 + G_k) = 0  (G_k = V0 u^{2^k-1}/Z_k^2) by Newton doubling.
  The Frechet derivative is an honest multiplication series (all V0-dependence is multiplicative
  or via the linear recurrence M_k = dZ_k/dV0, M_1=1, M_{k+1}=2 Z_k M_k + u^{2^k-1}), so Newton
  converges quadratically: O(log N) iterations, each O(N log N) with FFT.  b_m = [u^{m+1}] V0.

VALIDATION (the __main__ below; verify.mjs drives the page in a browser and does not do this):
  * b_m match the published exact rationals (b_0=-1/2, b_1=1/8, b_2=-1/4, b_3=15/128,
    b_5=-47/1024, b_6=-1/16, b_7=987/32768, b_9=-3673/262144, b_10=1/32) to ~1e-13.
  * A_n reproduces every published checkpoint to 5-decimal agreement, from scratch: Huddleston's
    1990 table 72..115232 (the doubling entries plus 16384 and 115232 are checked below) and
    Ewing-Schober 1992's 240000 -> 1.7274 (as quoted by Bittner et al.; checked when N>=240000).

FLOAT64 REACH.  The Jacobian's series M_k grows ~2^{1.5k}~N^1.5, so float64 walls near
N=2^18=262144 (there the bound is ~1.7262).  Higher N needs extended precision; the field's
frontier value 1.68288 at 5,000,000 terms (Bittner-Cheong-Gates-Nguyen 2017) uses the same
class of stable recursion carried further.
"""
import numpy as np, sys, math, json, time

# ---- FFT power-series ops over R[[u]] truncated to length L (orders 0..L-1) ----
def mul(a,b,L):
    n=len(a)+len(b)-1
    if n<=64: return np.convolve(a,b)[:L]
    size=1<<((n-1).bit_length())
    return np.fft.irfft(np.fft.rfft(a,size)*np.fft.rfft(b,size),size)[:L]
def inv(a,L):
    b=np.array([1.0/a[0]]); m=1
    while m<L:
        m2=min(2*m,L); at=a[:m2] if len(a)>=m2 else np.concatenate([a,np.zeros(m2-len(a))])
        t=-mul(at,b,m2); t[0]+=2.0; b=mul(b,t,m2); m=m2
    return b[:L]
def deriv(a):
    L=len(a); return np.zeros(1) if L<=1 else a[1:]*np.arange(1,L)
def integ(a):
    out=np.zeros(len(a)+1); out[1:]=a/np.arange(1,len(a)+1); return out
def slog(a,L):   # a[0]=1
    return integ(mul(deriv(a),inv(a,L),L))[:L]
def sexp(a,L):   # a[0]=0
    b=np.array([1.0]); m=1
    while m<L:
        m2=min(2*m,L); lb=slog(b[:m2] if len(b)>=m2 else b,m2); t=a[:m2]-lb; t[0]+=1.0
        b=mul(b,t,m2); m=m2
    return b[:L]
def shift(a,s,L):
    out=np.zeros(L)
    if s<L: out[s:]=a[:L-s]
    return out

def _E_and_J(V,L):
    Vinv=inv(V,L); S=np.zeros(L); dS=np.zeros(L)
    Z=V.copy(); M=np.zeros(L); M[0]=1.0
    k=1
    while (1<<k)-1<L:
        s=(1<<k)-1; Zinv=inv(Z,L); Zinv2=mul(Zinv,Zinv,L)
        Gk=shift(mul(V,Zinv2,L),s,L); onePG=Gk.copy(); onePG[0]+=1.0
        opinv=inv(onePG,L); wk=1.0/(1<<k)
        S=S+wk*slog(onePG,L)
        fac=Vinv-2.0*mul(M,Zinv,L)
        dS=dS+wk*mul(mul(Gk,fac,L),opinv,L)
        us=np.zeros(L); us[s]=1.0
        Mnew=2.0*mul(Z,M,L)+us; Znew=mul(Z,Z,L)+shift(V,s,L); Z=Znew; M=Mnew; k+=1
    return slog(V,L)+S, Vinv+dS

def compute_b(N):
    """Return numpy array b with b[m]=b_m for m=0..N (via V0 length N+2)."""
    Ltot=N+2; V=np.array([1.0]); L=1
    while L<Ltot:
        L=min(2*L,Ltot); Vp=np.zeros(L); Vp[:len(V)]=V; V=Vp
        for _ in range(2):
            E,J=_E_and_J(V,L); V=V-mul(E,inv(J,L),L)
    return V[1:N+2]   # b[m]=V[m+1]

PUB={72:2.09288,128:2.02781,256:1.97752,512:1.92751,1024:1.89534,2048:1.85461,
     4096:1.83452,8192:1.80616,16384:1.78636,32768:1.76683,65536:1.75337,115232:1.73847,
     240000:1.7274,5000000:1.68288}
EXACT={0:-0.5,1:0.125,2:-0.25,3:15/128,5:-47/1024,6:-1/16,7:987/32768,9:-3673/262144,10:1/32}

def area_ladder(b, ns):
    acc=0.0; out={}; ns=set(ns); mx=max(ns)
    for m in range(1,mx+1):
        acc+=m*b[m]*b[m]
        if m in ns: out[m]=math.pi*(1-acc)
    return out

if __name__=="__main__":
    # Verifier / recompute-check. Recomputes b_m and the Gronwall ladder FROM SCRATCH and
    # asserts (a) b_m match the published exact rationals, (b) A_n reproduces every published
    # checkpoint up to 65536 to 5-decimal agreement. Exit 1 on any failure.
    N=int(sys.argv[1]) if len(sys.argv)>1 else 65536
    t=time.time(); b=compute_b(N); dt=time.time()-t
    fails=0
    print(f"# series_newton verifier -- N={N}  time={dt:.1f}s")
    print("# (1) b_m vs published exact rationals:")
    for m,v in sorted(EXACT.items()):
        if m>N: continue
        ok=abs(b[m]-v)<1e-10
        fails+=0 if ok else 1
        print(f"  {'ok ' if ok else 'FAIL'}  b_{m}={b[m]:+.12f}  exact {v:+.12f}")
    print("# (2) Gronwall bound A_n = pi(1 - sum m b_m^2) vs published checkpoints:")
    lad=area_ladder(b,[n for n in PUB if n<=N]+[N])
    for n in sorted(lad):
        p=PUB.get(n)
        if p is None:
            print(f"  new   n={n:>8}  A_n={lad[n]:.6f}   (computed here, beyond the published table)")
            continue
        # published values are printed to 4-5 decimals; agreement to < 5e-5 is "to the digit"
        ok = abs(lad[n]-p) < 5e-5
        fails+=0 if ok else 1
        print(f"  {'ok ' if ok else 'FAIL'}  n={n:>8}  A_n={lad[n]:.6f}   published {p:.5f}   |diff|={abs(lad[n]-p):.1e}")
    print(f"# RESULT: {'PASS' if fails==0 else str(fails)+' FAILURE(S)'}")
    sys.exit(1 if fails else 0)