dimers.py
23.3 kB · python · 534 lines
1import json2import os3import subprocess4import sys5import time6from math import ceil, comb, floor, isqrt, log78HERE = os.path.dirname(os.path.abspath(__file__))9DATA_DIR = os.path.join("data", os.path.relpath(HERE))1011DENSITY = 0.71213# DESIGN1415def digits(base, code):16 return [(i, j) for i in range(base) for j in range(base) if code >> (base * i + j) & 1]1718def level(base, code, n):19 F = digits(base, code)20 cells = {(0, 0)}21 for _ in range(n):22 cells = {(base * r + i, base * c + j) for (r, c) in cells for (i, j) in F}23 return cells2425def signed(base, code):26 return sum((-1) ** (i + j) for (i, j) in digits(base, code))2728def imbalance_law(base, code, n):29 s, fill = signed(base, code), len(digits(base, code))30 return s ** n if base % 2 else fill ** (n - 1) * s3132def neighbours(cells, cell):33 r, c = cell34 return [q for q in ((r + 1, c), (r - 1, c), (r, c + 1), (r, c - 1)) if q in cells]3536# MATCHING3738def has_perfect_matching(cells):39 black = [p for p in cells if (p[0] + p[1]) % 2 == 0]40 if 2 * len(black) != len(cells):41 return False42 mate = {}43 for root in black:44 seen = set()45 stack = [(root, iter(neighbours(cells, root)))]46 found = False47 while stack and not found:48 u, it = stack[-1]49 advanced = False50 for v in it:51 if v in seen:52 continue53 seen.add(v)54 if v not in mate:55 w = v56 for x, _ in reversed(stack):57 nxt = mate.get(x)58 mate[w] = x59 mate[x] = w60 w = nxt61 found = True62 break63 stack.append((mate[v], iter(neighbours(cells, mate[v]))))64 advanced = True65 break66 if not advanced and not found:67 stack.pop()68 if not found:69 return False70 return True7172def brute_count(cells, side):73 states = {0: 1}74 for p in range(side * side):75 r, c = divmod(p, side)76 nxt = {}77 for mask, ways in states.items():78 covered = mask & 179 rest = mask >> 180 if (r, c) not in cells or covered:81 if (r, c) not in cells and covered:82 continue83 nxt[rest] = nxt.get(rest, 0) + ways84 continue85 if c + 1 < side and (r, c + 1) in cells and not rest & 1:86 key = rest | 187 nxt[key] = nxt.get(key, 0) + ways88 if (r + 1, c) in cells:89 key = rest | 1 << (side - 1)90 nxt[key] = nxt.get(key, 0) + ways91 states = nxt92 return states.get(0, 0)9394# KASTELEYN9596def kasteleyn(cells, side, marked=None, holes=True):97 black = sorted(p for p in cells if (p[0] + p[1]) % 2 == 0)98 white = sorted(p for p in cells if (p[0] + p[1]) % 2 == 1)99 if len(black) != len(white):100 return None101 col = {p: k for k, p in enumerate(white)}102 below = {}103 for c in range(side):104 count = 0105 for r in range(side - 1, -1, -1):106 below[(r, c)] = count107 if (r, c) not in cells and holes:108 count += 1109 rows = []110 for b in black:111 row = {}112 for w in neighbours(cells, b):113 if w[0] == b[0]:114 left = min(b, w, key=lambda q: q[1])115 sign = (-1) ** below[left]116 else:117 sign = (-1) ** b[1]118 row[col[w]] = (sign, bool(marked and marked(b, w)))119 rows.append(row)120 return rows121122def gp_matrix(rows, value="x"):123 triples = []124 for i, row in enumerate(rows):125 for k, (sign, mark) in row.items():126 entry = f"{sign}*{value}" if mark else str(sign)127 triples.append(f"[{i + 1},{k + 1},{entry}]")128 n = len(rows)129 return f"E=[{','.join(triples)}];M=matrix({n},{n});for(k=1,#E,M[E[k][1],E[k][2]]=E[k][3]);"130131def gp(script):132 out = subprocess.run(["gp", "-q", "-D", "parisizemax=4000000000", "-D", "threadsizemax=4000000000"], input=script, capture_output=True, text=True, check=True)133 return out.stdout.split()134135def kasteleyn_count(cells, side, holes=True):136 rows = kasteleyn(cells, side, holes=holes)137 if rows is None:138 return 0139 if not rows:140 return 1141 return abs(int(gp(gp_matrix(rows) + "print(matdet(M));")[0]))142143# STRUCTURE144145def components(cells):146 seen, out = set(), []147 for p in cells:148 if p in seen:149 continue150 comp, k = [p], 0151 seen.add(p)152 while k < len(comp):153 for q in neighbours(cells, comp[k]):154 if q not in seen:155 seen.add(q)156 comp.append(q)157 k += 1158 out.append(comp)159 return out160161def sealed_cells(F, side):162 wide = any((i, j + 1) in F for (i, j) in F)163 tall = any((i + 1, j) in F for (i, j) in F)164 top = side - 1165 for C in components(F):166 if sum((-1) ** (i + j) for (i, j) in C) == 0:167 continue168 reach = any(169 (wide and j == top and (i, 0) in F) or (wide and j == 0 and (i, top) in F)170 or (tall and i == top and (0, j) in F) or (tall and i == 0 and (top, j) in F)171 for (i, j) in C172 )173 if not reach:174 return True175 return False176177def sealed(base, code, m=1):178 return sealed_cells(level(base, code, m), base ** m)179180def one_way(base, code):181 F, top = set(digits(base, code)), base - 1182 for flip in (False, True):183 G = {(j, i) for (i, j) in F} if flip else F184 tall = any((i + 1, j) in G for (i, j) in G) and any((0, j) in G and (top, j) in G for j in range(base))185 if not tall:186 return G187 return None188189def run_returns(base, code):190 G, top = one_way(base, code), base - 1191 rows = [i for i in range(base) if (i, 0) in G and (i, top) in G] if any((i, j + 1) in G for (i, j) in G) else []192 subsets = [frozenset(r for k, r in enumerate(rows) if mask >> k & 1) for mask in range(2 ** len(rows))]193 step = {A: [B for B in subsets if has_perfect_matching(G - {(i, 0) for i in A} - {(i, top) for i in B})] for A in subsets}194 seen, frontier = set(), [frozenset()]195 while frontier:196 nxt = []197 for A in frontier:198 for B in step[A]:199 if not B:200 return True201 if B not in seen:202 seen.add(B)203 nxt.append(B)204 frontier = nxt205 return False206207def orbit(base, code):208 F, top, out = digits(base, code), base - 1, set()209 for k in range(8):210 image = []211 for i, j in F:212 for _ in range(k % 4):213 i, j = j, top - i214 image.append((i, top - j) if k >= 4 else (i, j))215 out.add(sum(1 << (base * i + j) for i, j in image))216 return min(out)217218def picture(base, code):219 F = set(digits(base, code))220 return " ".join("".join("#" if (i, j) in F else "." for j in range(base)) for i in range(base))221222# VERBS223224LATE = ((5, 15571455), (5, 19627890), (5, 19757811), (5, 19920882), (5, 32709486), (5, 32715747), (5, 32715771), (5, 32912238), (7, 136308971855667))225226def first_level(base, code, top, cap):227 for n in range(1, top + 1):228 if len(digits(base, code)) ** n > cap:229 return "cap"230 if has_perfect_matching(level(base, code, n)):231 return n232 return None233234def imbalance():235 for base, top in ((2, 6), (3, 4)):236 values, wrong = {}, 0237 for code in range(1, 2 ** (base * base)):238 s = signed(base, code)239 values[s] = values.get(s, 0) + 1240 for n in range(1, top + 1):241 cells = level(base, code, n)242 direct = sum(1 if (r + c) % 2 == 0 else -1 for (r, c) in cells)243 wrong += direct != imbalance_law(base, code, n)244 zero = values.get(0, 0)245 print(f"base {base}: {2 ** (base * base) - 1} codes, levels 1..{top}, law against the coordinate count wrong {wrong} times")246 print(f" s = 0 at {zero} codes, C({base * base},{base * base // 2}) - 1 = {comb(base * base, base * base // 2) - 1}")247 print(f" codes per s: {dict(sorted(values.items()))}")248249def census():250 for base, top in ((2, 4), (3, 4), (4, 3)):251 tally, never = {}, []252 for code in range(1, 2 ** (base * base)):253 if signed(base, code) != 0:254 continue255 n0 = first_level(base, code, top, 7000)256 tally[n0] = tally.get(n0, 0) + 1257 if n0 is None:258 never.append(code)259 print(f"base {base}: first tileable level over the s = 0 codes, levels 1..{top}: {tally}")260 print(f" untileable through level {top}: {len(never)} codes in {len({orbit(base, c) for c in never})} orbits")261 rest = never262 for m in (1, 2, 3):263 hit = [c for c in rest if len(digits(base, c)) ** m <= 70000 and sealed(base, c, m)]264 rest = [c for c in rest if c not in set(hit)]265 print(f" a sealed unbalanced component at level {m}: {len(hit)} more codes; {len(rest)} codes in {len({orbit(base, c) for c in rest})} orbits left{', e.g. ' + str(rest[:6]) if rest and m == 3 else ''}")266 if base == 3:267 print(f" base 3 untileable codes: {never}")268 if rest:269 print(f" the {len(rest)} left tile level {top + 1}: {sum(has_perfect_matching(level(base, c, top + 1)) for c in rest)}, carry a sealed unbalanced component at level {top + 1}: {sum(sealed(base, c, top + 1) for c in rest)}")270 ways = [c for c in range(1, 2 ** (base * base)) if one_way(base, c) is not None]271 tiling = [c for c in ways if has_perfect_matching(level(base, c, 1))]272 print(f" blocks meeting in one direction only: {len(ways)} codes, the {len(tiling)} that tile level 1 all have a run walk back to the empty set: {all(run_returns(base, c) for c in tiling)}")273 dead = [c for c in rest if one_way(base, c) is not None and not run_returns(base, c)]274 rest = [c for c in rest if c not in set(dead)]275 print(f" one direction only and no run walk back to the empty set: {len(dead)} codes in {len({orbit(base, c) for c in dead})} orbits {sorted({orbit(base, c) for c in dead})}; {len(rest)} codes in {len({orbit(base, c) for c in rest})} orbits left {sorted({orbit(base, c) for c in rest})}")276 print("late codes: no tiling at level 1, a tiling at level 2")277 for base, code in LATE:278 one, two = level(base, code, 1), level(base, code, 2)279 T2 = kasteleyn_count(two, base * base)280 brute = brute_count(two, base * base) if base == 5 else "-"281 print(f" base {base} code {code} [{picture(base, code)}] fill {len(one)} s {signed(base, code)} sealed {sealed(base, code)} T(1) {kasteleyn_count(one, base)} {brute_count(one, base)} match(2) {has_perfect_matching(two)} T(2) {T2} brute {brute} orbit {orbit(base, code)}")282283def deficiency(cells):284 mate, size = {}, 0285 for root in (p for p in cells if (p[0] + p[1]) % 2 == 0):286 seen, stack, found = set(), [(root, iter(neighbours(cells, root)))], False287 while stack and not found:288 u, it = stack[-1]289 advanced = False290 for v in it:291 if v in seen:292 continue293 seen.add(v)294 if v not in mate:295 w = v296 for x, _ in reversed(stack):297 nxt = mate.get(x)298 mate[w], mate[x] = x, w299 w = nxt300 found = True301 break302 stack.append((mate[v], iter(neighbours(cells, mate[v]))))303 advanced = True304 break305 if not advanced and not found:306 stack.pop()307 size += found308 return len(cells) - 2 * size309310def search(draws=None):311 import random312 plan = draws or ((5, 1500000), (6, 300000), (7, 100000))313 for base, count in plan:314 rng = random.Random(base)315 codes = {code for code in (sum(1 << k for k in range(base * base) if rng.random() < DENSITY) for _ in range(count)) if code and signed(base, code) == 0}316 untileable = unsealed = tested3 = third = 0317 late = set()318 for code in sorted(codes):319 one = level(base, code, 1)320 if has_perfect_matching(one):321 continue322 untileable += 1323 if sealed(base, code):324 continue325 unsealed += 1326 d2 = deficiency(level(base, code, 2))327 if d2 == 0:328 late.add(code)329 continue330 if d2 < len(one) * deficiency(one) and len(one) ** 3 <= 16000:331 tested3 += 1332 third += has_perfect_matching(level(base, code, 3))333 orbits = sorted({orbit(base, c) for c in late})334 print(f"base {base}: {count} draws, {len(codes)} distinct codes with s = 0, {untileable} untileable at level 1, {unsealed} of those unsealed, {len(late)} tile first at level 2 ({len(orbits)} orbits {orbits}), {tested3} tried at level 3 (fill^3 <= 16000, level-2 deficiency below fill times level-1 deficiency), {third} tile first at level 3")335 if base == 5:336 neighbourhood(sorted(set(orbits) | {orbit(b, code) for b, code in LATE if b == 5}))337338def neighbourhood(seeds, base=5, radius=4):339 from itertools import combinations340 seen, late, tested, third = set(), set(), 0, 0341 for seed in seeds:342 for r in range(1, radius + 1):343 for flip in combinations(range(base * base), r):344 code = seed345 for k in flip:346 code ^= 1 << k347 key = orbit(base, code)348 if key in seen:349 continue350 seen.add(key)351 if signed(base, code) != 0 or has_perfect_matching(level(base, code, 1)) or sealed(base, code):352 continue353 if has_perfect_matching(level(base, code, 2)):354 late.add(key)355 continue356 if len(digits(base, code)) ** 3 <= 16000:357 tested += 1358 third += has_perfect_matching(level(base, code, 3))359 print(f" neighbourhood of {len(seeds)} late base-{base} orbits to Hamming distance {radius}: {len(seen)} orbits, {len(late)} late orbits, {tested} orbits untileable at levels 1 and 2 with fill^3 <= 16000 tried at level 3, {third} tile there")360361CODES = ((2, 15, 6), (3, 63, 4), (3, 495, 4), (3, 255, 4))362363def remember(key, T):364 path = os.path.join(DATA_DIR, "counts.json")365 store = json.load(open(path)) if os.path.exists(path) else {}366 store[key] = str(T)367 os.makedirs(DATA_DIR, exist_ok=True)368 json.dump(store, open(path, "w"))369370def cached_count(base, code, n):371 path = os.path.join(DATA_DIR, "counts.json")372 store = json.load(open(path)) if os.path.exists(path) else {}373 key = f"{base},{code},{n}"374 if key in store:375 return int(store[key])376 T = kasteleyn_count(level(base, code, n), base ** n)377 remember(key, T)378 return T379380def rectangle(side):381 script = f"default(realprecision,{side * side // 3 + 50});print(round(prod(j=1,{side // 2},prod(k=1,{side // 2},4*cos(Pi*j/{side + 1})^2+4*cos(Pi*k/{side + 1})^2))));"382 return int(gp(script)[0])383384def fibonacci(n):385 a, b = 0, 1386 for _ in range(n):387 a, b = b, a + b388 return a389390def control():391 import random392 rng = random.Random(1)393 agree = textbook = tileable = 0394 for _ in range(300):395 side = rng.randint(2, 7)396 cells = {(r, c) for r in range(side) for c in range(side) if rng.random() < 0.8}397 T = brute_count(cells, side)398 tileable += T > 0399 agree += kasteleyn_count(cells, side) == T and has_perfect_matching(cells) == (T > 0)400 textbook += kasteleyn_count(cells, side, holes=False) == T401 print(f"control: 300 random cell sets, side 2..7, density 0.8: Kasteleyn with the hole signs and the matching test agree with brute force at {agree}, the column signs alone at {textbook}; {tileable} sets tile")402 for base, code, n in ((3, 495, 1), (3, 495, 2)):403 print(f" carpet level {n}: column signs alone give |det K| = {kasteleyn_count(level(base, code, n), base ** n, holes=False)}")404405def count():406 control()407 for base, code, top in CODES:408 print(f"base {base} code {code} [{picture(base, code)}]")409 for n in range(1, top + 1):410 start = time.time()411 cells = level(base, code, n)412 T = kasteleyn_count(cells, base ** n)413 checks = []414 if base ** n <= 9:415 checks.append(f"brute {brute_count(cells, base ** n) == T}")416 if code == 15:417 checks.append(f"product formula {rectangle(2 ** n) == T}")418 if code == 63:419 checks.append(f"F(3^n+1)^(2^(n-1)) {fibonacci(3 ** n + 1) ** (2 ** (n - 1)) == T}")420 root = isqrt(T)421 two = (T & -T).bit_length() - 1422 odd = isqrt(T >> two)423 shown = str(T) if len(str(T)) <= 60 else f"{str(T)[:20]}...{str(T)[-20:]} ({len(str(T))} digits)"424 square = odd * odd == T >> two425 print(f" level {n}: cells {len(cells)}, T = {shown}, square {root * root == T}, T = 2^{two} times {'the square of ' + (str(odd) if len(str(odd)) <= 40 else str(len(str(odd))) + '-digit ' + str(odd)[:12] + '...') if square else 'a non-square'}, {', '.join(checks)}, {time.time() - start:.1f} s")426 remember(f"{base},{code},{n}", T)427428def crossing_polynomial(base, code, n):429 cells, unit = level(base, code, n), base ** (n - 1)430 rows = kasteleyn(cells, base ** n, lambda p, q: (p[0] // unit, p[1] // unit) != (q[0] // unit, q[1] // unit))431 edges = sum(1 for row in rows for _, mark in row.values() if mark)432 script = gp_matrix(rows, "x") + f"v=vector({edges + 1},t,matdet(subst(M,x,t-1)));P=polinterpolate(vector({edges + 1},t,t-1),v);P=P*sign(subst(P,x,1));c=content(P);print(vector(poldegree(P)+1,k,polcoeff(P,k-1)));print(c);print(issquare(P/c));"433 out = gp(script)434 return edges, [int(a) for a in "".join(out[:-2]).strip("[]").split(",")], int(out[-2]), out[-1] == "1"435436def cross():437 for base, code, n in ((3, 495, 2), (3, 495, 3), (3, 255, 2), (3, 255, 3), (2, 15, 2), (2, 15, 3), (2, 15, 4)):438 start = time.time()439 fill = len(digits(base, code))440 edges, N, content, square = crossing_polynomial(base, code, n)441 T, below, unit = cached_count(base, code, n), cached_count(base, code, n - 1), cached_count(base, code, 1)442 mean = sum(k * a for k, a in enumerate(N)) / sum(N)443 support = [k for k, a in enumerate(N) if a]444 mode = max(range(len(N)), key=lambda k: N[k])445 print(f"base {base} code {code} level {n}: {edges} edges cross the {fill} blocks of level {n - 1}")446 print(f" sum N_k = T({n}) {sum(N) == T}, N_0 = T({n - 1})^{fill} {N[0] == below ** fill}, small blocks only T(1)^({fill}^{n - 1}) = {unit ** fill ** (n - 1)}")447 print(f" N_k nonzero exactly at the even k from 0 to {support[-1]}: {support == list(range(0, support[-1] + 1, 2))}, mode {mode}, mean {mean:.6f}, mean per crossing edge {mean / edges:.6f}")448 print(f" the polynomial is its content {content} times a square: {square}")449 if (base, code, n) == (3, 495, 2):450 inner = [0] * 9451 for k in range(5):452 inner[2 * k] += comb(4, k) * 2 ** (4 - k)453 inner[8] += 1454 outer = [sum(inner[i] * inner[k - i] for i in range(max(0, k - 8), min(k, 8) + 1)) for k in range(17)]455 print(f" P_2(x) = ((x^2 + 2)^4 + x^8)^2: {outer == N}")456 print(f" share of tilings that respect the blocks N_0/T = {N[0] / T:.6e}, T/N_0 = {T / N[0]:.6f}")457 print(f" N_k = {N if len(N) <= 30 else N[:6]} {'' if len(N) <= 30 else '...'}, {time.time() - start:.1f} s")458459def exposure(base, code, depth):460 F, top = set(digits(base, code)), base - 1461 step = {}462 for E in range(16):463 out = []464 for i, j in F:465 child = 0466 child |= 1 if (i > 0 and (i - 1, j) not in F) or (i == 0 and (E & 1 or (top, j) not in F)) else 0467 child |= 2 if (i < top and (i + 1, j) not in F) or (i == top and (E & 2 or (0, j) not in F)) else 0468 child |= 4 if (j > 0 and (i, j - 1) not in F) or (j == 0 and (E & 4 or (i, top) not in F)) else 0469 child |= 8 if (j < top and (i, j + 1) not in F) or (j == top and (E & 8 or (i, 0) not in F)) else 0470 out.append(child)471 step[E] = out472 share, history = {15: 1.0}, []473 for _ in range(depth):474 nxt = {}475 for E, p in share.items():476 for child in step[E]:477 nxt[child] = nxt.get(child, 0.0) + p / len(F)478 share = nxt479 history.append(dict(share))480 return history, all(step[child][k] == child for E in range(16) for k, child in enumerate(step[E]))481482def hadamard(share):483 total = 0.0484 for E, p in share.items():485 d = 4 - bin(E).count("1")486 if d == 0:487 return float("-inf")488 total += p * log(d)489 return total / 4490491def contact(base, code):492 F, top = set(digits(base, code)), base - 1493 rows = sum((i, 0) in F and (i, top) in F for i in range(base)) if any((i, j + 1) in F for i, j in F) else 0494 cols = sum((0, j) in F and (top, j) in F for j in range(base)) if any((i + 1, j) in F for i, j in F) else 0495 return max(rows, cols)496497def down(x, places=9):498 return f"{floor(x * 10 ** places) / 10 ** places:.{places}f}"499500def up(x, places=6):501 return f"{ceil(x * 10 ** places) / 10 ** places:.{places}f}"502503def growth():504 G = 0.915965594177219015054603514932384110774505 phi = (1 + 5 ** 0.5) / 2506 exact = {15: ("G/pi", G / 3.141592653589793), 63: ("log(phi)/2", log(phi) / 2)}507 for base, code, top in CODES:508 fill = len(digits(base, code))509 print(f"base {base} code {code} [{picture(base, code)}] fill {fill}")510 L = [log(cached_count(base, code, n)) / fill ** n for n in range(1, top + 1)]511 for n in range(1, top + 1):512 inc = L[n - 1] - L[n - 2] if n > 1 else None513 ratio = (L[n - 1] - L[n - 2]) / (L[n - 2] - L[n - 3]) if n > 2 else None514 print(f" level {n}: log T / cells = {down(L[n - 1])}" + (f", increment {inc:.6f}" if inc else "") + (f", increment ratio {ratio:.4f}" if ratio else ""))515 history, idempotent = exposure(base, code, 400)516 H = [hadamard(history[n - 1]) for n in (1, 2, 4, 10, 20, 40, 60, 400)]517 print(f" each digit's exposure map is idempotent: {idempotent}; Hadamard bound H_n at n = 1, 2, 4, 10, 20, 40, 60, 400: {', '.join(f'{h:.13f}' for h in H)}")518 print(f" bracket: {down(L[-1], 6)} <= theta <= {up(H[-1] + 1e-9)}")519 q, last = contact(base, code) / fill, (L[-1] - L[-2]) / (L[-2] - L[-3])520 tail = lambda r: (L[-1] - L[-2]) * r / (1 - r)521 print(f" contact ratio q = {contact(base, code)}/{fill}; geometric tail from level {top}: with q {L[-1] + tail(q):.6f}, with the last ratio {L[-1] + tail(last):.6f}")522 if code in exact:523 name, value = exact[code]524 print(f" closed form {name} = {value:.12f}")525 print(f"code 27 [{picture(3, 27)}] has no crossing edge, T(n) = 2^(4^(n-1)): {all(cached_count(3, 27, n) == 2 ** 4 ** (n - 1) for n in (1, 2, 3))}, theta = log(2)/4 = {log(2) / 4:.12f}")526527VERBS = {"imbalance": imbalance, "census": census, "search": search, "count": count, "cross": cross, "growth": growth}528529if __name__ == "__main__":530 for name in sys.argv[1:] or list(VERBS):531 start = time.time()532 print(f"# {name}")533 VERBS[name]()534 print(f"# {name} {time.time() - start:.1f} s")