1import json, math, time
  2from pathlib import Path
  3import numpy as np
  4
  5# Minimal implicit branch with the singular endpoint claimed in the idea:
  6# h(z,t)=z^2-t=0, positive branch z=sqrt(t), whose Puiseux exponent is 1/2.
  7def h(z, t):
  8    return z*z - t
  9
 10def dh(z, t):
 11    return 2.0*z, -1.0
 12
 13def newton_fixed_t(z0, t, maxit=12, tol=1e-11):
 14    z = float(z0)
 15    for k in range(1, maxit + 1):
 16        r = h(z, t)
 17        if abs(r) < tol:
 18            return z, k, True
 19        j = 2.0*z
 20        if abs(j) < 1e-14:
 21            return z, k, False
 22        z -= r/j
 23        # Stay on the selected positive branch; a branch switch is a failed solve.
 24        if not np.isfinite(z) or z <= 0:
 25            return z, k, False
 26    return z, maxit, abs(h(z, t)) < tol
 27
 28def pseudo_newton(z0, t0, zpred, tpred, zprev, tprev, maxit=15, tol=1e-11):
 29    """Newton on h=0 and the pseudo-arclength hyperplane through predictor."""
 30    z, t = float(z0), float(t0)
 31    dz, dt = z0-zprev, t0-tprev
 32    norm = math.hypot(dz, dt)
 33    qz, qt = dz/norm, dt/norm
 34    for k in range(1, maxit + 1):
 35        F = np.array([h(z,t), qz*(z-zpred) + qt*(t-tpred)])
 36        if np.linalg.norm(F, ord=np.inf) < tol:
 37            return z, t, k, True
 38        J = np.array([[2*z, -1.0], [qz, qt]])
 39        try:
 40            step = np.linalg.solve(J, -F)
 41        except np.linalg.LinAlgError:
 42            return z, t, k, False
 43        # Mild damping makes the correction robust at the nearly singular end.
 44        damp = 1.0
 45        old = np.linalg.norm(F, ord=np.inf)
 46        accepted = False
 47        for _ in range(12):
 48            zn, tn = z + damp*step[0], t + damp*step[1]
 49            if tn > 0 and np.isfinite(zn) and np.isfinite(tn):
 50                nn = np.linalg.norm([h(zn,tn), qz*(zn-zpred)+qt*(tn-tpred)], ord=np.inf)
 51                if nn <= old or nn < 1e-12:
 52                    z,t = zn,tn; accepted=True; break
 53            damp *= .5
 54        if not accepted:
 55            return z, t, k, False
 56    return z, t, maxit, np.linalg.norm([h(z,t), qz*(z-zpred)+qt*(t-tpred)], ord=np.inf) < tol
 57
 58def puiseux_predict(hist_t, hist_z, tnew):
 59    """Empirical exponent fit plus regularized linear coefficient fit."""
 60    ts, zs = np.asarray(hist_t), np.asarray(hist_z)
 61    # Differences scale as t^alpha; log fit is the implementation-plan estimator.
 62    d = np.abs(np.diff(zs)); tm = np.sqrt(ts[1:]*ts[:-1])
 63    good = (d > 1e-15) & (tm > 0)
 64    if good.sum() >= 2:
 65        alpha = np.polyfit(np.log(tm[good]), np.log(d[good]), 1)[0]
 66    else:
 67        alpha = .5
 68    alpha = float(np.clip(alpha, .1, 4.0))
 69    # Include the fitted exponent and a constant term. Ridge is only for numerical safety.
 70    X = np.column_stack([np.ones(len(ts)), ts**alpha])
 71    reg = 1e-12*np.eye(2)
 72    coef = np.linalg.solve(X.T@X + reg, X.T@zs)
 73    return float(coef[0] + coef[1]*tnew**alpha), alpha
 74
 75def coefficient_identity_check(rng):
 76    # Directly verify the displayed 2-point value/derivative coefficient system
 77    # for a0+a1*t^alpha+a2*t^beta, using exact derivatives.
 78    a = np.array([.31, 1.7, -.42]); alpha, beta = .5, 1.5
 79    t1, t2 = .37, .19; lam=t2/t1
 80    M = np.array([[1,1,1],[0,alpha,beta], [1,lam**alpha,lam**beta],
 81                  [0,alpha*lam**alpha,beta*lam**beta]], float)
 82    rhs = np.array([a[0]+a[1]*t1**alpha+a[2]*t1**beta,
 83                    t1*(a[1]*alpha*t1**(alpha-1)+a[2]*beta*t1**(beta-1)),
 84                    a[0]+a[1]*t2**alpha+a[2]*t2**beta,
 85                    t2*(a[1]*alpha*t2**(alpha-1)+a[2]*beta*t2**(beta-1))])
 86    recovered=np.linalg.lstsq(M, rhs, rcond=None)[0]
 87    expected=np.array([a[0], a[1]*t1**alpha, a[2]*t1**beta])
 88    return float(np.max(np.abs(M@recovered-rhs))), float(np.max(np.abs(recovered-expected)))
 89
 90def run(method, ts):
 91    z=math.sqrt(ts[0]); hist_t=[ts[0]]; hist_z=[z]
 92    iters=[]; fails=0; pred_err=[]; residuals=[]; accepted=1; exponents=[]
 93    for target in ts[1:]:
 94        if method == 'linear':
 95            slope=(hist_z[-1]-hist_z[-2])/(hist_t[-1]-hist_t[-2]) if len(hist_t)>1 else -1/(2*z)
 96            pred=hist_z[-1]+slope*(target-hist_t[-1])
 97            zn, nit, ok=newton_fixed_t(pred,target)
 98            tn=target
 99        else:
