1import json
  2import math
  3import numpy as np
  4
  5
  6def fejer_coeffs(K):
  7    """p(t)=(1/(K+1))*sum_{j=0}^K t^j."""
  8    c = np.zeros(K + 1)
  9    c[:] = 1.0 / (K + 1)
 10    return c
 11
 12
 13def jackson_coeffs(K):
 14    """A Jackson-type polynomial: normalized square of a Fejer amplitude.
 15
 16    For N=floor(K/2)+1, q(t)=sum_{j=0}^{N-1}t^j/N and p=q^2.
 17    It has p(1)=1, degree <=K, and on the unit circle its magnitude is
 18    the squared Fejer-kernel amplitude, giving quadratic suppression away
 19    from t=1.
 20    """
 21    N = K // 2 + 1
 22    q = np.ones(N) / N
 23    c = np.convolve(q, q)
 24    assert len(c) - 1 <= K
 25    return c
 26
 27
 28def polyval(c, z):
 29    out = np.zeros_like(z, dtype=np.complex128)
 30    for a in c[::-1]:
 31        out = out * z + a
 32    return out
 33
 34
 35def rotation_spectrum(s, n=4096, phase_max=math.pi):
 36    """J eigenvalues=(1+exp(i phi))/2, with dist(1,sigma(J))=s.
 37    Dense phases expose the worst-case envelope rather than lucky zeros.
 38    """
 39    phi0 = 2.0 * math.asin(min(1.0, s))
 40    phis = np.linspace(phi0, phase_max, n)
 41    # Include conjugate pairs implicitly since magnitudes are symmetric.
 42    return (1.0 + np.exp(1j * phis)) / 2.0, phis
 43
 44
 45def residual_ratio(lam, coeffs):
 46    t = 2.0 * lam - 1.0
 47    p = polyval(coeffs, t)
 48    return float(np.max(np.abs((1.0 - lam) * p)))
 49
 50
 51def weighted_residual(lam, coeffs, seed=0):
 52    rng = np.random.default_rng(seed)
 53    y = rng.normal(size=len(lam)) + 1j * rng.normal(size=len(lam))
 54    r0 = np.linalg.norm((1-lam)*y)
 55    p = polyval(coeffs, 2*lam-1)
 56    return float(np.linalg.norm((1-lam)*p*y) / r0)
 57
 58
 59def apply_filter(lam, y, coeffs):
 60    """Horner application of p(2J-I) to a vector, using J diagonalized here."""
 61    t = 2*lam - 1
 62    z = np.zeros_like(y, dtype=np.complex128)
 63    for a in coeffs[::-1]:
 64        z = t*z + a*y
 65    return z
 66
 67
 68def math_checks():
 69    checks = {}
 70    for K in [3, 7, 15, 31]:
 71        f, j = fejer_coeffs(K), jackson_coeffs(K)
 72        checks[str(K)] = {
 73            "fejer_p1_error": abs(polyval(f, np.array([1+0j]))[0]-1),
 74            "jackson_p1_error": abs(polyval(j, np.array([1+0j]))[0]-1),
 75            "jackson_degree": len(j)-1,
 76            "jackson_nonnegative_coeffs": bool(np.min(j) >= 0),
 77        }
 78    # Direct polynomial-vs-iterates identity on a random diagonal operator.
 79    rng = np.random.default_rng(3)
 80    lam = (1 + np.exp(1j*rng.uniform(-math.pi, math.pi, 20))) / 2
 81    y = rng.normal(size=20)+1j*rng.normal(size=20)
 82    K=15; c=fejer_coeffs(K)
 83    direct = apply_filter(lam, y, c)
 84    z=y.copy(); avg=np.zeros_like(y)
 85    for _ in range(K+1):
 86        avg += z/(K+1)
 87        z=(2*lam-1)*z
 88    checks["iterate_identity_error"] = float(np.linalg.norm(direct-avg))
 89    return checks
 90
 91
 92def run_sweep():
 93    rows=[]
 94    # Include both critical sK=O(1) and increasingly gapped regimes.
 95    for K in [7, 15, 31, 63]:
 96        for s in [0.25/K, 0.5/K, 1.0/K, 1.5/K, 2.0/K, 4.0/K, 8.0/K, 16.0/K]:
 97            s=min(s, .95)
 98            lam,_=rotation_spectrum(s)
 99            f=fejer_coeffs(K); j=jackson_coeffs(K)
