verify.py
14.2 kB · python · 433 lines
1# AVATARS23import itertools4import math5from decimal import Decimal, getcontext6from fractions import Fraction78getcontext().prec = 40910GAMMA = Decimal("0.5772156649015328606065120900824024310422")111213def die(what, got, want):14 raise AssertionError("%s: got %r, want %r" % (what, got, want))151617def check(what, got, want):18 if got != want:19 die(what, got, want)202122def factor(m):23 f = {}24 d = 225 while d * d <= m:26 while m % d == 0:27 f[d] = f.get(d, 0) + 128 m //= d29 d += 130 if m > 1:31 f[m] = f.get(m, 0) + 132 return f333435def esym(b):36 poly = [1]37 for v in b:38 new = [0] * (len(poly) + 1)39 for i, c in enumerate(poly):40 new[i] += c41 new[i + 1] += c * v42 poly = new43 return poly444546def corners(D):47 return list(itertools.product((0, 1), repeat=D))484950def cell_counts(D, n):51 s = 2 * n + 152 counts = {c: 0 for c in corners(D)}53 for cell in itertools.product(range(1, s + 1), repeat=D):54 counts[tuple(1 if v % 2 == 0 else 0 for v in cell)] += 155 return counts565758def signature(F, D):59 f = [0] * (D + 1)60 for c in F:61 f[sum(c)] += 162 return tuple(f)636465def sig_fill(f, D, n):66 return sum(f[w] * (n + 1) ** (D - w) * n ** w for w in range(D + 1))676869def divisor_power(x, n):70 return math.prod(a * n + 1 for a in factor(x).values())717273def partitions(m, top=None):74 if top is None:75 top = m76 if m == 0:77 return 178 return sum(partitions(m - k, k) for k in range(1, min(m, top) + 1))798081def multisets(D, total):82 if D == 0:83 yield ()84 return85 for v in range(total + 1):86 for rest in multisets(D - 1, total - v):87 yield (v,) + rest888990def interpolate(pts):91 k = len(pts)92 coeffs = [Fraction(0)] * k93 for i in range(k):94 xi, yi = pts[i]95 num = [Fraction(1)]96 den = Fraction(1)97 for j in range(k):98 if i == j:99 continue100 xj = pts[j][0]101 den *= xi - xj102 new = [Fraction(0)] * (len(num) + 1)103 for t, c in enumerate(num):104 new[t + 1] += c105 new[t] += c * (-xj)106 num = new107 scale = Fraction(yi) / den108 for t, c in enumerate(num):109 coeffs[t] += c * scale110 while len(coeffs) > 1 and coeffs[-1] == 0:111 coeffs.pop()112 return coeffs113114115def slopes(coeffs):116 if len(coeffs) == 1:117 return [] if coeffs[0] == 1 else None118 if coeffs[0] != 1:119 return None120 lead = coeffs[-1]121 if lead.denominator != 1 or lead <= 0:122 return None123 lead = int(lead)124 deg = len(coeffs) - 1125 divs = [d for d in range(1, lead + 1) if lead % d == 0]126 for combo in itertools.combinations_with_replacement(divs, deg):127 if math.prod(combo) != lead:128 continue129 poly = [Fraction(1)]130 for a in combo:131 new = [Fraction(0)] * (len(poly) + 1)132 for t, c in enumerate(poly):133 new[t] += c134 new[t + 1] += c * a135 poly = new136 if poly == coeffs:137 return sorted(combo, reverse=True)138 return None139140141PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23]142143144def minimal_avatar(sl):145 return math.prod(PRIMES[i] ** a for i, a in enumerate(sl))146147148def observables(F, n):149 s = 2 * n + 1150 odd = [i for i in range(1, s + 1) if i % 2 == 1]151 even = [i for i in range(1, s + 1) if i % 2 == 0]152 cells = set()153 for c in F:154 axes = [even if c[k] else odd for k in range(3)]155 for cell in itertools.product(*axes):156 cells.add(cell)157 fill = len(cells)158 verts = set()159 edges = set()160 faces = {}161 for (x, y, z) in cells:162 for dx in (0, 1):163 for dy in (0, 1):164 for dz in (0, 1):165 verts.add((x - dx, y - dy, z - dz))166 for a in (0, 1):167 for b in (0, 1):168 edges.add((0, x - 1, y - a, z - b))169 edges.add((1, x - a, y - 1, z - b))170 edges.add((2, x - a, y - b, z - 1))171 for key in ((0, x, y, z), (0, x - 1, y, z), (1, x, y, z), (1, x, y - 1, z), (2, x, y, z), (2, x, y, z - 1)):172 faces[key] = faces.get(key, 0) + 1173 parent = {c: c for c in cells}174175 def find(a):176 while parent[a] != a:177 parent[a] = parent[parent[a]]178 a = parent[a]179 return a180181 for (x, y, z) in cells:182 for nb in ((x + 1, y, z), (x, y + 1, z), (x, y, z + 1)):183 if nb in cells:184 ra, rb = find((x, y, z)), find(nb)185 if ra != rb:186 parent[ra] = rb187 comps = len({find(c) for c in cells}) if cells else 0188 graph_edges = sum(1 for v in faces.values() if v == 2)189 return {190 "fill": fill,191 "voids": s ** 3 - fill,192 "surface": sum(1 for v in faces.values() if v == 1),193 "vertices": len(verts),194 "edges": len(edges),195 "faces": len(faces),196 "euler": len(verts) - len(edges) + len(faces) - fill,197 "components": comps,198 "cycle_rank": graph_edges - fill + comps,199 }200201202def orbit_reps():203 cs = corners(3)204 index = {c: i for i, c in enumerate(cs)}205 maps = []206 for perm in itertools.permutations(range(3)):207 for flip in itertools.product((0, 1), repeat=3):208 maps.append([index[tuple(cs[i][perm[k]] ^ flip[k] for k in range(3))] for i in range(8)])209 seen = set()210 reps = []211 for code in range(256):212 if code in seen:213 continue214 orbit = set()215 for m in maps:216 v = 0217 for i in range(8):218 if code >> i & 1:219 v |= 1 << m[i]220 orbit.add(v)221 seen |= orbit222 reps.append((code, len(orbit)))223 return reps224225226def code_design(code):227 cs = corners(3)228 return [cs[i] for i in range(8) if code >> i & 1]229230231def check_fill_law():232 D = 3233 cs = corners(D)234 for n in range(7):235 counts = cell_counts(D, n)236 check("cell parity total at n=%d" % n, sum(counts.values()), (2 * n + 1) ** D)237 for mask in range(256):238 F = [cs[i] for i in range(8) if mask >> i & 1]239 f = signature(F, D)240 got = sum(counts[c] for c in F)241 want = sig_fill(f, D, n)242 if got != want:243 die("fill law mask=%d n=%d" % (mask, n), got, want)244 print("fill law: all 256 designs literal at side %d" % (2 * n + 1))245246247def fill_tables():248 tables = {}249 for D in (1, 2, 3):250 cs = corners(D)251 counts = [cell_counts(D, n) for n in range(7)]252 table = {}253 for mask in range(1 << (1 << D)):254 F = [cs[i] for i in range(1 << D) if mask >> i & 1]255 key = tuple(sum(counts[n][c] for c in F) for n in range(7))256 table.setdefault(key, []).append(mask)257 tables[D] = table258 return tables259260261def check_criterion(tables):262 agreed = 0263 for x in range(2, 3001):264 fac = factor(x)265 D = len(fac)266 a = sorted(fac.values(), reverse=True)267 b = [v - 1 for v in a]268 e = esym(b)269 caps = [math.comb(D, w) for w in range(D + 1)]270 realizable = all(0 <= e[w] <= caps[w] for w in range(D + 1))271 inequality = sum(a) <= 2 * D272 check("criterion agreement at x=%d" % x, realizable, inequality)273 if D <= 3:274 key = tuple(divisor_power(x, n) for n in range(7))275 found = len(tables[D].get(key, []))276 want = math.prod(math.comb(caps[w], e[w]) for w in range(D + 1)) if realizable else 0277 check("design count at x=%d" % x, found, want)278 agreed += 1279 print("criterion: %d integers 2..3000, realizability and Omega<=2omega agree everywhere" % agreed)280281282def check_counterexamples(tables):283 for x, D in ((8, 1), (72, 2)):284 key = tuple(divisor_power(x, n) for n in range(7))285 check("no avatar for x=%d" % x, tables[D].get(key, []), [])286 check("distinct fill polynomials in dimension 2", len(tables[2]), 12)287 print("counterexamples: d(8^n) absent in dimension 1, d(72^n) absent among the 12 fills of dimension 2")288289290def check_census():291 want = [1, 2, 4, 7, 12, 19, 30, 45, 67]292 for D in range(9):293 sigs = {tuple(esym(sorted(b, reverse=True))) for b in multisets(D, D)}294 got = len(sigs)295 target = sum(partitions(m) for m in range(D + 1))296 check("census at D=%d" % D, got, target)297 check("census value at D=%d" % D, got, want[D])298 print("census: D=%d qualifying polynomials %d" % (D, got))299300301def check_sponge():302 F = [c for c in corners(3) if sum(c) <= 1]303 for n in range(11):304 s = 2 * n + 1305 removed = 0306 kept = 0307 for cell in itertools.product(range(1, s + 1), repeat=3):308 if sum(1 for v in cell if v % 2 == 0) >= 2:309 removed += 1310 else:311 kept += 1312 check("void at n=%d" % n, removed, n * n * (4 * n + 3))313 check("fill at n=%d" % n, kept, divisor_power(240, n))314 check("complement at n=%d" % n, removed + kept, s ** 3)315 check("signature of the Menger design", signature(F, 3), (1, 3, 0, 0))316 check("a(1) is the centre and six face centres", 1 * 1 * (4 * 1 + 3), 7)317 fac = factor(240)318 check("240 sits on the boundary", (sum(fac.values()), 2 * len(fac)), (6, 6))319 fac = factor(480)320 check("480 is past the boundary", sum(fac.values()) > 2 * len(fac), True)321 print("sponge: fill d(240^n) and void n^2(4n+3) literal for n=0..10, 240 on the boundary, 480 past it")322323324def check_robin():325 ladders = {30: {2: 1, 3: 1, 5: 1}, 60: {2: 2, 3: 1, 5: 1}, 120: {2: 3, 3: 1, 5: 1},326 180: {2: 2, 3: 2, 5: 1}, 240: {2: 4, 3: 1, 5: 1}, 360: {2: 3, 3: 2, 5: 1},327 900: {2: 2, 3: 2, 5: 2}}328 eg = GAMMA.exp()329 check("e^gamma to twelve places", eg.quantize(Decimal("1.000000000000")), Decimal("1.781072417990"))330 threshold = (Decimal(15) / (4 * eg)).exp().exp()331 check("threshold below 5040", threshold < 5040, True)332 check("threshold to two places", threshold.quantize(Decimal("1.00")), Decimal("3681.17"))333 total = 0334 above = 0335 seen = set()336 seen_above = set()337 best = None338 for x, fac in ladders.items():339 for n in range(1, 21):340 total += 1341 N = x ** n342 seen.add(N)343 s = math.prod((p ** (e * n + 1) - 1) // (p - 1) for p, e in fac.items())344 check("abundancy below 15/4 at x=%d n=%d" % (x, n), Decimal(s) / Decimal(N) < Decimal(15) / 4, True)345 if N <= 5040:346 continue347 above += 1348 seen_above.add(N)349 ratio = Decimal(s) / (Decimal(N) * Decimal(N).ln().ln())350 check("Robin at N=%d" % N, ratio < eg, True)351 if best is None or ratio > best[0]:352 best = (ratio, N)353 check("powers in the domain", total, 140)354 check("distinct integers in the domain", len(seen), 130)355 check("powers above 5040", above, 131)356 check("distinct integers above 5040", len(seen_above), 122)357 check("argmax of the Robin ratio", best[1], 14400)358 check("max Robin ratio", best[0].quantize(Decimal("1.000000000000")), Decimal("1.573259905933"))359 print("robin: 140 powers, 131 above 5040, all satisfy Robin, max %s at N=14400" % best[0].quantize(Decimal("1.000000000000")))360361362def check_ca():363 ca = [2, 6, 12, 60, 120, 360, 2520, 5040, 55440, 720720, 1441440, 4324320, 21621600]364 verdicts = []365 for m in ca:366 fac = factor(m)367 verdicts.append(sum(fac.values()) <= 2 * len(fac))368 check("first twelve colossally abundant numbers are avatars", verdicts[:12], [True] * 12)369 check("the thirteenth is not", verdicts[12], False)370 fac = factor(21621600)371 check("factorisation of 21621600", fac, {2: 5, 3: 3, 5: 2, 7: 1, 11: 1, 13: 1})372 check("Omega and omega of 21621600", (sum(fac.values()), len(fac)), (13, 6))373 print("colossally abundant: first 12 of A004490 are avatars, 21621600 has Omega=13 > 12=2omega")374375376def check_scan():377 reps = orbit_reps()378 check("orbit count", len(reps), 22)379 check("orbit sizes sum to 256", sum(r[1] for r in reps), 256)380 names = ["fill", "voids", "surface", "vertices", "edges", "faces", "euler", "components", "cycle_rank"]381 nonfill = []382 constants = []383 fills = []384 late = []385 for code, _ in reps:386 F = code_design(code)387 vals = [observables(F, n) for n in range(11)]388 for obs in names:389 seq = [v[obs] for v in vals]390 coeffs = interpolate([(n, seq[n]) for n in range(11)])391 if len(coeffs) > 4:392 tail = interpolate([(n, seq[n]) for n in range(1, 5)])393 if not all(sum(c * n ** t for t, c in enumerate(tail)) == seq[n] for n in range(1, 11)):394 die("law shape code=%d obs=%s" % (code, obs), len(coeffs) - 1, "degree at most 3")395 check("late law is zero at n=0 code=%d obs=%s" % (code, obs), seq[0], 0)396 late.append((code, obs))397 sl = slopes(coeffs)398 if sl is None:399 continue400 if sl == []:401 constants.append((code, obs))402 elif obs == "fill":403 fills.append((code, minimal_avatar(sl)))404 else:405 nonfill.append((code, obs, minimal_avatar(sl)))406 print("scan: code %d done" % code)407 check("late-starting laws", sorted(late), [(30, "components"), (30, "cycle_rank"), (126, "components"), (126, "cycle_rank")])408 check("non-fill divisor avatars", nonfill, [409 (0, "voids", 900), (1, "euler", 30), (1, "components", 30),410 (3, "euler", 6), (3, "components", 6), (7, "components", 2),411 (15, "euler", 2), (15, "components", 2)])412 check("constant identities", sorted(constants), sorted([413 (23, "components"), (27, "components"), (31, "components"), (61, "components"),414 (63, "components"), (111, "components"), (127, "components"), (255, "euler"), (255, "components")]))415 check("fill avatars", fills, [(1, 30), (3, 60), (7, 120), (15, 180), (23, 240), (27, 180), (63, 360), (255, 900)])416 print("scan: 22 representatives, 9 observables, n=0..10, exactly 8 non-fill rows and 9 constants")417418419def main():420 check_fill_law()421 tables = fill_tables()422 check_criterion(tables)423 check_counterexamples(tables)424 check_census()425 check_sponge()426 check_robin()427 check_ca()428 check_scan()429 print("all green")430431432if __name__ == "__main__":433 main()