1import json
 2from pathlib import Path
 3import numpy as np
 4from scipy.linalg import eig
 5
 6
 7def pencil(y,n,s=0):
 8    H0=np.array([[y[s+i+j] for j in range(n)] for i in range(n)],float)
 9    H1=np.array([[y[s+i+j+1] for j in range(n)] for i in range(n)],float)
10    u,sv,vh=np.linalg.svd(H0); r=max(1,int(np.sum(sv>1e-9*sv[0])))
11    U=u[:,:r]; V=vh.T[:,:r]
12    z=eig(U.T@H1@V,U.T@H0@V,right=False)
13    return z[np.isfinite(z)]
14
15
16def winding(values):
17    ph=np.unwrap(np.angle(values))
18    return int(np.rint((ph[-1]-ph[0])/(2*np.pi)))
19
20
21def certify(y,n=4,shifts=(0,3,6,9,12),tol=.10):
22    rec=[]
23    for s in shifts:
24        z=pencil(y,n,s); rec.append(z[np.abs(z)<1.05])
25    ref=rec[0]; chosen=[]
26    for z in ref:
27        support=sum(bool(len(a)) and np.min(np.abs(a-z))<tol for a in rec)
28        if support>=max(2,len(shifts)//2): chosen.append(z)
29    return np.asarray(chosen,complex),rec
30
31
32def fit_predict(y,z,train):
33    if len(z)==0: return np.inf
34    tt=np.arange(train); X=np.stack([zj**tt for zj in z],1)
35    a=np.linalg.lstsq(X,y[:train],rcond=None)[0]
36    all_t=np.arange(len(y)); Xa=np.stack([zj**all_t for zj in z],1)
37    return float(np.mean((np.real(Xa@a)[train:]-y[train:])**2))
38
39
40def one_trial(noise,seed):
41    rng=np.random.default_rng(seed); T=260; t=np.arange(T)
42    true=np.array([.96*np.exp(1j*.31),.88*np.exp(1j*.73)])
43    amp=np.array([1+.2j,.65-.1j]); clean=np.real(sum(a*z**t for a,z in zip(amp,true)))
44    y=clean+noise*rng.normal(size=T)
45    cert,rec=certify(y)
46    # A real sequence requires conjugate completion; reject incomplete candidates.
47    pos=[z for z in cert if z.imag>1e-7 and abs(z)<1]
48    idea=np.concatenate([np.array(pos),np.conj(np.array(pos))]) if pos else np.array([],complex)
49    base=np.array([.75*np.exp(1j*.2),.75*np.exp(1j*.6),.75*np.exp(1j*.9),.75*np.exp(1j*1.2)])
50    persistence=[int(sum(bool(a.size) and np.min(np.abs(a-z))<.10 for a in rec)) for z in pos]
51    return {"certified_pairs":len(pos),"persistence":persistence,
52            "idea_mse":fit_predict(y,idea,100),"baseline_mse":fit_predict(y,base,100)}
53
54
55def run():
56    # Pure sequence verifies the generalized-eigenvalue claim to numerical precision.
57    t=np.arange(40); true=np.array([.96*np.exp(1j*.31),.88*np.exp(1j*.73)])
58    clean=np.real((1+.2j)*true[0]**t+(.65-.1j)*true[1]**t)
59    z=pencil(clean,4,0)
60    eig_err=max(float(np.min(np.abs(z-a))) for a in true)
61    p=.8*np.exp(2j*np.pi*np.arange(512)/512)
62    q=np.prod(1-true[:,None]*p[None,:],axis=0)
63    math={"noiseless_max_eigenvalue_error":eig_err,"winding_radius_0.8":winding(q),"expected_winding":0,"contour_margin":float(np.min(np.abs(q)))}
64    trials={}
65    for noise in (.01,.035,.10):
66        vals=[one_trial(noise,s) for s in (7,17,27,37,47)]
67        trials[str(noise)]={"mean":{k:float(np.mean([v[k] if np.isscalar(v[k]) else np.mean(v[k]) for v in vals])) for k in ("certified_pairs","idea_mse","baseline_mse")},"trials":vals}
68    out={"math_check":math,"noise_sweep":trials,"seed_set":[7,17,27,37,47]}
69    Path('results.json').write_text(json.dumps(out,indent=2)); print(json.dumps(out,indent=2))
70
71if __name__=='__main__': run()