"""n35-A — analytic no-go, numerical sweep.
Claim: a single local, stable, canonical scalar sector (tau, eps) with the
relational shift symmetry (tau -> tau + const), read by matter through ANY
local symmetry-respecting coupling, cannot have both a Newtonian static
response chi ~ 1/q^2 and z=1 propagation Re omega ~ q.
Proof sketch: shift symmetry forces a(q)=a2 q^2, c(q)=i c1 q, g_tau(q)=O(q);
z=1 requires b(0)>0 (stability: c1^2 <= a2 b0); chi ~ 1/q^2 requires
det K ~ q^4 i.e. b0=0 AND c1=0. Contradiction. The b0=c1=0 corner has
p=2 and z=2 exactly — it is the n33 model.
Sweep: 5000 random draws of the stable class, all three coupling cases;
count points with p in [1.95,2.05] AND z in [0.95,1.05]. Expected: zero.
Also verifies the identity xi * c_eff = 2 Omega eta to machine precision."""
import numpy as np
rng=np.random.default_rng(7)
qs=np.geomspace(1e-3,1e-2,7)
def chi_and_z(a2,b0,b2,c1,c2,gt1,ge):
    chi=[];rw=[]
    for q in qs:
        a=a2*q*q; b=b0+b2*q*q; c=1j*c1*q+c2*q*q
        K=np.array([[a,c],[np.conj(c),b]])
        if np.linalg.eigvalsh(K).min()<=0: return None
        g=np.array([1j*gt1*q,ge])
        chi.append(float(np.real(np.conj(g)@np.linalg.inv(K)@g)))
        M=np.array([[np.conj(c),b],[-a,-c]])
        w=1j*np.linalg.eigvals(M)
        rw.append(np.abs(np.real(w)).max())
    chi=np.array(chi); rw=np.array(rw)
    p=-np.polyfit(np.log(qs),np.log(np.abs(chi)+1e-300),1)[0]
    z= np.polyfit(np.log(qs),np.log(rw+1e-300),1)[0] if rw.max()>1e-14 else 99
    return p,z
hits=0
for trial in range(5000):
    a2,b2=rng.uniform(0.1,5,2); b0=rng.choice([0.0,rng.uniform(0.01,5)])
    c1=(rng.uniform(-0.99,0.99)*np.sqrt(a2*b0)) if b0>0 else 0.0
    c2=rng.uniform(-1,1)*np.sqrt(a2*b2)*0.9
    gt1,ge=[(0,1),(1,0),(rng.uniform(-1,1),rng.uniform(-1,1))][trial%3]
    r=chi_and_z(a2,b0,b2,c1,c2,gt1,ge)
    if r and 1.95<=r[0]<=2.05 and 0.95<=r[1]<=1.05: hits+=1
print(f"points with Newton AND z=1: {hits} / 5000  (no-go: expect 0)")
r=chi_and_z(2.0,0.0,1.0,0.0,0.0,0,1)
print(f"the n33 corner (b0=c1=0): p={r[0]:.4f}  z={r[1]:.4f}  (Newton, non-propagating)")
Om,eta=0.09,0.5
for U in (1e-3,1e-1):
    xi=np.sqrt(2*Om*eta**2/U); c=np.sqrt(2*Om*U)
    print(f"identity check U={U}: xi*c_eff={xi*c:.12f} vs 2*Omega*eta={2*Om*eta:.12f}")