zeta_shadow.py
13.3 kB · python · 329 lines
1import os2import sys3import time45import mpmath as mp67HERE = os.path.dirname(os.path.abspath(__file__))8sys.path.insert(0, os.path.join(HERE, "..", "design-zeta"))910from design_zeta import Design1112# THE LADDER1314LADDER = [(5, (0, 1)), (4, (0, 1)), (3, (0, 1)), (4, (0, 1, 2)), (5, (0, 1, 2, 3)),15 (10, tuple(range(8))), (10, tuple(range(9))), (20, tuple(range(19))),16 (50, tuple(range(49))), (2, (0, 1))]1718NEW = [(20, tuple(range(19))), (50, tuple(range(49)))]1920OLD = [d for d in LADDER if d not in NEW]2122SETS = {"ladder": LADDER, "new": NEW, "old": OLD}2324NZ = 122526CAPZ = {(20, tuple(range(19))): 6, (50, tuple(range(49))): 6}2728DIFF = mp.mpf("1e-6")2930def line(*a):31 print(" ".join(str(x) for x in a), flush=True)3233def tag(q, F):34 miss = tuple(d for d in range(q) if d not in F)35 if not miss:36 return "base " + str(q) + " full"37 if len(miss) <= len(F):38 return "base " + str(q) + " missing {" + ",".join(str(d) for d in miss) + "}"39 return "base " + str(q) + " {" + ",".join(str(d) for d in F) + "}"4041def med(xs):42 v = sorted(xs)43 n = len(v)44 if n % 2:45 return v[n // 2]46 return (v[n // 2 - 1] + v[n // 2]) / 24748def words(q, F, L):49 cur = [0]50 for _ in range(L):51 cur = [q * m + a for m in cur for a in F]52 return cur5354# THE PRINCIPAL FIBRE5556def fibre(q, F, L):57 return (mp.mpf(len(F)) / q) ** L5859def arc(q, F, L):60 N = q ** L61 tot = mp.mpf(1) / N62 for n in words(q, F, L):63 if n:64 tot += mp.sin(mp.pi * mp.mpf(n) / N) / (mp.pi * n)65 return tot6667def check(q, F, L):68 mp.mp.dps = 4069 N = q ** L70 D = set(words(q, F, L))71 G = [mp.mpf(0)] * N72 for a in range(N):73 G[a] = mp.fsum([mp.e ** (2j * mp.pi * m * a / N) for m in D])74 worst = mp.mpf(0)75 for n in range(N):76 v = mp.fsum([G[a] * mp.e ** (-2j * mp.pi * n * a / N) for a in range(N)]) / N77 worst = max(worst, abs(v - (1 if n in D else 0)))78 return worst, G[0] / N7980def mass(which):81 line("MASS the principal fibre of the position identity. zeta_(F,L)(s) = q^(-L) sum_(a mod q^L)")82 line(" G_L(a/q^L) S_L(s, a/q^L) splits the LEVEL-L polynomial against the partial sum of zeta to")83 line(" q^L, and the a = 0 fibre carries G_L(0)/q^L = (k/q)^L exactly. The level weight falls to 0,")84 line(" so no level is forced: the L = 1 reading c = k/q is a choice and the sweep tests it.")85 line(" ARC is the continuous form, int over abs(t) < 1/(2 q^L) of G_L, exact as a sinc sum;")86 line(" the full set's own readings run toward Si(pi)/pi =", mp.nstr(mp.si(mp.pi) / mp.pi, 10))87 for q, F in which:88 k = len(F)89 c = mp.mpf(k) / q90 if k ** 2 <= 4096 and q ** 2 <= 4096:91 w, f = check(q, F, 2)92 line(" IDENTITY", tag(q, F), "L 2 largest abs 1_(D_L) - inverse transform",93 mp.nstr(w, 4), "a = 0 fibre", mp.nstr(f, 12), "against (k/q)^2",94 mp.nstr(c ** 2, 12))95 row = [tag(q, F), "k/q", mp.nstr(c, 10), "alpha", mp.nstr(mp.log(k) / mp.log(q), 10)]96 for L in (1, 2, 3):97 if k ** L > 200000:98 break99 a = arc(q, F, L)100 row += ["| L", L, "fibre", mp.nstr(fibre(q, F, L), 8), "arc", mp.nstr(a, 8),101 "arc/fibre", mp.nstr(a / fibre(q, F, L), 8)]102 line(" ", *row)103104# THE FIRST-ORDER SHADOW105106def deriv(d, s):107 return (d.zeta(s + DIFF)[0] - d.zeta(s - DIFF)[0]) / (2 * DIFF)108109def pole_gap(d, s):110 per = 2 * mp.pi / mp.log(d.q)111 best = mp.inf112 for m in range(0, 3):113 x = d.alpha - m114 j = mp.nint(mp.im(s) * mp.log(d.q) / (2 * mp.pi))115 for jj in (j - 1, j, j + 1):116 best = min(best, abs(s - mp.mpc(x, jj * per)))117 return best118119def newton(d, s0, trust=mp.mpf("0.6"), cap=14):120 s = s0121 for _ in range(cap):122 v, _ = d.zeta(s)123 g = deriv(d, s)124 if g == 0:125 return None126 step = v / g127 if abs(step) > trust:128 step = step * trust / abs(step)129 s = s - step130 if abs(step) < mp.mpf("1e-28"):131 break132 v, _ = d.zeta(s)133 if abs(v) > mp.mpf("1e-16"):134 return None135 if abs(s - s0) > mp.mpf("1.5") or pole_gap(d, s) < mp.mpf("0.02"):136 return None137 return s138139def predict(which, nz):140 line("PREDICT at a zeta zero rho_0 one has zeta(rho_0) = 0, so for ANY constant c the split")141 line(" zeta_F = c zeta + E_F gives E_F(rho_0) = zeta_F(rho_0), and a zero of zeta_F near rho_0")142 line(" sits at rho_0 - zeta_F(rho_0)/(c zeta'(rho_0)) to first order. STEP is the constant-free")143 line(" reading c zeta' -> zeta_F', that is Newton's own first step -zeta_F(rho_0)/zeta_F'(rho_0),")144 line(" and it is the primary column. PRED is the same law at the L = 1 fibre reading c = k/q.")145 line(" Both see no design zero. FOUND is the zero Newton reaches from rho_0, accepted only at")146 line(" abs(zeta_F) < 1e-16, within 1.5 of rho_0 and 0.02 clear of the pole lattice; a miss is a")147 line(" zero the trust region does not reach and its rung's medians are conditioned on that.")148 out = {}149 for q, F in which:150 d = Design(q, F)151 c = mp.mpf(len(F)) / q152 n_use = CAPZ.get((q, F), nz)153 line(" ", tag(q, F), "alpha", mp.nstr(d.alpha, 10), "k/q", mp.nstr(c, 8),154 "zeros", n_use, "top gamma", mp.nstr(mp.im(mp.zetazero(n_use)), 12))155 rows = []156 for n in range(1, n_use + 1):157 g = mp.im(mp.zetazero(n))158 rho = mp.mpc(mp.mpf("0.5"), g)159 v, eb = d.zeta(rho)160 if abs(v) < mp.mpf("1e-20"):161 line(" n", n, "gamma", mp.nstr(g, 12), "abs zeta_F", mp.nstr(abs(v), 4),162 "CONTROL E_F = 0 and both predicted offsets are 0")163 continue164 zp = mp.zeta(rho, derivative=1)165 fp = deriv(d, rho)166 cpl = fp / zp167 step = -v / fp168 pred = -v / (c * zp)169 s = newton(d, rho)170 if s is None:171 line(" n", n, "gamma", mp.nstr(g, 12), "absE", mp.nstr(abs(v), 8),172 "coupling", mp.nstr(cpl, 8), "step", mp.nstr(abs(step), 8),173 "pred", mp.nstr(abs(pred), 8), "FOUND none")174 continue175 off = s - rho176 rows.append((step, pred, off, cpl, abs(v), eb, s))177 line(" n", n, "gamma", mp.nstr(g, 12), "absE", mp.nstr(abs(v), 8),178 "coupling", mp.nstr(cpl, 8), "off", mp.nstr(off, 8),179 "step ratio", mp.nstr(abs(step) / abs(off), 8),180 "pred ratio", mp.nstr(abs(pred) / abs(off), 8),181 "abs Im off / abs Re off", mp.nstr(abs(mp.im(off)) / abs(mp.re(off)), 8))182 out[(q, F)] = rows183 if rows:184 r0 = sorted(abs(st) / abs(o) for st, p, o, cl, av, eb, s in rows)185 r1 = sorted(abs(p) / abs(o) for st, p, o, cl, av, eb, s in rows)186 iz = sorted(abs(mp.im(o)) / abs(mp.re(o)) for st, p, o, cl, av, eb, s in rows)187 line(" SUMMARY located", len(rows), "of", n_use,188 "| STEP median ratio", mp.nstr(med(r0), 8), "band [", mp.nstr(r0[0], 6), ",",189 mp.nstr(r0[-1], 6), "] largest abs ratio - 1",190 mp.nstr(max(abs(x - 1) for x in r0), 6),191 "| PRED median ratio", mp.nstr(med(r1), 8), "band [", mp.nstr(r1[0], 6), ",",192 mp.nstr(r1[-1], 6), "] largest abs ratio - 1",193 mp.nstr(max(abs(x - 1) for x in r1), 6),194 "| median abs coupling - k/q",195 mp.nstr(med([abs(cl - c) for st, p, o, cl, av, eb, s in rows]), 8),196 "median abs coupling - 1",197 mp.nstr(med([abs(cl - 1) for st, p, o, cl, av, eb, s in rows]), 8),198 "| median abs off",199 mp.nstr(med([abs(o) for st, p, o, cl, av, eb, s in rows]), 8),200 "median absE", mp.nstr(med([av for st, p, o, cl, av, eb, s in rows]), 8),201 "largest ladder bound", mp.nstr(max(eb for st, p, o, cl, av, eb, s in rows), 4),202 "| abs Im off / abs Re off band [", mp.nstr(iz[0], 6), ",",203 mp.nstr(iz[-1], 6), "]")204 return out205206# THE RUNGS207208ZZ = None209210def ordinates():211 global ZZ212 if ZZ is None:213 ZZ = [mp.im(mp.zetazero(n)) for n in range(1, 40)]214 return ZZ215216def shadow_stat(rows):217 zz = ordinates()218 ims = sorted(mp.im(s) for st, p, o, cl, av, eb, s in rows)219 band = [t for t in zz if ims[0] - 2 < t < ims[-1] + 2]220 if len(band) < 2 or len(ims) < 3:221 return None, None222 null = (band[-1] - band[0]) / (len(band) - 1) / 4223 ds = [min(abs(t - y) for y in zz) for t in ims]224 return sum(ds) / len(ds), null225226def loglog(xs, ys):227 lx = [mp.log(x) for x in xs]228 ly = [mp.log(y) for y in ys]229 n = len(lx)230 mx = sum(lx) / n231 my = sum(ly) / n232 sxy = sum((a - mx) * (b - my) for a, b in zip(lx, ly))233 sxx = sum((a - mx) ** 2 for a in lx)234 b = sxy / sxx235 a = my - b * mx236 res = sum((y - (a + b * x)) ** 2 for x, y in zip(lx, ly))237 tot = sum((y - my) ** 2 for y in ly)238 return b, 1 - res / tot239240def rungs(which, nz):241 got = predict(which, nz)242 line("RATE the median paired offset against the missing-digit density m/q = 1 - k/q and against")243 line(" 1 - alpha. The two normalisations differ only by log q, so the ladder discriminates them")244 line(" only over the range of (1 - alpha)/(m/q), which is printed. Every median is over the")245 line(" zeros the trust region reaches, and the rungs do not share one height.")246 xs, ys, zs = [], [], []247 for q, F in which:248 rows = got.get((q, F))249 if not rows:250 continue251 a = mp.log(len(F)) / mp.log(q)252 mq = 1 - mp.mpf(len(F)) / q253 if mq == 0:254 continue255 m = med([abs(o) for st, p, o, cl, av, eb, s in rows])256 xs.append(mq)257 ys.append(m)258 zs.append(1 - a)259 line(" ", tag(q, F), "alpha", mp.nstr(a, 10), "n", len(rows), "1-alpha",260 mp.nstr(1 - a, 8), "m/q", mp.nstr(mq, 8), "(1-alpha)/(m/q)",261 mp.nstr((1 - a) / mq, 8), "median absE",262 mp.nstr(med([av for st, p, o, cl, av, eb, s in rows]), 8),263 "median abs off", mp.nstr(m, 8), "/(m/q)", mp.nstr(m / mq, 8),264 "/(1-alpha)", mp.nstr(m / (1 - a), 8), "| median abs Im off",265 mp.nstr(med([abs(mp.im(o)) for st, p, o, cl, av, eb, s in rows]), 8),266 "median abs Re s - 1/2",267 mp.nstr(med([abs(mp.re(s) - mp.mpf("0.5")) for st, p, o, cl, av, eb, s in rows]), 8))268 if len(xs) > 2:269 b1, r1 = loglog(xs, ys)270 b2, r2 = loglog(zs, ys)271 c1 = sorted(y / x for x, y in zip(xs, ys))272 c2 = sorted(y / x for x, y in zip(zs, ys))273 g = sorted(z / x for x, z in zip(xs, zs))274 line(" FIT a least squares in the logs, a fit and not a theorem: median abs off scales as")275 line(" (m/q)^", mp.nstr(b1, 6), "with R2", mp.nstr(r1, 6), "and as (1-alpha)^",276 mp.nstr(b2, 6), "with R2", mp.nstr(r2, 6))277 line(" the m/q column spans", mp.nstr(c1[-1] / c1[0], 6), "and the 1-alpha column",278 mp.nstr(c2[-1] / c2[0], 6), "while (1-alpha)/(m/q) itself spans",279 mp.nstr(g[-1] / g[0], 6), "over the ladder, so the discrimination is",280 mp.nstr((c2[-1] / c2[0]) / (c1[-1] / c1[0]), 6), "inside that gap")281 line("LADDER the shadow statistic of the family row on the zeros this study pairs: mean distance")282 line(" from a located design ordinate to the nearest zeta ordinate over a quarter of the mean")283 line(" gap between consecutive zeta ordinates in the range. The pairing is zeta-zero-first where")284 line(" the family row's is design-zero-first, so this is a parallel ladder and not that row.")285 for q, F in which:286 rows = got.get((q, F))287 if not rows:288 continue289 a = mp.log(len(F)) / mp.log(q)290 md, null = shadow_stat(rows)291 if md is None:292 continue293 line(" ", tag(q, F), "alpha", mp.nstr(a, 10), "n", len(rows), "mean dist",294 mp.nstr(md, 8), "null gap/4", mp.nstr(null, 8), "ratio to null",295 mp.nstr(md / null, 8))296 line("SHARPNESS the constant-free step against the k/q reading, pooled by the size of the offset")297 pool = []298 for q, F in which:299 for st, p, o, cl, av, eb, s in got.get((q, F)) or []:300 pool.append((abs(o), abs(st) / abs(o), abs(p) / abs(o)))301 for lo, hi in ((0, mp.mpf("0.05")), (mp.mpf("0.05"), mp.mpf("0.1")),302 (mp.mpf("0.1"), mp.mpf("0.2")), (mp.mpf("0.2"), mp.mpf("0.4")),303 (mp.mpf("0.4"), mp.inf)):304 b = [(r0, r1) for u, r0, r1 in pool if lo <= u < hi]305 if not b:306 continue307 line(" abs off in [", mp.nstr(lo, 3), ",", mp.nstr(hi, 3), ") n", len(b),308 "STEP median", mp.nstr(med([x[0] for x in b]), 8), "largest abs ratio - 1",309 mp.nstr(max(abs(x[0] - 1) for x in b), 6), "| PRED median",310 mp.nstr(med([x[1] for x in b]), 8), "largest abs ratio - 1",311 mp.nstr(max(abs(x[1] - 1) for x in b), 6))312313# STUDY314315def main():316 argv = sys.argv[1:]317 verb = argv[0] if argv else "mass"318 which = SETS[argv[1]] if len(argv) > 1 else SETS["ladder"]319 t0 = time.time()320 if verb == "mass":321 mass(which)322 if verb == "predict":323 predict(which, NZ)324 if verb == "rungs":325 rungs(which, NZ)326 line("seconds", round(time.time() - t0, 1))327328if __name__ == "__main__":329 main()