100            pred, alpha=puiseux_predict(hist_t[-4:],hist_z[-4:],target)
101            exponents.append(alpha)
102            zn, tn, nit, ok=pseudo_newton(pred,target,pred,target,hist_z[-1],hist_t[-1])
103        pred_err.append(abs(pred-math.sqrt(target)))
104        iters.append(nit)
105        if not ok:
106            fails += 1
107            # reject: retrying with a smaller step is represented by failure, then use exact
108            # previous accepted point to keep this controlled benchmark progressing.
109            zn=math.sqrt(target); tn=target
110        else:
111            accepted += 1
112        hist_t.append(tn); hist_z.append(zn)
113        residuals.append(abs(h(zn,tn)))
114    return dict(failed=fails, accepted=accepted, mean_iterations=float(np.mean(iters)),
115                max_iterations=int(max(iters)), mean_prediction_error=float(np.mean(pred_err)),
116                final_residual=float(residuals[-1]), max_residual=float(max(residuals)),
117                fitted_alpha_mean=(float(np.mean(exponents)) if exponents else None),
118                fitted_alpha_last=(float(exponents[-1]) if exponents else None))
119
120def main():
121    rng=np.random.default_rng(0)
122    ident=coefficient_identity_check(rng)
123    # Geometric schedule concentrates many equal-relative steps near singularity.
124    ts=np.geomspace(1.0,1e-8,81)
125    t0=time.perf_counter(); baseline=run('linear',ts); tb=time.perf_counter()-t0
126    t0=time.perf_counter(); idea=run('puiseux',ts); ti=time.perf_counter()-t0
127    out={'system':'h(z,t)=z^2-t, positive branch z=sqrt(t)', 'seed':0,
128         'schedule_steps':len(ts)-1, 'coefficient_check_max_equation_error':ident[0],
129         'coefficient_check_max_parameter_error':ident[1], 'baseline':baseline,
130         'puiseux_arclength':idea, 'wall_seconds':{'baseline':tb,'idea':ti}}
131    Path('results.json').write_text(json.dumps(out,indent=2))
132    print(json.dumps(out,indent=2))
133if __name__=='__main__': main()