verify.py
24.5 kB · python · 685 lines
1import time2from decimal import Decimal, getcontext3from fractions import Fraction4from itertools import permutations, product5from math import comb, factorial, gcd, prod67# POLYNOMIALS89def pmul(a, b):10 out = [0] * (len(a) + len(b) - 1)11 for i, x in enumerate(a):12 if x:13 for j, y in enumerate(b):14 if y:15 out[i + j] += x * y16 return out1718def ppow(a, e):19 r = [1]20 for _ in range(e):21 r = pmul(r, a)22 return r2324def padd(a, b):25 out = [0] * max(len(a), len(b))26 for i, x in enumerate(a):27 out[i] += x28 for i, x in enumerate(b):29 out[i] += x30 return trim(out)3132def pscale(a, s):33 return trim([s * x for x in a])3435def trim(a):36 while len(a) > 1 and a[-1] == 0:37 a.pop()38 return a3940def peval(a, x):41 v = 042 for c in reversed(a):43 v = v * x + c44 return v4546def psub(a, b):47 return padd(a, pscale(b, -1))4849def compose_affine(a, u, v):50 out = [0]51 for c in reversed(a):52 out = padd(pmul(out, [v, u]), [c])53 return out5455# DESIGNS5657def corners(D):58 return list(product((0, 1), repeat=D))5960def code(F, D):61 return sum(1 << sum(c[j] << j for j in range(D)) for c in F)6263def signature(F, D):64 sig = [0] * (D + 1)65 for c in F:66 sig[sum(c)] += 167 return tuple(sig)6869def designs(D):70 cs = corners(D)71 for mask in range(1 << len(cs)):72 yield [c for i, c in enumerate(cs) if mask >> i & 1]7374def signatures(D):75 ranges = [range(comb(D, w) + 1) for w in range(D + 1)]76 return [tuple(s) for s in product(*ranges)]7778def fillpoly(sig):79 D = len(sig) - 180 out = [0]81 for w, f in enumerate(sig):82 if f:83 out = padd(out, pscale(pmul(ppow([0, 1], D - w), ppow([-1, 1], w)), f))84 return out8586def grid_fill(F, D, k):87 n = 2 * k - 188 want = set(F)89 total = 090 for cell in product(range(n), repeat=D):91 if tuple(x & 1 for x in cell) in want:92 total += 193 return total9495# CLASSICAL FAMILIES9697def polygonal(m, k):98 return ((m - 2) * k * k - (m - 4) * k) // 299100def centered(m, k):101 return m * k * (k - 1) // 2 + 1102103def centered_hex(m):104 return 3 * m * m - 3 * m + 1105106def is_prime(v):107 if v < 2:108 return False109 d = 2110 while d * d <= v:111 if v % d == 0:112 return False113 d += 1114 return True115116# RECORDS117118RECORDS = {119 "A000290": (0, [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121]),120 "A000384": (0, [0, 1, 6, 15, 28, 45, 66, 91, 120, 153, 190, 231]),121 "A000567": (0, [0, 1, 8, 21, 40, 65, 96, 133, 176, 225, 280, 341]),122 "A001844": (0, [1, 5, 13, 25, 41, 61, 85, 113, 145, 181, 221, 265]),123 "A003215": (0, [1, 7, 19, 37, 61, 91, 127, 169, 217, 271, 331, 397]),124 "A016754": (0, [1, 9, 25, 49, 81, 121, 169, 225, 289, 361, 441, 529]),125 "A000578": (0, [0, 1, 8, 27, 64, 125, 216, 343, 512, 729, 1000, 1331]),126 "A103532": (0, [1, 20, 81, 208, 425, 756, 1225, 1856, 2673, 3700, 4961, 6480]),127 "A395241": (0, [0, 7, 44, 135, 304, 575, 972, 1519, 2240, 3159, 4300, 5687]),128 "A005898": (0, [1, 9, 35, 91, 189, 341, 559, 855, 1241, 1729, 2331, 3059]),129 "A016755": (0, [1, 27, 125, 343, 729, 1331, 2197, 3375, 4913, 6859, 9261, 12167]),130 "A001018": (0, [1, 8, 64, 512, 4096, 32768, 262144, 2097152, 16777216]),131 "A016185": (0, [0, 1, 17, 217, 2465, 26281, 269297, 2685817, 26269505]),132 "A381517": (0, [4, 16, 80, 496, 3536, 26992, 212048, 1684720, 13442768]),133 "A009964": (0, [1, 20, 400, 8000, 160000, 3200000, 64000000, 1280000000]),134 "A332705": (0, [6, 72, 1056, 18048, 336384, 6531072, 129048576, 2568388608]),135 "A000616": (-1, [1, 2, 3, 6, 22, 402, 1228158, 400507806843728]),136 "A129824": (0, [2, 4, 12, 64, 700, 17424, 1053696, 160579584, 62856336636]),137 "A396934": (0, [0, 2, 4, 12, 34, 122, 362, 1130, 3406, 10506, 31550, 95260]),138 "A398348": (1, [2, 22, 111618, 6005363762644688,139 7089215977519836239803174210135872]),140 "A154105": (0, [7, 37, 91, 169, 271, 397, 547, 721, 919, 1141, 1387, 1657]),141 "A299916": (0, [1, 6, 42, 306, 2250, 16578, 122202, 900882, 6641514, 48963042]),142 "A056040": (0, [1, 1, 2, 6, 6, 30, 20, 140, 70, 630, 252, 2772, 924, 12012, 3432,143 51480, 12870]),144}145146def record(name, index):147 offset, data = RECORDS[name]148 slot = index - offset149 assert 0 <= slot < len(data), "%s index %d outside the stored terms" % (name, index)150 return data[slot]151152# THE SIX ROWS OF THE PLANE153154PLANE = [155 (1, (1, 0, 0), "A000290", 0, ("polygonal", 4)),156 (3, (1, 1, 0), "A000384", 0, ("polygonal", 6)),157 (7, (1, 2, 0), "A000567", 0, ("polygonal", 8)),158 (9, (1, 0, 1), "A001844", -1, ("centered", 4)),159 (11, (1, 1, 1), "A003215", -1, ("centered", 6)),160 (15, (1, 2, 1), "A016754", -1, ("centered", 8)),161]162163SOLID = [164 (1, (1, 0, 0, 0), "A000578", 0),165 (23, (1, 3, 0, 0), "A103532", -1),166 (232, (0, 0, 3, 1), "A395241", -1),167 (129, (1, 0, 0, 1), "A005898", -1),168 (255, (1, 3, 3, 1), "A016755", -1),169]170171# CHECKS172173def check_fill_law():174 for D in (1, 2, 3):175 for F in designs(D):176 sig = signature(F, D)177 poly = fillpoly(sig)178 assert len(poly) - 1 <= D, (D, sig)179 if F:180 assert poly[D] == len(F), (D, sig, poly)181 else:182 assert poly == [0], (D, sig)183 for k in range(1, 8):184 want = grid_fill(F, D, k)185 got = peval(poly, k)186 closed = sum(k ** (D - sum(c)) * (k - 1) ** sum(c) for c in F)187 assert want == got == closed, (D, code(F, D), k, want, got, closed)188 return "D = 1,2,3, all 4/16/256 designs, k = 1..7"189190def check_endpoints():191 for D in range(1, 5):192 for sig in signatures(D):193 poly = fillpoly(sig)194 assert peval(poly, 1) == sig[0], (D, sig)195 assert peval(poly, 0) == (-1) ** D * sig[D], (D, sig)196 rev = fillpoly(tuple(reversed(sig)))197 mirror = compose_affine(poly, -1, 1)198 assert mirror == pscale(rev, (-1) ** D), (D, sig)199 diff = poly200 for _ in range(D):201 diff = psub(compose_affine(diff, 1, 1), diff)202 assert diff == [factorial(D) * sum(sig)] or (sum(sig) == 0 and diff == [0]), (D, sig)203 return "D = 1..4, every weight signature"204205def check_plane():206 seen = {}207 for F in designs(2):208 seen.setdefault(signature(F, 2), []).append(code(F, 2))209 assert len(seen) == 12, len(seen)210 for sig, codes in seen.items():211 poly = fillpoly(sig)212 p = sum(sig)213 f0, f1, f2 = sig214 for k in range(0, 31):215 got = peval(poly, k)216 if f0 == 1 and f2 == 0:217 assert got == polygonal(2 * p + 2, k), (sig, k)218 elif f0 == 1 and f2 == 1:219 assert got == centered(2 * p, k), (sig, k)220 elif f2 == 1:221 assert got == polygonal(2 * p + 2, 1 - k), (sig, k)222 assert got == peval(fillpoly(tuple(reversed(sig))), 1 - k), (sig, k)223 else:224 assert got == f1 * k * (k - 1), (sig, k)225 assert got == 2 * f1 * ((k - 1) * k // 2), (sig, k)226 assert got - peval(poly, k - 1) == 2 * p * (k - 1) + f0 - f2, (sig, k)227 for c, sig, name, shift, family in PLANE:228 F = [x for x in corners(2) if (c >> (x[0] + 2 * x[1])) & 1]229 assert signature(F, 2) == sig, (c, signature(F, 2))230 poly = fillpoly(sig)231 kind, m = family232 for k in range(2, 10):233 got = peval(poly, k)234 assert got == grid_fill(F, 2, k), (c, k)235 assert got == record(name, k + shift), (c, name, k, got)236 if kind == "polygonal":237 assert got == polygonal(m, k) and m == 2 * sum(sig) + 2, (c, k)238 else:239 assert got == centered(m, k) and m == 2 * sum(sig), (c, k)240 return "all 12 plane signatures at k = 0..30, the six records at k = 2..9"241242def check_solid():243 for c, sig, name, shift in SOLID:244 F = [x for x in corners(3) if (c >> (x[0] + 2 * x[1] + 4 * x[2])) & 1]245 assert signature(F, 3) == sig, (c, signature(F, 3))246 poly = fillpoly(sig)247 for k in range(2, 10):248 got = peval(poly, k)249 assert got == grid_fill(F, 3, k), (c, k)250 assert got == record(name, k + shift), (c, name, k, got, record(name, k + shift))251 solid = fillpoly((1, 3, 3, 1))252 sponge = fillpoly((1, 3, 0, 0))253 void = fillpoly((0, 0, 3, 1))254 assert padd(sponge, void) == solid, (sponge, void, solid)255 return "D = 3 records at k = 2..9, complement identity as polynomials"256257def check_census():258 for D in range(1, 5):259 polys = set()260 for F in designs(D):261 polys.add(tuple(fillpoly(signature(F, D))))262 closed = prod(1 + comb(D, w) for w in range(D + 1))263 assert len(polys) == closed == record("A129824", D), (D, len(polys), closed)264 for D in range(0, 9):265 assert prod(1 + comb(D, w) for w in range(D + 1)) == record("A129824", D), D266 return "distinct fill polynomials enumerated at D = 1..4, closed form to D = 8"267268def burnside_cube(D):269 cs = corners(D)270 index = {c: i for i, c in enumerate(cs)}271 total = 0272 for perm in permutations(range(D)):273 for t in range(1 << D):274 img = [index[tuple(c[perm[i]] ^ (t >> i & 1) for i in range(D))] for c in cs]275 total += 1 << cycles(img)276 return total // ((1 << D) * factorial(D))277278def cycles(img):279 seen = [False] * len(img)280 count = 0281 for s in range(len(img)):282 if not seen[s]:283 count += 1284 j = s285 while not seen[j]:286 seen[j] = True287 j = img[j]288 return count289290def orbit_count(D):291 cs = corners(D)292 index = {c: i for i, c in enumerate(cs)}293 maps = []294 for perm in permutations(range(D)):295 for t in range(1 << D):296 maps.append([index[tuple(c[perm[i]] ^ (t >> i & 1) for i in range(D))] for c in cs])297 reps = set()298 for mask in range(1 << len(cs)):299 best = mask300 for img in maps:301 moved = 0302 for i in range(len(cs)):303 if mask >> i & 1:304 moved |= 1 << img[i]305 best = min(best, moved)306 reps.add(best)307 return len(reps)308309def check_shapes():310 for D in range(1, 7):311 got = burnside_cube(D)312 assert got == record("A000616", D), (D, got, record("A000616", D))313 for D in (1, 2, 3):314 assert orbit_count(D) == record("A000616", D), D315 seq = [prod(1 + comb(D, w) for w in range(D + 1)) for D in range(0, 7)]316 shapes = [record("A000616", D) for D in range(0, 7)]317 for D in range(1, 5):318 assert seq[D] > shapes[D], (D, seq[D], shapes[D])319 for D in (5, 6):320 assert seq[D] < shapes[D], (D, seq[D], shapes[D])321 assert shapes[6] // seq[6] > 380000000, shapes[6] // seq[6]322 return "Burnside D = 1..6, orbit walk D = 1..3, crossover at D = 5"323324def burnside_torus3(n):325 cells = [(x, y, z) for x in range(n) for y in range(n) for z in range(n)]326 index = {c: i for i, c in enumerate(cells)}327 line = [(s, b) for s in (1, -1) for b in range(n)]328 total = 0329 for perm in permutations(range(3)):330 for m in product(line, repeat=3):331 img = []332 for c in cells:333 p = (c[perm[0]], c[perm[1]], c[perm[2]])334 img.append(index[tuple((m[i][0] * p[i] + m[i][1]) % n for i in range(3))])335 total += 1 << cycles(img)336 return total // (48 * n ** 3)337338def check_torus():339 for n in range(1, 6):340 got = burnside_torus3(n)341 assert got == record("A398348", n), (n, got)342 assert burnside_torus3(3) == 111618343 return "A398348 recomputed at n = 1..5, group order 48 n^3"344345def tile(F, D, q):346 keep = set(F)347 return [c for c in product(range(q), repeat=D) if tuple(x & 1 for x in c) in keep]348349def fractal(F, D, q, L):350 cells = {tuple([0] * D)}351 base = tile(F, D, q)352 for _ in range(L):353 cells = {tuple(c[i] * q + b[i] for i in range(D)) for c in cells for b in base}354 return cells355356def surface(cells, D):357 total = 0358 for c in cells:359 for i in range(D):360 for step in (-1, 1):361 nb = list(c)362 nb[i] += step363 if tuple(nb) not in cells:364 total += 1365 return total366367def check_level():368 carpet2 = [c for c in corners(2) if sum(c) <= 1]369 carpet3 = [c for c in corners(3) if sum(c) <= 1]370 void2 = [c for c in corners(2) if sum(c) in (0, 2)]371 solid3 = corners(3)372 for L in range(1, 5):373 cells = fractal(carpet2, 2, 3, L)374 assert len(cells) == 8 ** L == record("A001018", L), L375 assert 9 ** L - len(cells) == record("A016185", L), L376 assert surface(cells, 2) == record("A381517", L), (L, surface(cells, 2))377 assert surface(cells, 2) == (4 * 8 ** L + 16 * 3 ** L) // 5, L378 for L in range(1, 4):379 cells = fractal(carpet3, 3, 3, L)380 assert len(cells) == 20 ** L == record("A009964", L), L381 assert surface(cells, 3) == record("A332705", L), (L, surface(cells, 3))382 assert surface(cells, 3) == 2 * 20 ** L + 4 * 8 ** L, L383 for L in range(1, 4):384 assert len(fractal(void2, 2, 3, L)) == 5 ** L, L385 assert len(fractal(solid3, 3, 3, L)) == 27 ** L, L386 assert surface(fractal(solid3, 3, 3, L), 3) == 6 * 9 ** L, L387 for L in range(3, 5):388 assert record("A381517", L) == 11 * record("A381517", L - 1) - 24 * record("A381517", L - 2), L389 for L in range(3, 8):390 assert record("A332705", L) == 28 * record("A332705", L - 1) - 160 * record("A332705", L - 2), L391 return "carpet to L = 4, sponge to L = 3, cells and surface counted face by face"392393def check_gasket():394 for n in range(0, 12):395 total = 0396 for i in range(1 << n):397 free = ((1 << n) - 1) ^ i398 j = free399 while True:400 if gcd(i, j) == 1:401 total += 1402 if j == 0:403 break404 j = (j - 1) & free405 assert total == record("A396934", n), (n, total)406 pairs = 0407 for i in range(1 << n):408 free = ((1 << n) - 1) ^ i409 pairs += 1 << bin(free).count("1")410 assert pairs == 3 ** n, n411 return "A396934 counted pair by pair at n = 0..11, support 3^n"412413def check_mesh():414 tree = []415 for k in range(1, 21):416 R = 2 * k - 1417 pts = 0418 for x in range(-R, R + 1):419 for y in range(-R, R + 1):420 z = -x - y421 if abs(z) <= R and abs(x) <= R and abs(y) <= R:422 pts += 1423 assert pts == 12 * k * k - 6 * k + 1, (k, pts)424 assert pts == centered_hex(2 * k), k425 assert pts == 3 * R * R + 3 * R + 1, k426 assert pts % 3 == 1, k427 if k <= 12:428 assert pts == record("A154105", k - 1), k429 if is_prime(pts):430 tree.append(pts)431 assert tree == [7, 37, 271, 397, 547, 919, 1657, 1951, 2269, 4219], tree432 for m in range(1, 60):433 assert centered_hex(m) == m ** 3 - (m - 1) ** 3, m434 return "hexagon lattice points at k = 1..20, ten prime vertex counts"435436def check_pigeonhole():437 for a in range(1, 40):438 p = [0, 1 - a, a]439 c = [1, -a, a]440 assert peval(p, 0) == 0 and peval(p, 1) == 1, a441 assert peval(c, 0) == 1 and peval(c, 1) == 1, a442 for k in range(0, 20):443 assert peval(p, k) == polygonal(2 * a + 2, k), (a, k)444 assert peval(c, k) == centered(2 * a, k), (a, k)445 for a in range(1, 12):446 for b in range(-30, 31):447 for c0 in range(-3, 4):448 q = [c0, b, a]449 if peval(q, 0) == 0 and peval(q, 1) == 1:450 assert q == [0, 1 - a, a], q451 if peval(q, 0) == 1 and peval(q, 1) == 1:452 assert q == [1, -a, a], q453 return "the two normalisations pin a quadratic outright, leading coefficient 1..11"454455def slab_data(F, D, q):456 T = tile(F, D, q)457 seen = set(T)458 c = len(T)459 ls, Ws = [], []460 for a in range(D):461 ls.append(sum(1 for t in T if t[a] == 0))462 assert ls[a] == sum(1 for t in T if t[a] == q - 1), (F, a)463 Ws.append(sum(1 for t in T if tuple(t[i] + (i == a) for i in range(D)) in seen))464 return c, ls, Ws465466def occupancy(F, D, q, L):467 keep = set(F)468 side = q ** L469 grid = bytearray(side ** D)470 for cell in product(range(side), repeat=D):471 ok = True472 for j in range(L):473 if tuple((x // q ** j) % q & 1 for x in cell) not in keep:474 ok = False475 break476 if ok:477 idx = 0478 for x in cell:479 idx = idx * side + x480 grid[idx] = 1481 return grid, side482483def faces(grid, side, D):484 strides = [side ** (D - 1 - i) for i in range(D)]485 total = 0486 for idx in range(len(grid)):487 if not grid[idx]:488 continue489 rest = idx490 coord = []491 for stride in strides:492 coord.append(rest // stride)493 rest %= stride494 for i in range(D):495 for step in (-1, 1):496 v = coord[i] + step497 if v < 0 or v >= side or not grid[idx + step * strides[i]]:498 total += 1499 return total500501def check_surface():502 split = 0503 witness = None504 for D, top in ((2, 4), (3, 3)):505 for mask in range(1, 1 << (1 << D)):506 cs = corners(D)507 F = [c for i, c in enumerate(cs) if mask >> i & 1]508 c, ls, Ws = slab_data(F, D, 3)509 for a in range(D):510 assert ls[a] < c, (D, mask, ls, c)511 sur = []512 for L in range(top + 1):513 grid, side = occupancy(F, D, 3, L)514 sur.append(faces(grid, side, D))515 assert sur[0] == 2 * D, (D, mask, sur)516 for L in range(top):517 want = c * sur[L] - 2 * sum(Ws[a] * ls[a] ** L for a in range(D))518 assert sur[L + 1] == want, (D, mask, L, sur, want)519 coef = {}520 for a in range(D):521 if Ws[a]:522 coef[ls[a]] = coef.get(ls[a], 0) + Fraction(2 * Ws[a], c - ls[a])523 lead = Fraction(2 * D) - sum(coef.values())524 for L in range(top + 1):525 got = lead * c ** L + sum(b * v ** L for v, b in coef.items())526 assert got == sur[L], (D, mask, L, got, sur[L])527 if D == 3 and len(coef) > 1:528 split += 1529 if D == 2 and mask == 11:530 witness = (c, ls, Ws, sur)531 assert split == 141, split532 c, ls, Ws, sur = witness533 assert (c, ls, Ws) == (7, [3, 2], [2, 4]), witness534 grid, side = occupancy([x for x in corners(2) if (11 >> (x[0] + 2 * x[1])) & 1], 2, 3, 5)535 sur = sur + [faces(grid, side, 2)]536 assert sur == [4, 16, 84, 520, 3468, 23824], sur537 for L in range(3, 6):538 assert sur[L] == 12 * sur[L - 1] - 41 * sur[L - 2] + 42 * sur[L - 3], L539 det = sur[1] * sur[1] - sur[0] * sur[2]540 assert det != 0, det541 alpha = Fraction(sur[3] * sur[1] - sur[2] * sur[2], det)542 beta = Fraction(sur[2] * sur[0] - sur[1] * sur[1], det)543 assert alpha * sur[3] + beta * sur[2] != sur[4], (alpha, beta)544 return "all 15 plane and 255 solid designs, faces counted literally to L = 4 and L = 3"545546def menger_analog(D):547 return [v for v in product(range(3), repeat=D) if sum(1 for x in v if x == 1) <= 1]548549def digit_weights(D):550 weights = {}551 for v in menger_analog(D):552 weights[sum(v)] = weights.get(sum(v), 0) + 1553 return weights554555def ladder(D, L):556 weights = digit_weights(D)557 state = {0: 1}558 for _ in range(L):559 nxt = {}560 for c, count in state.items():561 for s, w in weights.items():562 m = s - D563 if (c + m) % 3 == 0:564 key = (c + m) // 3565 nxt[key] = nxt.get(key, 0) + count * w566 state = nxt567 return state.get(0, 0)568569def check_ladder():570 for L in range(0, 9):571 assert ladder(3, L) == record("A299916", L), (L, ladder(3, L))572 for L in range(2, 9):573 assert ladder(3, L) == 9 * ladder(3, L - 1) - 12 * ladder(3, L - 2), L574 four = [ladder(4, L) for L in range(0, 9)]575 assert four[:7] == [1, 6, 132, 1848, 29040, 441408, 6772128], four576 for L in range(2, 9):577 assert four[L] == 11 * four[L - 1] + 66 * four[L - 2], L578 assert four[2] * four[2] - four[1] * four[3] != 0579 for D in range(2, 11):580 cells = [v for v in menger_analog(D) if sum(v) == D]581 closed = comb(D, D // 2) if D % 2 == 0 else comb(D, (D - 1) // 2) * (D + 1) // 2582 assert len(cells) == closed == ladder(D, 1), (D, len(cells), closed)583 tree = [2, 6, 6, 30, 20, 140, 70, 630, 252]584 assert [ladder(D, 1) for D in range(2, 11)] == tree585 for D in range(1, 17):586 swing = factorial(D) // factorial(D // 2) ** 2587 assert swing == record("A056040", D), (D, swing)588 if D >= 2:589 closed = comb(D, D // 2) if D % 2 == 0 else comb(D, (D - 1) // 2) * (D + 1) // 2590 assert closed == swing, (D, closed, swing)591 assert 121 + 4 * 66 == 385592 getcontext().prec = 40593 root = (Decimal(11) + Decimal(385).sqrt()) / 2594 assert abs(root * root - 11 * root - 66) < Decimal("1e-30"), root595 assert str(root)[:11] == "15.31070843", root596 exponent = root.ln() / Decimal(3).ln()597 assert str(exponent)[:12] == "2.4836355003", exponent598 return "carry ladder at D = 3,4 to L = 8, level-one slice at D = 2..10 and A056040 to 16"599600# TABLES601602def show(poly):603 parts = []604 for e in range(len(poly) - 1, -1, -1):605 c = poly[e]606 if not c:607 continue608 term = "k^%d" % e if e > 1 else ("k" if e == 1 else "")609 head = "" if abs(c) == 1 and e else str(abs(c))610 parts.append(("- " if c < 0 else "+ ") + head + term)611 if not parts:612 return "0"613 body = " ".join(parts)614 return body[2:] if body.startswith("+ ") else "-" + body[2:]615616def tables():617 print("")618 print("the twelve fill sequences of the plane")619 print(" %-9s %-10s %-16s %s" % ("signature", "codes", "fill at n = 2k-1", "family"))620 rows = {}621 for F in designs(2):622 rows.setdefault(signature(F, 2), []).append(code(F, 2))623 for sig in sorted(rows):624 f0, f1, f2 = sig625 p = sum(sig)626 if p == 0:627 family = "empty"628 elif f0 == 1 and f2 == 0:629 family = "polygonal m = %d" % (2 * p + 2)630 elif f0 == 1 and f2 == 1:631 family = "centered m = %d" % (2 * p)632 elif sig == tuple(reversed(sig)):633 family = "self-mirror"634 else:635 family = "mirror of %s" % "".join(map(str, reversed(sig)))636 codes = ",".join(str(c) for c in sorted(rows[sig]))637 print(" %-9s %-10s %-16s %s" % ("".join(map(str, sig)), codes, show(fillpoly(sig)), family))638 print("")639 print("the two censuses")640 print(" %-3s %-22s %-12s %s" % ("D", "designs", "sequences", "shapes"))641 for D in range(0, 7):642 seq = prod(1 + comb(D, w) for w in range(D + 1))643 print(" %-3d %-22d %-12d %d" % (D, 1 << (1 << D), seq, record("A000616", D)))644 print("")645 print("the records this census reads")646 for name, shift, sig, D in [("A000290", 0, (1, 0, 0), 2), ("A000384", 0, (1, 1, 0), 2),647 ("A000567", 0, (1, 2, 0), 2), ("A001844", -1, (1, 0, 1), 2),648 ("A003215", -1, (1, 1, 1), 2), ("A016754", -1, (1, 2, 1), 2),649 ("A000578", 0, (1, 0, 0, 0), 3), ("A103532", -1, (1, 3, 0, 0), 3),650 ("A395241", -1, (0, 0, 3, 1), 3), ("A005898", -1, (1, 0, 0, 1), 3),651 ("A016755", -1, (1, 3, 3, 1), 3)]:652 terms = [peval(fillpoly(sig), k) for k in range(2, 8)]653 print(" %-8s shift %-3d %-24s %s" % (name, shift, show(fillpoly(sig)),654 ", ".join(map(str, terms))))655656# DOOR657658def main():659 t0 = time.time()660 checks = [661 ("fill law", check_fill_law),662 ("endpoint, mirror and difference laws", check_endpoints),663 ("the plane", check_plane),664 ("the solid", check_solid),665 ("sequence census", check_census),666 ("shape census", check_shapes),667 ("toroidal census", check_torus),668 ("level axis", check_level),669 ("the surface law in general", check_surface),670 ("gasket coprimality", check_gasket),671 ("slice mesh", check_mesh),672 ("the diagonal ladder", check_ladder),673 ("classical pigeonhole", check_pigeonhole),674 ]675 for name, fn in checks:676 t = time.time()677 domain = fn()678 print("%-38s PASS %-58s %5.1f s" % (name, domain, time.time() - t))679 print("total %.1f s" % (time.time() - t0))680 tables()681 print("")682 print("all green")683684if __name__ == "__main__":685 main()