transient.py
8.1 kB · python · 252 lines
1import argparse2import json3import os4import resource5import sys6import time7from math import comb89from mpmath import mp, mpf, cos, log, pi, floor, nint1011# DIGIT POLYNOMIAL1213def digit_poly(D):14 base = [0] * (2 * D - 1)15 for k in range(D):16 base[2 * k] = comb(D - 1, k)17 out = [0] * (2 * D + 1)18 for j, v in enumerate(base):19 out[j] += v20 out[j + 1] += D * v21 out[j + 2] += v22 return out2324def fill_value(D):25 return 2 ** (D - 1) * (D + 2)2627def half_width(D):28 return (D - 1) // 22930def coefficient(P, x):31 return P[x] if 0 <= x < len(P) else 03233def sign_vector(D):34 return [(3 if c % 3 == 0 else 0) - 1 for c in range(-half_width(D), half_width(D) + 1)]3536# STRUCTURAL ASSERTIONS3738def assert_nonnegative(D):39 assert all(v >= 0 for v in digit_poly(D)), "digit poly has a negative coefficient at D=%d" % D4041def assert_window_closed(D):42 P = digit_poly(D)43 H = half_width(D)44 for c in range(-H, H + 1):45 for s, val in enumerate(P):46 if not val or (c + D - s) % 3:47 continue48 assert abs((c + D - s) // 3) <= H, "carry window escapes at D=%d" % D4950def assert_palindrome(D):51 P = digit_poly(D)52 assert all(P[j] == P[2 * D - j] for j in range(2 * D + 1)), "digit poly not palindromic at D=%d" % D5354def assert_colsum(D):55 P = digit_poly(D)56 H = half_width(D)57 f = fill_value(D)58 for c in range(-H, H + 1):59 cs = sum(coefficient(P, c + D - 3 * cp) for cp in range(-H, H + 1))60 v = (3 if c % 3 == 0 else 0) - 161 assert f - 3 * cs == (D - 1) * v, "colsum identity fails at D=%d c=%d" % (D, c)6263# FOLDED CORE, COLUMNS OF THE TRANSPOSE6465def folded_columns(D):66 P = digit_poly(D)67 H = half_width(D)68 cols = []69 for c in range(H + 1):70 col = []71 for cp in range(H + 1):72 v = coefficient(P, c + D - 3 * cp)73 if c:74 v += coefficient(P, -c + D - 3 * cp)75 if v:76 col.append((cp, v))77 cols.append(col)78 return cols7980def functional(D):81 out = []82 for c in range(half_width(D) + 1):83 v = (3 if c % 3 == 0 else 0) - 184 out.append(v if c == 0 else 2 * v)85 return out8687# UNFOLDED CENSUS, THE INDEPENDENT ROUTE8889def census_V(D, top):90 P = digit_poly(D)91 H = half_width(D)92 n = 2 * H + 193 step = []94 for cp in range(-H, H + 1):95 step.append([(c + H, coefficient(P, c + D - 3 * cp))96 for c in range(-H, H + 1) if coefficient(P, c + D - 3 * cp)])97 u = [0] * n98 u[H] = 199 v = sign_vector(D)100 out = [sum(v[i] * u[i] for i in range(n))]101 for _ in range(top):102 w = [sum(x * u[j] for j, x in row) for row in step]103 u = w104 out.append(sum(v[i] * u[i] for i in range(n)))105 return out106107# UNFOLDED ROW SWEEP, THE INDEPENDENT ROUTE108109def unfolded_rows(D, top):110 P = digit_poly(D)111 H = half_width(D)112 n = 2 * H + 1113 cols = []114 for c in range(-H, H + 1):115 cols.append([(cp + H, coefficient(P, c + D - 3 * cp))116 for cp in range(-H, H + 1) if coefficient(P, c + D - 3 * cp)])117 r = sign_vector(D)118 out = [r]119 for _ in range(top):120 r = [sum(x * r[j] for j, x in col) for col in cols]121 out.append(r)122 return out123124# ROW SWEEP125126def row_sweep(D, cap):127 assert_nonnegative(D)128 assert_window_closed(D)129 assert_palindrome(D)130 cols = folded_columns(D)131 r = functional(D)132 n = len(r)133 last_neg = -1134 for k in range(1, cap + 1):135 nxt = [0] * n136 for c in range(n):137 s = 0138 for cp, w in cols[c]:139 x = r[cp]140 if x:141 s += w * x142 nxt[c] = s143 r = nxt144 if r[0] < 0:145 last_neg = k146 lo = min(r)147 if lo >= 0:148 return k, last_neg, lo > 0149 raise AssertionError("no nonnegative row within cap at D=%d" % D)150151# CROSSING LEVEL152153def k_star(D, levels=400):154 log_r = mpf(0)155 s = mpf(0)156 for i in range(2, levels):157 a = cos(pi / mpf(3) ** i)158 b = cos(2 * pi / mpf(3) ** i)159 log_r += log(a / b)160 s += log((D - 2 * a) / (D - 2)) - log((D + 2 * b) / (D + 2))161 return ((D - 1) * log_r + s) / log(mpf(D + 2) / mpf(D - 2))162163def l_zero(D, guard):164 ks = k_star(D)165 near = abs(ks - nint(ks))166 assert near > guard, "K* within %s of an integer at D=%d" % (guard, D)167 L0 = int(floor(ks))168 return (L0 if L0 % 2 else L0 - 1), ks, near169170# SELF TEST171172def selftest():173 for D in range(4, 62, 2):174 assert_colsum(D)175 print("colsum identity fill - 3 colsum(c) = (D-1) v_c bites at even D = 4..60; the theorem is prop:mass")176 for D in range(4, 37, 2):177 H = half_width(D)178 cols = folded_columns(D)179 r = functional(D)180 n = len(r)181 top = 4 * D + 8182 mine = [r[0]]183 folded = [r]184 for _ in range(top):185 r = [sum(w * r[cp] for cp, w in cols[c]) for c in range(n)]186 folded.append(r)187 mine.append(r[0])188 assert mine == census_V(D, top), "folded sweep disagrees with the census at D=%d" % D189 full = unfolded_rows(D, top)190 for k in range(top + 1):191 a, b = folded[k], full[k]192 assert all(b[H - c] == b[H + c] for c in range(H + 1)), "unfolded row not symmetric at D=%d k=%d" % (D, k)193 assert all(a[c] == b[H + c] * (1 if c == 0 else 2) for c in range(H + 1)), \194 "folded row disagrees with the unfolded row at D=%d k=%d" % (D, k)195 print("folded sweep V(L) = 3 m0(L) - b(L) against the unfolded census: even D = 4..36, L <= 4D+8, exact")196 print("folded row against the unfolded row on all of S, entrywise: even D = 4..36, L <= 4D+8, exact")197198# RUN199200def run(lo, hi, cap_slope, out_path, budget):201 rows = json.load(open(out_path)) if out_path and os.path.exists(out_path) else []202 done = {row["D"] for row in rows}203 guard = mpf(10) ** -20204 start = time.time()205 for D in range(lo, hi + 1, 2):206 if D in done:207 continue208 L0, ks, near = l_zero(D, guard)209 t0 = time.time()210 t, last_neg, strict = row_sweep(D, int(cap_slope * D * D) + 200)211 dt = time.time() - t0212 rows.append({"D": D, "Lstar": last_neg, "t": t, "L0": L0, "strict": strict,213 "kstar": mp.nstr(ks, 14), "gap": mp.nstr(near, 6),214 "agree": last_neg == L0, "seconds": round(dt, 3)})215 print("D=%3d L*=%5d L0=%5d t=%5d K*=%-16s %s %s %.2fs"216 % (D, last_neg, L0, t, mp.nstr(ks, 12), "ok" if last_neg == L0 else "MISS",217 "strict" if strict else "SLACK", dt), flush=True)218 if out_path:219 json.dump(rows, open(out_path, "w"), indent=1)220 if budget and time.time() - start > budget:221 print("budget wall after D=%d" % D, flush=True)222 break223 return rows224225def main():226 ap = argparse.ArgumentParser()227 ap.add_argument("--lo", type=int, default=6)228 ap.add_argument("--hi", type=int, default=120)229 ap.add_argument("--cap-slope", type=float, default=0.08)230 ap.add_argument("--out", default="")231 ap.add_argument("--budget", type=float, default=0.0)232 ap.add_argument("--selftest", action="store_true")233 a = ap.parse_args()234 mp.dps = 60235 if a.selftest:236 selftest()237 rows = run(a.lo, a.hi, a.cap_slope, a.out, a.budget)238 bad = [r["D"] for r in rows if not r["agree"]]239 slack = [r["D"] for r in rows if not r["strict"]]240 gap = min(rows, key=lambda r: float(r["gap"]))241 steps = sorted({r["t"] - r["Lstar"] for r in rows})242 print("rows %d agree %d misses %s" % (len(rows), sum(r["agree"] for r in rows), bad or "none"))243 print("t - L* values %s non-strict rows %s" % (steps, slack or "none"))244 print("least K* distance to an integer: %s at D = %d" % (gap["gap"], gap["D"]))245 peak = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss246 peak = peak / 1048576.0 if sys.platform == "darwin" else peak / 1024.0247 print("total sweep seconds %.1f, deepest row D = %d at %.1f s, peak resident %.1f MB"248 % (sum(r["seconds"] for r in rows), rows[-1]["D"], rows[-1]["seconds"], peak))249 return 0 if not bad and not slack else 1250251if __name__ == "__main__":252 sys.exit(main())