100            rf=residual_ratio(lam,f); rj=residual_ratio(lam,j)
101            # The rule in the proposal, with candidate K fixed for this row.
102            chosen = "fejer" if s*K < 2 else "jackson"
103            ra = rf if chosen=="fejer" else rj
104            rows.append({"K":K,"s":s,"sK":s*K,"fejer":rf,"jackson":rj,
105                         "adaptive":ra,"choice":chosen,
106                         "jackson_scaled":rj*K*K*s,
107                         "fejer_scaled":rf*(K+1)})
108    return rows
109
110
111def summary(rows):
112    # Ratios are robust envelope quantities on the same dense spectrum.
113    critical=[r for r in rows if r["sK"] <= 1.5]
114    gapped=[r for r in rows if r["sK"] >= 4]
115    # Fit log slope of Jackson residual vs s at each K in gapped range.
116    slopes=[]
117    for K in sorted(set(r["K"] for r in gapped)):
118        q=[r for r in gapped if r["K"]==K]
119        slopes.append(float(np.polyfit(np.log([r["s"] for r in q]), np.log([r["jackson"] for r in q]),1)[0]))
120    # Candidate adaptive degree selection, measuring best of {3,7,15} with rule.
121    adaptive_gain=[]
122    for s in [0.01,0.03,0.06,0.12,0.25,0.5]:
123        vals=[]
124        for K in [3,7,15]:
125            lam,_=rotation_spectrum(s)
126            f=fejer_coeffs(K); j=jackson_coeffs(K)
127            vals.append((residual_ratio(lam,f) if s*K<2 else residual_ratio(lam,j),K))
128        adaptive_gain.append({"s":s,"selected":min(vals)[1],"residual":min(vals)[0]})
129    return {
130      "critical_mean_fejer_times_Kplus1": float(np.mean([r["fejer_scaled"] for r in critical])),
131      "critical_max_fejer_times_Kplus1": float(np.max([r["fejer_scaled"] for r in critical])),
132      "gapped_jackson_K2s_mean": float(np.mean([r["jackson_scaled"] for r in gapped])),
133      "gapped_jackson_log_s_slopes": slopes,
134      "adaptive_candidates": adaptive_gain,
135      "rows":rows
136    }
137
138if __name__ == "__main__":
139    out={"math_checks":math_checks()}
140    rows=run_sweep(); out["summary"]=summary(rows)
141    with open("results.json","w") as f: json.dump(out,f,indent=2)
142    print(json.dumps(out["math_checks"],indent=2))
143    print(json.dumps({k:v for k,v in out["summary"].items() if k!="rows"},indent=2))
144
145
146def adaptive_coeffs(K, gap_hat, baseline_residual=None, candidate_residual=None,
147                    reject_ratio=1.10):
148    """Select Fejer below s*K=2, otherwise Jackson, with residual safeguard."""
149    f = fejer_coeffs(K)
150    if gap_hat * K < 2.0:
151        return f, "fejer", False
152    j = jackson_coeffs(K)
153    rejected = (baseline_residual is not None and candidate_residual is not None
154                and candidate_residual > reject_ratio * baseline_residual)
155    return (f, "fejer", True) if rejected else (j, "jackson", False)
156
157
158def adaptive_apply(lam, y, K, gap_hat, reject_ratio=1.10):
159    """Apply the adaptive rule on a diagonal toy operator."""
160    f, j = fejer_coeffs(K), jackson_coeffs(K)
161    r0 = np.linalg.norm((1-lam)*y)
162    yf = apply_filter(lam, y, f)
163    rf = np.linalg.norm((1-lam)*yf) / r0
164    if gap_hat*K < 2:
165        return yf, "fejer", False, rf
166    yj = apply_filter(lam, y, j)
167    rj = np.linalg.norm((1-lam)*yj) / r0
168    if rj > reject_ratio*rf:
169        return yf, "fejer", True, rf
170    return yj, "jackson", False, rj