carry_free.py
9.1 kB · python · 305 lines
1import math2import os3import sys45sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "mrly-pairing"))67import pairing89# MONOID1011B = 321213DEEP = 181415WIDE = 141617BASES = [(3, 16), (4, 12), (5, 11), (6, 9), (7, 8)]1819RHO3 = 2.2075122021def gens(L):22 out = []23 for m in range(2, 1 << L):24 v = 025 for i in range(L):26 if (m >> i) & 1:27 v |= 1 << (i * B)28 out.append((m.bit_length() - 1, v))29 return out3031def deg(p):32 return (p.bit_length() - 1) // B3334def coeffs(p):35 c = []36 while p:37 c.append(p & ((1 << B) - 1))38 p >>= B39 return c4041def value(p, q):42 v = 043 for a in reversed(coeffs(p)):44 v = v * q + a45 return v4647def show(p):48 t = []49 for i, a in enumerate(coeffs(p)):50 if a == 0:51 continue52 h = "" if a == 1 and i else str(a)53 b = "" if i == 0 else ("x" if i == 1 else "x^" + str(i))54 t.append(h + b)55 return " + ".join(t) if t else "0"5657def monoid(L):58 g = gens(L)59 seen = {1}60 cur = [1]61 while cur:62 nxt = []63 for p in cur:64 d = deg(p)65 for dd, v in g:66 if dd + d >= L:67 break68 r = p * v69 if r not in seen:70 seen.add(r)71 nxt.append(r)72 cur = nxt73 return sorted(seen), g7475# INVERSE7677def nustar(L):78 els, g = monoid(L)79 nu = dict.fromkeys(els, 0)80 nu[1] = 181 for p in els:82 v = nu[p]83 if v == 0:84 continue85 d = deg(p)86 for dd, w in g:87 if dd + d >= L:88 break89 nu[p * w] -= v90 return els, nu9192def power(l):93 return 1 if l == 0 else (-1 if l == 1 else 0)9495def ladder(L):96 els, nu = nustar(L)97 for l in range(L):98 assert nu[1 << (l * B)] == power(l)99 run = 0100 best = 0101 mx, sums, cens = [], [], [0] * L102 top = 0103 d = 0104 for p in els:105 while deg(p) > d:106 sums.append(run + power(d + 1))107 mx.append(max(best, abs(run + power(d + 1))))108 d += 1109 run += nu[p]110 cens[d] += 1111 top = max(top, max(coeffs(p)))112 best = max(best, abs(run))113 sums.append(run + power(L))114 mx.append(max(best, abs(run + power(L))))115 return mx, sums, cens, top116117# BOUND118119def qset(L, els):120 top = [p for p in els if deg(p) < L]121 for q in range(2, 4096):122 hi = max(value(p, q) for p in top)123 if hi < q ** L:124 return q, max(top, key=lambda p: value(p, q - 1))125 return None, None126127def closed(L):128 return next(q for q in range(2, 4096) if (q + 1) ** (L - 1) < q ** L)129130def digits(n, q):131 d = []132 while n:133 d.append(n % q)134 n //= q135 return "".join(str(a) for a in reversed(d))136137def qord(L, els):138 lst = [p for p in els if deg(p) <= L]139 for q in range(2, 4096):140 prev = -1141 bad = None142 for i, p in enumerate(lst):143 v = value(p, q)144 if v <= prev:145 bad = (lst[i - 1], p)146 break147 prev = v148 if bad is None:149 return q, wit(lst, q - 1)150 return None, None151152def wit(lst, q):153 prev = -1154 for i, p in enumerate(lst):155 v = value(p, q)156 if v <= prev:157 return (lst[i - 1], p)158 prev = v159 return None160161def window(q):162 return max(l for l in range(2, 64) if all(closed(j) <= q for j in range(2, l + 1)))163164def ordwindow(tab, q):165 return max([l for l in sorted(tab) if tab[l] <= q] or [1])166167# VERBS168169def line(*a):170 print(*a, flush=True)171172def lemma():173 els, _ = monoid(WIDE + 1)174 line("THE CARRY BOUND for F = {0,1}, the monoid M* of products of 0/1 polynomials")175 line("")176 line(" max coefficient over degree < L, against the exact cap binomial(L-1, floor((L-1)/2)) and the crude cap 2^(L-1)")177 cmax = [0] * (WIDE + 2)178 for p in els:179 for a in coeffs(p):180 if a > cmax[deg(p)]:181 cmax[deg(p)] = a182 run = 0183 for L in range(2, WIDE + 2):184 run = max(run, cmax[L - 1])185 line(" L %2d max coef %6d exact cap %6d crude cap %8d"186 % (L, run, math.comb(L - 1, (L - 1) // 2), 2 ** (L - 1)))187 line("")188 line(" q_set(L): least q with every P of degree < L below q^L, and the last P over the cut")189 line(" q_ord(L): least q with evaluation increasing on all of M* to degree L")190 tset, tord = {}, {}191 for L in range(3, WIDE + 1):192 qs, ps = qset(L, els)193 qo, w = qord(L, els)194 tset[L], tord[L] = qs, qo195 line(" L %2d q_set %2d q_ord %2d closed %2d %s" % (L, qs, qo, closed(L),196 "" if closed(L) == qs else "SPLIT"))197 line(" cut at q = %d: %s = %d over q^L = %d, %d digits base %d"198 % (qs - 1, show(ps), value(ps, qs - 1), (qs - 1) ** L,199 len(digits(value(ps, qs - 1), qs - 1)), qs - 1))200 line(" order at q = %d: %s = %d before %s = %d"201 % (qo - 1, show(w[1]), value(w[1], qo - 1), show(w[0]), value(w[0], qo - 1)))202 line("")203 mx, sums, _, _ = ladder(WIDE + 2)204 line(" the pushforward term by term, nu_F(n) against the sum of nu*(P) over P(q) = n")205 for q, L in [(2, 15), (3, 10), (4, 8), (5, 7)]:206 bad, coll, dis = push(q, L)207 line(" q %2d to n <= q^L = %10d mismatches %d colliding values %d distinct values %d"208 % (q, q ** L, bad, coll, dis))209 line("")210 line(" the prediction against nu_F, sibling generator lab/py/mrly-pairing verb inverse")211 line(" base-free ladder: M(q^L) %s" % sums)212 line(" base-free ladder: maxima %s" % mx)213 for q, L in BASES:214 nu = pairing.dirichlet_inverse(q, 0b11, L)215 sig, run = pairing.ladder_stats(nu, q, L)216 ds = next((l + 1 for l in range(len(sig)) if sig[l] != sums[l]), None)217 dm = next((l + 1 for l in range(len(run)) if run[l] != mx[l]), None)218 line(" q %2d to level %2d carry-free window %2d order window %2d sum departs %s maxima depart %s"219 % (q, L, window(q), ordwindow(tord, q), ds, dm))220 line(" M(q^L) %s" % sig)221 line(" maxima %s" % run)222 del nu223224def push(q, L):225 els, nu = nustar(L + 1)226 g = {}227 for p in els:228 if nu[p] == 0:229 continue230 v = value(p, q)231 if v <= q ** L:232 g[v] = g.get(v, 0) + nu[p]233 nf = pairing.dirichlet_inverse(q, 0b11, L)234 bad = 0235 coll = 0236 seen = {}237 for p in els:238 v = value(p, q)239 if v <= q ** L:240 seen[v] = seen.get(v, 0) + 1241 coll += seen[v] > 1242 for n in range(1, q ** L + 1):243 if int(nf[n]) != g.get(n, 0):244 bad += 1245 return bad, coll, len(seen)246247def grep(seq):248 path = os.environ.get("OEIS_STRIPPED")249 if not path or not os.path.exists(path):250 return "no local dump named by OEIS_STRIPPED"251 key = "," + ",".join(str(v) for v in seq) + ","252 hits = []253 with open(path, encoding="utf-8", errors="ignore") as f:254 for row in f:255 if key in row:256 hits.append(row.split()[0])257 return ("absent" if not hits else ", ".join(hits[:8]))258259def sequence():260 mx, sums, cens, top = ladder(DEEP)261 line("THE BASE-FREE LADDER of nu*, levels 1..%d" % DEEP)262 line("")263 line(" level running max M(q^L) monoid census by degree")264 for l in range(1, DEEP + 1):265 line(" %2d %8d %6d %8d" % (l, mx[l - 1], sums[l - 1], cens[l - 1]))266 line("")267 line(" running maxima %s" % mx)268 line(" M(q^L) %s" % sums)269 line(" monoid census %s" % cens)270 line(" partial census %s" % [sum(cens[:l]) for l in range(1, DEEP + 1)])271 line("")272 line(" max coefficient at degree %d is %d, under the packing width 2^%d and the crude cap 2^%d"273 % (DEEP - 1, top, B, DEEP - 1))274 line("")275 line(" OEIS running maxima %s" % grep(mx[:12]))276 line(" OEIS monoid census %s" % grep(cens[:12]))277 line(" OEIS partial census %s" % grep([sum(cens[:l]) for l in range(1, 13)]))278279def exponent():280 from math import log281 mx, _, _, _ = ladder(DEEP)282 line("THE BASE-FREE MERTENS EXPONENT against the design mass 2^L")283 line("")284 line(" L max 2^L max/2^L log2(max)/L step log2 step")285 for l in range(1, DEEP + 1):286 st = mx[l - 1] / mx[l - 2] if l > 1 and mx[l - 2] else 0.0287 line(" %2d %8d %8d %9.6f %9.6f %6.4f %9.6f"288 % (l, mx[l - 1], 2 ** l, mx[l - 1] / 2 ** l,289 log(mx[l - 1], 2) / l if mx[l - 1] else 0.0,290 st, log(st, 2) if st else 0.0))291 line("")292 for r in (0, 1):293 cl = [l for l in range(1, DEEP + 1) if l % 2 == r]294 line(" L = %d mod 2 max/2^L %s" % (r, ["%.6f" % (mx[l - 1] / 2 ** l) for l in cl]))295 line("")296 for w in (4, 6, 8):297 g = (mx[DEEP - 1] / mx[DEEP - 1 - w]) ** (1.0 / w)298 line(" geometric mean step over the last %2d levels %.6f log2 %.6f" % (w, g, log(g, 2)))299 line("")300 line(" the design mass rate is 2 and q^(Re rho) at base 3 {0,1} is %.6f, lab/py/mrly-pairing verb box" % RHO3)301302VERBS = {"lemma": lemma, "sequence": sequence, "exponent": exponent}303304if __name__ == "__main__":305 VERBS[sys.argv[1]]()