fills.py
14.6 kB · python · 443 lines
1import math2import operator3import os4import time5from concurrent.futures import ProcessPoolExecutor6from fractions import Fraction7from itertools import permutations, product89from sympy import Poly10from sympy.abc import n as VAR1112KMAX = 1013DISC_DIMS = (2, 3, 4, 5, 6)14PARALLEL_FROM = 515CHUNKS = 9616LADDER_DEPTH = 4001718def corners(d):19 return [tuple((i >> (d - 1 - a)) & 1 for a in range(d)) for i in range(1 << d)]2021def block_tables(values, combine):22 tables = []23 for start in range(0, len(values), 8):24 block = values[start:start + 8]25 table = [0] * (1 << len(block))26 for mask in range(1, 1 << len(block)):27 low = mask & -mask28 table[mask] = combine(table[mask ^ low], block[low.bit_length() - 1])29 tables.append(table)30 return tables3132def block_apply(tables, code, combine):33 out = 034 for table in tables:35 out = combine(out, table[code & (len(table) - 1)])36 code >>= 837 return out3839def sum_tables(values):40 return block_tables(values, operator.add)4142def or_tables(mapping):43 return block_tables([1 << target for target in mapping], operator.or_)4445def block_sum(tables, code):46 return block_apply(tables, code, operator.add)4748def block_or(tables, code):49 return block_apply(tables, code, operator.or_)5051def fill_tables(d, k):52 return sum_tables([k ** (d - sum(c)) * (k - 1) ** sum(c) for c in corners(d)])5354def polymul(a, b):55 out = [a[0] * b[0] * 0] * (len(a) + len(b) - 1)56 for i, left in enumerate(a):57 for j, right in enumerate(b):58 out[i + j] = out[i + j] + left * right59 return out6061def scaled_basis(xs):62 basis = []63 for i, xi in enumerate(xs):64 term = [Fraction(1)]65 for j, xj in enumerate(xs):66 if i != j:67 term = polymul(term, [Fraction(-xj), Fraction(1)])68 term = [value / Fraction(xi - xj) for value in term]69 basis.append(term)70 scale = 171 for term in basis:72 for value in term:73 scale = scale * value.denominator // math.gcd(scale, value.denominator)74 return scale, [[int(value * scale) for value in term] for term in basis]7576def group_maps(d):77 cs = corners(d)78 index = {c: i for i, c in enumerate(cs)}79 maps = []80 for perm in permutations(range(d)):81 for pattern in range(1 << d):82 flip = tuple((pattern >> (d - 1 - a)) & 1 for a in range(d))83 maps.append(tuple(index[tuple(c[perm[a]] ^ flip[a] for a in range(d))]84 for c in cs))85 return maps8687def classify(d, total):88 maps = [or_tables(mapping) for mapping in group_maps(d)]89 rep_of = [-1] * total90 reps = []91 for code in range(total):92 if rep_of[code] >= 0:93 continue94 reps.append(code)95 for tables in maps:96 rep_of[block_or(tables, code)] = code97 return reps, rep_of9899def a129824(d):100 out = 1101 for w in range(d + 1):102 out *= 1 + math.comb(d, w)103 return out104105def sweep(d):106 size = 1 << d107 xs = list(range(1, d + 2))108 scale, basis = scaled_basis(xs)109 fit = [fill_tables(d, k) for k in xs]110 held = [(k, fill_tables(d, k)) for k in range(d + 2, KMAX + 1)]111 pop = sum_tables([1] * size)112 polys = []113 failures = 0114 for code in range(1 << size):115 ys = [block_sum(tables, code) for tables in fit]116 coef = [sum(ys[i] * basis[i][j] for i in range(d + 1)) for j in range(d + 1)]117 if any(c % scale for c in coef):118 failures += 1119 coef = [0] * (d + 1)120 else:121 coef = [c // scale for c in coef]122 if coef[d] != block_sum(pop, code):123 failures += 1124 for k, tables in held:125 if sum(c * k ** j for j, c in enumerate(coef)) != block_sum(tables, code):126 failures += 1127 break128 polys.append(tuple(coef))129 return polys, failures130131def polynomial_report():132 out = []133 for d in range(1, 5):134 polys, failures = sweep(d)135 reps, rep_of = classify(d, len(polys))136 members = {}137 for code, rep in enumerate(rep_of):138 members.setdefault(rep, []).append(code)139 lower = sum(1 for r in reps if len({polys[c][:-1] for c in members[r]}) > 1)140 lead = sum(1 for r in reps if len({polys[c][-1] for c in members[r]}) > 1)141 out.append((d, len(polys), len(reps), len(set(polys)), a129824(d),142 lower, lead, failures))143 return out144145def sig_bounds(d):146 return [math.comb(d, w) + 1 for w in range(d + 1)]147148def unrank(index, bounds):149 f = [0] * len(bounds)150 for i in range(len(bounds) - 1, -1, -1):151 f[i] = index % bounds[i]152 index //= bounds[i]153 return f154155def poly_of_sig(f, d):156 coef = [0] * (d + 1)157 for w, fw in enumerate(f):158 if fw:159 for j in range(d - w + 1):160 coef[w + j] += fw * math.comb(d - w, j)161 return coef162163def factor_list(coef):164 parts = []165 for g, mult in Poly(list(reversed(coef)), VAR, domain="ZZ").factor_list()[1]:166 parts.extend([g] * mult)167 return parts168169def quadratic_discs(coef):170 if not any(coef):171 return []172 out = []173 for g in factor_list(coef):174 if g.degree() == 2:175 a, b, c = g.all_coeffs()176 out.append(int(b * b - 4 * a * c))177 return out178179def evaluate_scaled(coef, q):180 top = len(coef) - 1181 return sum(c * (-1) ** i * q ** (top - i) for i, c in enumerate(coef))182183def peel_unit_linear(coef):184 c = list(coef)185 while len(c) > 1:186 lead = c[-1]187 taken = 0188 for q in range(1, abs(lead) + 1):189 if lead % q == 0 and evaluate_scaled(c, q) == 0:190 taken = q191 break192 if not taken:193 break194 out = [0] * (len(c) - 1)195 out[0] = c[0]196 for i in range(1, len(out)):197 out[i] = c[i] - taken * out[i - 1]198 c = out199 return c200201def remainder_disc(coef):202 rest = peel_unit_linear(coef)203 if len(rest) == 3:204 return rest[1] * rest[1] - 4 * rest[2] * rest[0]205 return None206207def disc_chunk(job):208 d, lo, hi = job209 bounds = sig_bounds(d)210 every, remainder = set(), set()211 for index in range(lo, hi):212 coef = poly_of_sig(unrank(index, bounds), d)213 every.update(quadratic_discs(coef))214 value = remainder_disc(coef)215 if value is not None:216 remainder.add(value)217 return every, remainder218219def disc_sets(d, workers):220 total = 1221 for b in sig_bounds(d):222 total *= b223 if d < PARALLEL_FROM:224 every, remainder = disc_chunk((d, 0, total))225 return every, remainder, total226 edges = [total * i // CHUNKS for i in range(CHUNKS + 1)]227 jobs = [(d, edges[i], edges[i + 1]) for i in range(CHUNKS)]228 every, remainder = set(), set()229 with ProcessPoolExecutor(max_workers=workers) as pool:230 for part, rest in pool.map(disc_chunk, jobs):231 every |= part232 remainder |= rest233 return every, remainder, total234235def gapless_run(found):236 out = []237 for v in range(-3, -LADDER_DEPTH - 1, -1):238 if v % 4 not in (0, 1):239 continue240 if v not in found:241 break242 out.append(v)243 return out244245def polytext(coef):246 return " + ".join("{}n^{}".format(c, i) if i > 1 else247 ("{}n".format(c) if i == 1 else str(c))248 for i, c in reversed(list(enumerate(coef))) if c)249250def quartic_splits(d):251 hits = []252 for f in product(*[range(b) for b in sig_bounds(d)]):253 rest = peel_unit_linear(poly_of_sig(list(f), d))254 if len(rest) != 5:255 continue256 parts = factor_list(rest)257 if len(parts) == 2 and all(g.degree() == 2 for g in parts):258 hits.append((f, rest, parts))259 return hits260261def solid_cells(mask, s):262 return [(x, y, z) for x in range(s) for y in range(s) for z in range(s)263 if (mask >> ((x & 1) | ((y & 1) << 1) | ((z & 1) << 2))) & 1]264265def census_values(mask, s):266 cells = solid_cells(mask, s)267 solid = set(cells)268 m = s + 2269 shift = m * m * m270 exposed = 0271 vertices, edges, faces = set(), set(), set()272 for x, y, z in cells:273 for p in ((x - 1, y, z), (x + 1, y, z), (x, y - 1, z),274 (x, y + 1, z), (x, y, z - 1), (x, y, z + 1)):275 if p not in solid:276 exposed += 1277 for a in (0, 1):278 xa = (x + a) * m279 for b in (0, 1):280 yb = (xa + y + b) * m281 for c in (0, 1):282 vertices.add(yb + z + c)283 faces.add((x * m + y) * m + z + a)284 faces.add(shift + (x * m + y + a) * m + z)285 faces.add(2 * shift + ((x + a) * m + y) * m + z)286 for b in (0, 1):287 for c in (0, 1):288 edges.add((x * m + y + b) * m + z + c)289 edges.add(shift + ((x + b) * m + y) * m + z + c)290 edges.add(2 * shift + ((x + b) * m + y + c) * m + z)291 v, e, f, fills = len(vertices), len(edges), len(faces), len(cells)292 return [fills, s ** 3 - fills, exposed, v, e, f, v - e + f - fills]293294def census_polys(mask, fam, scale, basis):295 side = (lambda i: 2 * i + 1) if fam == "odd" else (lambda i: 2 * i)296 samples = {i: census_values(mask, side(i)) for i in range(1, 7)}297 out = []298 bad = 0299 for obs in range(7):300 ys = [samples[1 + i][obs] for i in range(4)]301 coef = [Fraction(sum(ys[i] * basis[i][j] for i in range(4)), scale)302 for j in range(4)]303 for h in (5, 6):304 if sum(c * h ** j for j, c in enumerate(coef)) != samples[h][obs]:305 bad += 1306 out.append(coef)307 return out, bad, samples308309def closed_fill(mask, k):310 return sum(k ** (3 - bin(j).count("1")) * (k - 1) ** bin(j).count("1")311 for j in range(8) if (mask >> j) & 1)312313def factor_multisets(value, count, low=1):314 if count == 1:315 return [[value]] if value >= low else []316 out = []317 a = low318 while a ** count <= value:319 if value % a == 0:320 for tail in factor_multisets(value // a, count - 1, a):321 out.append([a] + tail)322 a += 1323 return out324325def divisor_shape(coef):326 if all(c == 0 for c in coef):327 return False328 k = 0329 while coef[k] == 0:330 k += 1331 q = [c / coef[k] for c in coef[k:]]332 while len(q) > 1 and q[-1] == 0:333 q.pop()334 if len(q) == 1 or any(value.denominator != 1 for value in q):335 return False336 q = [int(value) for value in q]337 if q[-1] <= 0:338 return False339 for multiset in factor_multisets(q[-1], len(q) - 1):340 built = [1]341 for a in multiset:342 built = polymul(built, [1, a])343 if built == q:344 return True345 return False346347def census_locked(mask, scale, basis):348 bad = 0349 drift = 0350 unlocked = False351 for fam in ("odd", "even"):352 polys, misses, samples = census_polys(mask, fam, scale, basis)353 bad += misses354 if fam == "odd":355 drift = sum(1 for i in range(1, 7)356 if samples[i][0] != closed_fill(mask, i + 1))357 for coef in polys:358 if divisor_shape(coef):359 unlocked = True360 return not unlocked, bad, drift361362def lock_predicate(mask):363 if mask & 1:364 return 0365 size = bin(mask).count("1")366 inner = sum(1 for c in range(8) for i in range(3)367 if (mask >> c) & 1 and c < (c ^ (1 << i)) and (mask >> (c ^ (1 << i))) & 1)368 if size == 5 and inner == 4:369 return 1370 if inner == 0 and (mask >> 7) & 1 and size != 2:371 return 2372 return 0373374def lock_report():375 scale, basis = scaled_basis([1, 2, 3, 4])376 path, edgeless, locked = [], [], []377 holdout = 0378 agree = 0379 grid = 0380 for mask in range(256):381 clause = lock_predicate(mask)382 if clause == 1:383 path.append(mask)384 elif clause == 2:385 edgeless.append(mask)386 is_locked, bad, drift = census_locked(mask, scale, basis)387 holdout += bad388 grid += drift == 0389 if is_locked:390 locked.append(mask)391 agree += (clause != 0) == is_locked392 return path, edgeless, sorted(path + edgeless), locked, agree, holdout, grid393394def main():395 workers = min(8, os.cpu_count() or 1)396 started = time.time()397 print("DOMAIN")398 print(" every design at D = 1..4 for the polynomial sweep, 2^(2^D) of them")399 print(" every signature at D = 2..6 for the discriminants")400 print(" every design at D = 3 for the lock, all 256")401 print()402 print("FILL POLYNOMIALS")403 for d, designs, classes, distinct, closed, lower, lead, bad in polynomial_report():404 print(" D = {}: {} designs, {} classes, {} distinct polynomials, "405 "A129824 = {}".format(d, designs, classes, distinct, closed))406 print(" lower coefficients split {} of {} classes, leading splits {} of {}, "407 "failures {}".format(lower, classes, lead, classes, bad))408 print()409 print("QUARTIC REMAINDERS SPLITTING INTO TWO QUADRATICS AT D = 4")410 hits = quartic_splits(4)411 print(" count {}".format(len(hits)))412 for f, rest, parts in hits:413 body = " ".join("({})".format(polytext(list(reversed(g.all_coeffs()))))414 for g in parts)415 print(" f={}: {} = {}".format(tuple(f), polytext(rest), body))416 print()417 print("QUADRATIC-FACTOR DISCRIMINANTS")418 for d in DISC_DIMS:419 every, remainder, total = disc_sets(d, workers)420 wide = gapless_run(every)421 narrow = gapless_run(remainder)422 print(" D = {}: {} signatures".format(d, total))423 print(" all quadratic factors: gapless -3..{}, length {}, deepest {}".format(424 wide[-1], len(wide), min(every)))425 print(" peeled remainder only: gapless -3..{}, length {}, deepest {}".format(426 narrow[-1], len(narrow), min(remainder)))427 print()428 print("THE CENSUS LOCK AT D = 3")429 path, edgeless, predicted, locked, agree, holdout, grid = lock_report()430 print(" path clause, origin-free with 5 corners and 4 edges: {} designs {}".format(431 len(path), path))432 print(" edgeless clause, origin-free with 111 and size not 2: {} designs {}".format(433 len(edgeless), edgeless))434 print(" predicate union: {} designs {}".format(len(predicted), predicted))435 print(" census rule locked: {} designs {}".format(len(locked), locked))436 print(" predicate and census rule agree on {} of 256 designs".format(agree))437 print(" held-out census fits that missed: {}".format(holdout))438 print(" grid fills match the closed form at k = 2..7 on {} of 256 designs".format(grid))439 print()440 print("wall time {:.1f} s on {} workers".format(time.time() - started, workers))441442if __name__ == "__main__":443 main()