validate_newton.py — src/21-gravity

1"""Newtonian layer: lattice Green's function, inverse-square law, coherent additivity.23Validates: Lemma (stationary bias field)  u = m/(4 pi kappa r) on Z^3, via the exact4           Montroll representation G(x) = int_0^inf e^{-6t} I_x(2t) I_0(2t)^2 dt;5           Theorem (Newtonian limit)      U(r) ~ -G m1 m2 / r  =>  F ~ 1/r^2;6           Lemma (coherent additivity)    locked ensemble couples as m, unlocked as sqrt(m).7"""8import numpy as np9import mpmath as mp1011ok = True12mp.mp.dps = 201314def G(r):15    f = lambda t: mp.exp(-6*t)*mp.besseli(r, 2*t)*mp.besseli(0, 2*t)**216    return mp.quad(f, [0, mp.inf])1718print("  exact lattice Green vs 1/(4 pi r):")19rs = [3, 5, 8, 12]20Gv = [float(G(r)) for r in rs]21for r, g in zip(rs, Gv):22    ratio = 4*np.pi*r*g23    print(f"    r={r:2d}: 4*pi*r*G(r) = {ratio:.5f}")24    if abs(ratio - 1) > 1.5/r**2 + 0.005: ok = False   # manuscript claims 1 + O(r^-2)2526slope = np.polyfit(np.log(rs[1:]), np.log(Gv[1:]), 1)[0]   # r >= 5, beyond O(r^-2) zone27print(f"  potential power law: U ~ r^{slope:.4f}  (target -1)")28if abs(slope + 1) > 0.02: ok = False2930rng = np.random.default_rng(1)31m = 10**432locked = abs(np.sum(np.exp(1j*np.zeros(m))))33rand = np.mean([abs(np.sum(np.exp(1j*rng.uniform(0, 2*np.pi, m)))) for _ in range(200)])34print(f"  coupling: locked = {locked:.0f} (=m), unlocked ~ {rand:.0f} (target ~ sqrt(pi m)/2 = {np.sqrt(np.pi*m)/2:.0f})")35if not (abs(locked - m) < 1e-9 and 0.5*np.sqrt(m) < rand < 2*np.sqrt(m)): ok = False3637print("PASS" if ok else "FAIL", "- Newtonian layer")383940# ---- round-02: spanning-tree admissibility (W4) ----41# det L_red = tau(S) (Matrix-Tree); exact solvability mod Omega requires tau != 0 mod Omega.42def _tree_count(n):43    # cycle graph C_n: tau = n (known), verify via reduced-Laplacian determinant over Z44    import fractions45    L = [[0]*n for _ in range(n)]46    for i in range(n):47        L[i][i] = 2; L[i][(i+1)%n] -= 1; L[i][(i-1)%n] -= 148    # reduced: drop row/col 0; integer determinant by fraction-free elimination49    m = [row[1:] for row in L[1:]]50    from fractions import Fraction51    M = [[Fraction(v) for v in row] for row in m]52    det = Fraction(1)53    for c in range(n-1):54        piv = next(r for r in range(c, n-1) if M[r][c] != 0)55        if piv != c:56            M[c], M[piv] = M[piv], M[c]; det = -det57        det *= M[c][c]58        M[c] = [v / M[c][c] for v in M[c]]59        for r in range(c+1, n-1):60            if M[r][c] != 0:61                f = M[r][c]; M[r] = [vr - f*vc for vr, vc in zip(M[r], M[c])]62    return det63for _n in (5, 8, 12):64    assert _tree_count(_n) == _n, (_n, _tree_count(_n))65print("PASS spanning-tree admissibility: det L_red = tau(C_n) = n exactly (n = 5, 8, 12); tau != 0 mod Omega is the finite solvability predicate")66