multilayer.py
12.9 kB · python · 309 lines
1import itertools2import subprocess3import sys4import time56import mpmath7import numpy as np8import sympy910NB = 1.4511TWO_PI = 2 * np.pi12CONTRASTS = (1.6, 2.3, 3.5)13DEEP = 0.021415def lay(n, ph):16 c, s = np.cos(ph), np.sin(ph)17 return (c + 0j, 1j * s / n, 1j * n * s, c + 0j)1819def mul(p, q):20 a, b, c, d = p21 e, f, g, h = q22 return (a * e + b * g, a * f + b * h, c * e + d * g, c * f + d * h)2324def num_den(m):25 a, b, c, d = m26 x, y = NB * a + NB * NB * b, c + NB * d27 return x - y, x + y2829def peak(m):30 return np.maximum.reduce([np.abs(z) for z in m])3132def levels(base, digits, top, delta, na):33 m = lay(na, delta)34 dead = np.zeros(delta.shape, bool)35 out = [(m, dead.copy())]36 for k in range(top):37 spacer = lay(NB, np.mod(float(base) ** k * delta, TWO_PI))38 c = np.where(dead, 1.0, np.maximum(peak(m), 1e10) / 1e10)39 unit = tuple(z / c for z in m)40 n = None41 for d in range(base):42 x = unit if d in digits else spacer43 n = x if n is None else mul(n, x)44 big = peak(n)45 grow = len(digits) * np.log(c)46 dead |= np.log(big) + grow >= np.log(1e100)47 scale = np.where(dead, 1 / big, np.exp(np.minimum(grow, 700)))48 m = tuple(z * scale for z in n)49 out.append((m, dead.copy()))50 return out5152def log_rt(m, dead):53 num, _ = num_den(m)54 v = np.log(np.maximum(np.abs(num) / (2 * NB), 1e-14))55 return np.where(dead, np.inf, v)5657def phase_t(m):58 _, den = num_den(m)59 return -np.angle(den)6061def word(base, digits, level):62 w = [1]63 for _ in range(level):64 w = [x if d in digits else 0 for d in range(base) for x in (w if d in digits else [0] * len(w))]65 return w6667def brute(base, digits, level, delta, na):68 mpmath.mp.dps = 4069 out = []70 for dl in delta:71 dl = mpmath.mpf(float(dl))72 c, s = mpmath.cos(dl), mpmath.sin(dl)73 mats = {}74 for n in (mpmath.mpf(na), mpmath.mpf(NB)):75 mats[n] = (c, 1j * s / n, 1j * n * s, c)76 m = (1, 0, 0, 1)77 for x in word(base, digits, level):78 m = mul(m, mats[mpmath.mpf(na) if x else mpmath.mpf(NB)])79 out.append([complex(z) for z in m])80 return tuple(np.array(z) for z in zip(*out))8182def mahler_p(digits):83 c = np.zeros(max(digits) + 1)84 c[list(digits)] = 185 r = np.roots(c[::-1])86 return float(np.sum(np.log(np.maximum(1, np.abs(r)))))8788def mahler_q(digits, n=1 << 15):89 y = np.exp(2j * np.pi * (np.arange(n) + 0.5) / n)90 deg = len(digits) - 191 if deg == 0:92 return 0.093 coef = np.stack([y ** (d - i) for i, d in enumerate(digits)], axis=1)94 lead = coef[:, -1]95 comp = np.zeros((n, deg, deg), complex)96 comp[:, 0, :] = -coef[:, -2::-1] / lead[:, None]97 if deg > 1:98 comp[:, np.arange(1, deg), np.arange(deg - 1)] = 199 r = np.linalg.eigvals(comp)100 return float(np.mean(np.sum(np.log(np.maximum(1, np.abs(r))), axis=1)))101102def q_value(digits, x, y):103 return sum(x ** i * y ** (d - i) for i, d in enumerate(digits))104105X, Y, Z = sympy.symbols("x y z")106107def cyclotomic_one(coeffs):108 poly = sympy.Poly(sum(c * Z ** e for e, c in coeffs), Z)109 return poly.is_cyclotomic or (-poly).is_cyclotomic110111def boyd_class(digits):112 q = sum(X ** i * Y ** (d - i) for i, d in enumerate(digits))113 _, fs = sympy.factor_list(q, X, Y)114 for f, _ in fs:115 terms = sympy.Poly(f, X, Y).terms()116 if len(terms) == 1:117 continue118 e0 = np.array(terms[0][0])119 diffs = [np.array(t[0]) - e0 for t in terms]120 v = next(d for d in diffs if d.any())121 u = v // np.gcd(abs(int(v[0])), abs(int(v[1])))122 steps = []123 for d in diffs:124 if u[0] * d[1] - u[1] * d[0] != 0:125 return False126 steps.append(int(d @ u) // int(u @ u))127 lo = min(steps)128 if not cyclotomic_one([(s - lo, int(c)) for s, (_, c) in zip(steps, terms)]):129 return False130 return True131132def p_cyclotomic(digits):133 _, fs = sympy.factor_list(sum(Z ** d for d in digits), Z)134 return all(sympy.Poly(f, Z).is_cyclotomic for f, _ in fs)135136def fmt(xs, p=3):137 return " ".join(f"{x:.{p}f}" for x in xs)138139def check():140 rng = np.random.default_rng(1)141 delta = rng.uniform(0, np.pi, 8)142 worst = 0.0143 for base, digits, level in [(3, {0, 2}, 6), (4, {0, 1, 3}, 5), (5, {0, 1, 2, 4}, 4), (7, {0, 2, 3, 4, 6}, 4), (9, {0, 1, 4, 7, 8}, 3)]:144 err = 0.0145 for na in CONTRASTS:146 a = brute(base, digits, level, delta, na)147 b = levels(base, digits, level, delta, na)[level][0]148 scale = np.maximum.reduce([np.abs(z) for z in a])149 err = max(err, max(float(np.max(np.abs(x - y) / scale)) for x, y in zip(a, b)))150 worst = max(worst, err)151 print(f"check b={base} D={sorted(digits)} level={level} cells={base ** level} max rel err {err:.1e}")152 assert worst < 1e-12153 h = 1e-8154 for base, digits in [(4, [0, 1, 3]), (7, [0, 2, 3, 4, 6]), (3, [0, 2])]:155 delta = rng.uniform(0, np.pi, 64)156 lv = levels(base, set(digits), 4, delta, NB * (1 + h))157 r0 = num_den(lv[0][0])[0]158 errs = []159 for k in range(1, 5):160 born = np.prod([q_value(digits, np.exp(-2j * (base ** j) * delta), np.exp(-2j * (base ** j) * delta)) for j in range(k)], axis=0)161 got = num_den(lv[k][0])[0] / r0162 errs.append(float(np.max(np.abs(np.abs(got) - np.abs(born)) / np.maximum(1, np.abs(born)))))163 print(f"born derivative b={base} D={digits} levels 1-4 max rel err at h=1e-8: {fmt(errs, 9)}")164 assert max(errs) < 1e-5165166def born():167 rng = np.random.default_rng(2)168 delta = rng.uniform(0, np.pi, 200000)169 na = 1.455170 gaps, ses = [], []171 for base, digits, top in [(4, [0, 1, 3], 7), (5, [0, 1, 2, 4], 6), (7, [0, 1, 5], 6), (3, [0, 2], 8), (6, [0, 1, 3, 4], 6), (7, [0, 2, 3, 4, 6], 6)]:172 lv = levels(base, set(digits), top, delta, na)173 lr = [log_rt(*x) for x in lv]174 steps = [lr[k + 1] - lr[k] for k in range(top)]175 per = [float(np.mean(s)) for s in steps]176 se = float(np.std(lr[top] - lr[0]) / np.sqrt(delta.size) / top)177 slope = float(np.mean(lr[top] - lr[0]) / top)178 mp = mahler_p(digits)179 gaps.append(abs(slope - mp))180 ses.append(se)181 print(f"born nA={na} b={base} D={digits} m(P)={mp:.5f} m(Q)={mahler_q(digits):.5f} per-level {fmt(per, 4)} mean {slope:.4f} se {se:.4f} |mean - m(P)| {abs(slope - mp):.5f}")182 print(f"born: |mean - m(P)| at most {max(gaps[:3]):.5f} on the three with m(P) > 0 and {max(gaps):.5f} on all six designs, se at most {max(ses):.4f}")183184def deep_run(base, digits, top, delta, na):185 lv = levels(base, set(digits), top, delta, na)186 lr = [log_rt(*x) for x in lv]187 r = {key: [] for key in ("drift", "se", "lift", "lemma", "coup", "torus", "tse", "dead")}188 stop = False189 for k in range(1, top):190 a_all = phase_t(lv[k][0])191 b_all = phase_t(lay(NB, np.mod(float(base) ** k * delta, TWO_PI)))192 lt = np.log(np.maximum(np.abs(q_value(digits, np.exp(2j * a_all), np.exp(2j * b_all))), 1e-300))193 r["torus"].append(float(np.mean(lt)))194 r["tse"].append(float(np.std(lt) / np.sqrt(lt.size)))195 r["dead"].append(float(np.mean(lv[k][1])))196 sel = (lr[k] < np.log(DEEP)) & np.isfinite(lr[k + 1])197 stop = stop or sel.sum() < 300198 if stop:199 continue200 al = phase_t(lv[k][0])[sel]201 be = phase_t(lay(NB, np.mod(float(base) ** k * delta[sel], TWO_PI)))202 lq = np.log(np.abs(q_value(digits, np.exp(2j * al), np.exp(2j * be))))203 dd = (lr[k + 1] - lr[k])[sel]204 r["drift"].append(float(np.mean(dd)))205 r["se"].append(float(np.std(dd) / np.sqrt(dd.size)))206 r["lift"].append(float(np.mean(lq)))207 tiny = lr[k][sel] < np.log(1e-3)208 if tiny.sum() >= 50:209 r["lemma"].append(float(np.median(np.abs(dd - lq)[tiny])))210 r["coup"].append(float(np.abs(np.mean(np.exp(2j * (al - be))))))211 return r212213DEEP_DESIGNS = [(4, [0, 1, 3], 12), (5, [0, 1, 2, 4], 10), (7, [0, 1, 5], 10), (7, [0, 2, 3, 4, 6], 10), (9, [0, 1, 4, 7, 8], 9), (8, [0, 2, 3, 4, 5, 7], 9), (3, [0, 2], 14), (7, [0, 1, 3, 4], 10), (7, [0, 2, 4, 6], 10)]214215TABLE = DEEP_DESIGNS[2:]216217def deep():218 rng = np.random.default_rng(3)219 delta = rng.uniform(0, np.pi, 200000)220 cells = lemmas = 0221 worst = track = 0.0222 zero, pos = [], []223 for base, digits, top in DEEP_DESIGNS:224 mq = mahler_q(digits)225 print(f"deep b={base} D={digits} m(P)={mahler_p(digits):.4f} m(Q)={mq:.4f}")226 for na in CONTRASTS:227 r = deep_run(base, digits, top, delta, na)228 tail, tse = r["drift"][-4:], r["se"][-4:]229 mean, se = float(np.mean(tail)), float(np.mean(tse))230 gap = max(abs(a - b) for a, b in zip(r["drift"][-4:], r["lift"][-4:]))231 cells += 1232 track = max(track, gap)233 lemmas += bool(r["lemma"])234 worst = max([worst] + r["lemma"])235 law = f"{max(r['lemma']):.1e}" if r["lemma"] else "n/a, fewer than 50 frequencies below 1e-3"236 print(f" nA={na} deep drift {fmt(r['drift'])} | tail mean {mean:.3f} se {se:.3f} | (tail - m(Q))/se {(mean - mq) / se:.1f}")237 print(f" per-level se {fmt(r['se'])} | deep lift {fmt(r['lift'])} | tail max |drift - lift| {gap:.4f}")238 print(f" step law median err {law} | coupling {fmt(r['coup'], 2)}")239 print(f" torus mean of ln|Q_D| over all frequencies {fmt(r['torus'])} | last {r['torus'][-1]:.3f} se {r['tse'][-1]:.3f} | dark fraction {fmt(r['dead'], 2)}")240 if (base, digits, top) in TABLE:241 (zero if mq < 1e-6 else pos).append((mean, se))242 print(f"deep step law: {lemmas} of {cells} design-contrast cells have 50 or more frequencies below 1e-3, worst median err {worst:.1e}; deep set, last four levels: |drift - lift| at most {track:.4f}")243 print(f"deep zero set: m(Q) = 0 tail |mean| at most {max(abs(m) for m, _ in zero):.3f} (se at most {max(e for _, e in zero):.3f}); m(Q) > 0 tail mean at least {min(m for m, _ in pos):.3f} (se at most {max(e for _, e in pos):.3f})")244245def smyth():246 s = subprocess.run(["gp", "-q"], input="print(3*sqrt(3)/(4*Pi)*lfun(-3,2))", capture_output=True, text=True).stdout.strip()247 mq = mahler_q([0, 1, 3], 1 << 17)248 print(f"smyth constant from PARI {s[:14]} m(Q_(0,1,3)) numeric {mq:.10f} m(P_(0,1,3)) {mahler_p([0, 1, 3]):.10f}")249 assert abs(mq - float(s)) < 1e-8250251def passband(base, digits, top, delta, na):252 lv = levels(base, set(digits), top, delta, na)253 return np.array([float(np.mean(log_rt(*x) < 0)) for x in lv])254255def refute():256 rng = np.random.default_rng(4)257 delta = rng.uniform(0, np.pi, 200000)258 for base, digits in [(7, [0, 2, 3, 4, 6]), (9, [0, 1, 4, 7, 8]), (8, [0, 2, 3, 4, 5, 7]), (7, [0, 2, 4, 6]), (7, [0, 1, 3, 4])]:259 f = passband(base, digits, 12, delta, 2.3)260 ks = np.arange(1, 13)261 pc, bc = p_cyclotomic(digits), boyd_class(digits)262 print(f"refute b={base} D={digits} P cyclotomic {pc} m(P)={mahler_p(digits):.1e} Q generalised cyclotomic {bc} m(Q)={mahler_q(digits):.4f}")263 print(f" f_k sqrt(k) k=1..12 at nA=2.3: {fmt(f[1:] * np.sqrt(ks))}")264 if digits in ([0, 2, 3, 4, 6], [0, 1, 4, 7, 8], [0, 2, 3, 4, 5, 7]):265 assert pc and not bc266267def designs(lo, hi):268 seen = set()269 for base in range(lo, hi + 1):270 for size in range(2, base):271 for rest in itertools.combinations(range(1, base), size - 1):272 d = (0,) + rest273 mirror = tuple(sorted(max(d) - x for x in d))274 key = (base, min(d, mirror))275 if key not in seen:276 seen.add(key)277 yield key278279def census():280 rng = np.random.default_rng(5)281 delta = rng.uniform(0, np.pi, 100000)282 top, half = 12, 6283 agree = total = 0284 flats = {True: [], False: []}285 for base, digits in designs(3, 6):286 digits = list(digits)287 bc = boyd_class(digits)288 mq = mahler_q(digits)289 f = passband(base, digits, top, delta, 2.3)290 flat = f[top] * np.sqrt(top) / (f[half] * np.sqrt(half)) if f[half] > 0 else 0.0291 critical = flat > 0.9292 total += 1293 agree += critical == bc294 flats[bc].append(flat)295 tag = "critical" if critical else "decays"296 print(f"census b={base} D={digits} Q gen-cyclotomic {bc!s:5} m(Q)={mq:.4f} m(P)={mahler_p(digits):.4f} f6={f[half]:.4f} f12={f[top]:.5f} flatness {flat:.3f} {tag}")297 print(f"census designs {total} (bases 3 to 6, up to shift and mirror), flatness verdict matches the Boyd class on {agree}")298 for bc, name in ((True, "Boyd-cyclotomic"), (False, "other")):299 print(f"census {name}: {len(flats[bc])} designs, flatness {min(flats[bc]):.3f} to {max(flats[bc]):.3f}")300301VERBS = {"check": check, "born": born, "smyth": smyth, "deep": deep, "refute": refute, "census": census}302303if __name__ == "__main__":304 import warnings305 warnings.filterwarnings("ignore")306 for v in sys.argv[1:] or list(VERBS):307 t0 = time.time()308 VERBS[v]()309 print(f"[{v} {time.time() - t0:.1f}s]")