life.py
15.8 kB · python · 399 lines
1import time2from math import comb, log234import numpy as np56BINOM = np.array([comb(8, k) for k in range(9)], dtype=np.int64)7POP8 = np.array([i.bit_count() for i in range(256)], dtype=np.int64)8POP9 = np.array([i.bit_count() for i in range(512)], dtype=np.int64)9SYM = ((np.arange(512, dtype=np.int64)[:, None] >> POP8[None, :]) & 1).astype(np.uint8)10SUBSET = ((np.arange(512, dtype=np.int64)[:, None] >> np.arange(9)[None, :]) & 1)11FILLS = SUBSET @ BINOM12BITS = ((np.arange(512, dtype=np.int64)[:, None] >> np.arange(9)[None, :]) & 1)13EDGE = SUBSET @ np.array([comb(7, k - 1) if k else 0 for k in range(9)], dtype=np.int64)14ORDER = np.argsort(POP8, kind="stable")15STARTS = np.concatenate([[0], np.cumsum(BINOM)[:-1]])16NAMED = (17 ((3,), (2, 3)),18 ((3, 6), (2, 3)),19 ((2,), ()),20 ((3,), (0, 1, 2, 3, 4, 5, 6, 7, 8)),21 ((3, 6, 7, 8), (3, 4, 6, 7, 8)),22 ((3, 5, 6, 7, 8), (5, 6, 7, 8)),23 ((1, 3, 5, 7), (1, 3, 5, 7)),24 ((4, 6, 7, 8), (3, 5, 6, 7, 8)),25 ((3, 6), (1, 2, 5)),26 ((3, 6, 8), (2, 4, 5)),27 ((), ()),28 ((0, 1, 2, 3, 4, 5, 6, 7, 8), (0, 1, 2, 3, 4, 5, 6, 7, 8)),29 ((), (0, 1, 2, 3, 4, 5, 6, 7, 8)),30 ((0, 1, 2, 3, 4, 5, 6, 7, 8), ()),31)32LIFE = ((3,), (2, 3))33HORIZON = 6434POWERS = (1, 2, 4, 8, 16, 32, 64)35A071053 = (1, 3, 3, 5, 3, 9, 5, 11, 3, 9, 9, 15, 5, 15, 11, 21, 3, 9, 9, 15)36A246035 = (1, 9, 9, 25, 9, 81, 25, 121, 9, 81, 81, 225, 25, 225, 121, 441, 9, 81, 81, 225)37A160239 = (1, 8, 8, 24, 8, 64, 24, 112, 8, 64, 64, 192, 24, 192, 112, 416, 8, 64, 64, 192)3839def moore(d):40 mask = np.ones(3 ** d, dtype=np.uint8)41 center = 042 for _ in range(d):43 center = center * 3 + 144 mask[center] = 045 return mask.reshape([3] * d)4647def residue_corners(d, base=2):48 return [[(i // base ** (d - 1 - j)) % base for j in range(d)] for i in range(base ** d)]4950def tile_from_code(code, d, side=3):51 corners = residue_corners(d)52 filled = [(code >> i) & 1 for i in range(len(corners))]53 out = np.zeros([side] * d, dtype=np.uint8)54 for flat in range(side ** d):55 v = tuple((flat // side ** (d - 1 - j)) % side for j in range(d))56 corner = 057 for coord in v:58 corner = corner * 2 + coord % 259 out[v] = filled[corner]60 return out6162def subset_code(counts):63 code = 064 for k in counts:65 code |= 1 << k66 return code6768def rule_table(birth, survive):69 return np.concatenate([SYM[subset_code(birth)], SYM[subset_code(survive)]])7071def rule_name(birth, survive):72 return "B" + "".join(str(k) for k in sorted(birth)) + "/S" + "".join(str(k) for k in sorted(survive))7374def design_code(table):75 code = 076 for x in np.flatnonzero(table):77 code |= 1 << int(x)78 return code7980def block_table(table):81 out = np.zeros(512, dtype=np.uint8)82 for y in range(512):83 bits = [(y >> (8 - p)) & 1 for p in range(9)]84 centre = bits[4]85 count = sum(bits) - centre86 out[y] = table[centre * 256 + (2 ** count - 1)]87 return out8889def anf(tables):90 a = np.atleast_2d(tables).astype(np.uint8).copy()91 for b in range(a.shape[1].bit_length() - 1):92 step = 1 << b93 view = a.reshape(a.shape[0], -1, 2, step)94 view[:, :, 1, :] ^= view[:, :, 0, :]95 return a9697def degrees(tables):98 a = anf(tables)99 pops = POP9 if a.shape[1] == 512 else POP8100 return np.where(a.astype(bool), pops[None, :], -1).max(axis=1)101102def walsh(table):103 w = table.astype(np.int64).copy()104 for b in range(9):105 step = 1 << b106 view = w.reshape(-1, 2, step)107 low = view[:, 0, :].copy()108 high = view[:, 1, :].copy()109 view[:, 0, :] = low + high110 view[:, 1, :] = low - high111 return w112113def level_sums(table):114 sums = np.zeros(10, dtype=np.int64)115 np.add.at(sums, POP9, walsh(table))116 return sums117118def level_sums_closed(table):119 weights = np.bincount(POP9[table.astype(bool)], minlength=10)120 poly = np.zeros(10, dtype=np.int64)121 for w, count in enumerate(weights):122 if not count:123 continue124 term = np.array([1], dtype=np.int64)125 for _ in range(9 - w):126 term = np.convolve(term, [1, 1])127 for _ in range(w):128 term = np.convolve(term, [1, -1])129 poly += count * term130 return poly131132def is_pin(table):133 support = np.flatnonzero(table)134 if support.size == 0:135 return False136 bits = ((support[:, None] >> np.arange(9)[None, :]) & 1).sum(axis=0)137 fixed = int(((bits == 0) | (bits == support.size)).sum())138 return support.size == 1 << (9 - fixed)139140def is_level_set(table):141 for t in range(512):142 classes = POP9[np.arange(512) ^ t]143 lows = np.zeros(10, dtype=np.int64)144 highs = np.zeros(10, dtype=np.int64)145 np.add.at(lows, classes, table)146 np.add.at(highs, classes, 1 - table)147 if not np.any((lows > 0) & (highs > 0)):148 return True149 return False150151def pin_brute():152 out = np.zeros((512, 512), dtype=bool)153 for b in range(512):154 tables = np.concatenate([np.repeat(SYM[b][None, :], 512, axis=0), SYM], axis=1).astype(np.int64)155 sizes = tables.sum(axis=1)156 fixed = ((tables @ BITS == 0) | (tables @ BITS == sizes[:, None])).sum(axis=1)157 out[b] = (sizes > 0) & (sizes == (1 << (9 - fixed)))158 return out159160def pin_grid():161 fill = FILLS[:, None] + FILLS[None, :]162 outer = EDGE[:, None] + EDGE[None, :]163 centre = np.broadcast_to(FILLS[None, :], fill.shape)164 fixed = 8 * ((outer == 0) | (outer == fill)) + ((centre == 0) | (centre == fill))165 return fill == (1 << (9 - fixed))166167def level_set_grid():168 flags = np.zeros((512, 512), dtype=bool)169 xs = np.arange(512, dtype=np.int64)170 classes = POP9[xs[None, :] ^ xs[:, None]]171 for chosen in range(1024):172 levels = ((chosen >> classes) & 1).astype(np.int64)173 sums_b = np.add.reduceat(levels[:, :256][:, ORDER], STARTS, axis=1)174 sums_s = np.add.reduceat(levels[:, 256:][:, ORDER], STARTS, axis=1)175 good_b = ((sums_b == 0) | (sums_b == BINOM[None, :])).all(axis=1)176 good_s = ((sums_s == 0) | (sums_s == BINOM[None, :])).all(axis=1)177 keep = good_b & good_s178 if not keep.any():179 continue180 codes_b = ((sums_b > 0) << np.arange(9)[None, :]).sum(axis=1)181 codes_s = ((sums_s > 0) << np.arange(9)[None, :]).sum(axis=1)182 flags[codes_b[keep], codes_s[keep]] = True183 return flags184185def degree_grid():186 grid = np.zeros((512, 512), dtype=np.int8)187 for b in range(512):188 tables = np.concatenate([np.repeat(SYM[b][None, :], 512, axis=0), SYM], axis=1)189 grid[b] = degrees(tables)190 return grid191192def step(grid, birth, survive, mask):193 pad = np.pad(grid, 1)194 counts = np.zeros(grid.shape, dtype=np.int64)195 for di in range(3):196 for dj in range(3):197 if mask[di, dj]:198 counts += pad[di:di + grid.shape[0], dj:dj + grid.shape[1]]199 return np.where(grid == 1, np.isin(counts, survive), np.isin(counts, birth)).astype(np.uint8)200201def step150(row):202 pad = np.pad(row, 1)203 return ((pad[:-2] + pad[1:-1] + pad[2:]) % 2).astype(np.uint8)204205def t1_moore():206 print("T1 moore = carpet")207 for d in (1, 2, 3):208 code = (1 << (1 << d)) - 1 - (1 << ((1 << d) - 1))209 tile = tile_from_code(code, d)210 assert np.array_equal(tile, moore(d)), d211 print(f" D={d} code {code} level-1 side-3 tile = moore({d}), fill {int(tile.sum())} of {3 ** d}")212 assert (1 << 4) - 1 - (1 << 3) == 7213 print(" the plane case is bang dim 2, code 7, fill 8 of 9, dimension log(8)/log(3) = 1.892789")214215def t2_life():216 print("T2 life as a design")217 table = rule_table(*LIFE)218 fill = int(table.sum())219 code = design_code(table)220 lam = fill / 512221 print(f" {rule_name(*LIFE)} fill {fill} of 512 = {comb(8, 3)} + {comb(8, 2)} + {comb(8, 3)}, lambda {fill}/512")222 print(f" code {fill} bits, hex {code:0128x}")223 print(f" code under the 3x3 block order, hex {design_code(block_table(table)):0128x}")224 assert bin(code).count("1") == fill225 deg = int(degrees(table[None, :])[0])226 print(f" GF(2) degree {deg}, monomials {int(anf(table[None, :]).sum())}")227 sums = level_sums(table)228 closed = level_sums_closed(table)229 assert np.array_equal(sums, closed)230 print(" walsh level sums " + " ".join(f"S{k}={int(v)}" for k, v in enumerate(sums)))231 energy = int((walsh(table).astype(np.int64) ** 2).sum())232 assert energy == 512 * fill233 print(f" parseval sum of W(S)^2 over the 512 subsets {energy} = 512 * fill = {512 * fill}")234 print(f" popcount 4 both ways: f(c=1,|n|=3) = {int(table[256 + 7])}, f(c=0,|n|=4) = {int(table[15])}")235 print(f" level set {is_level_set(table)}, pin {is_pin(table)}, genus compound")236 print(f" dimension log2(fill) = {log2(fill):.6f} = 9 + log2({fill}/512) = {9 + log2(lam):.6f}")237 assert abs(log2(fill) - (9 + log2(lam))) < 1e-12238 print(f" outer-totalistic rules 2^18 = {1 << 18}, totalistic 2^10 = {1 << 10}")239240def t3_census(pins, levels):241 print("T3 the lambda census")242 dist = np.zeros(257, dtype=np.int64)243 dist[0] = 1244 for k in range(9):245 shifted = np.zeros_like(dist)246 shifted[BINOM[k]:] = dist[:257 - BINOM[k]]247 dist = dist + shifted248 hist = np.convolve(dist, dist)249 brute = np.bincount((FILLS[:, None] + FILLS[None, :]).ravel(), minlength=513)250 assert np.array_equal(hist, brute)251 assert int(hist.sum()) == 1 << 18252 seen = np.flatnonzero(hist)253 missing = [v for v in range(513) if hist[v] == 0]254 print(f" histogram total {int(hist.sum())}, distinct fills {seen.size} of 513, mirror symmetric {np.array_equal(hist, hist[::-1])}")255 print(f" unreachable fills {missing}")256 print(f" fill 0 and 512 count {int(hist[0])} each, peak fill {int(np.argmax(hist))} count {int(hist.max())}")257 print(f" rules sharing the fill of {rule_name(*LIFE)}: {int(hist[140])}")258 quarters = [int(hist[v]) for v in (64, 128, 140, 256, 384, 448)]259 print(f" counts at fills 64 128 140 256 384 448: {quarters}")260 return hist261262def t3_rules(grid, pins, levels, hist):263 print(" rule | fill | lambda | dimension | deg | genus | count at that fill")264 for birth, survive in NAMED:265 b, s = subset_code(birth), subset_code(survive)266 fill = int(FILLS[b] + FILLS[s])267 genus = "iso" if levels[b, s] else ("axis" if pins[b, s] else "comp")268 dim = f"{log2(fill):.6f}" if fill else "none"269 print(f" {rule_name(birth, survive)} | {fill} | {fill}/512 | {dim} | {int(grid[b, s])} | {genus} | {int(hist[fill])}")270271def t3_genus(pins, levels, grid):272 total = 1 << 18273 iso = int(levels.sum())274 axis = int((pins & ~levels).sum())275 comp = total - iso - axis276 print(f" genus over the 2^18: iso {iso}, axis only {axis}, compound {comp}")277 totalistic = set()278 for chosen in range(1024):279 b = subset_code(k for k in range(9) if (chosen >> k) & 1)280 s = subset_code(k for k in range(9) if (chosen >> (k + 1)) & 1)281 assert levels[b, s]282 totalistic.add((b, s))283 assert len(totalistic) == 1024284 print(f" totalistic rules {len(totalistic)} of the 2^18, every one a level set of the full popcount")285 hist = np.bincount(grid.ravel() + 1, minlength=11)286 print(" degree histogram " + " ".join(f"{k - 1}:{int(v)}" for k, v in enumerate(hist) if v))287288def t4_affine(grid):289 print("T4 the affine life-like rules")290 sym_deg = degrees(SYM)291 predicted = np.zeros((512, 512), dtype=np.int8)292 for b in range(512):293 diff = b ^ np.arange(512)294 alt = sym_deg[diff]295 predicted[b] = np.where(diff == 0, sym_deg[b], np.maximum(sym_deg[b], 1 + alt))296 assert np.array_equal(predicted, grid)297 print(" deg(B,S) = deg(B) when B = S, else max(deg(B), 1 + deg(B xor S)): holds on all 2^18")298 flat = np.argwhere(grid <= 1)299 names = []300 for b, s in flat:301 birth = [k for k in range(9) if (b >> k) & 1]302 survive = [k for k in range(9) if (s >> k) & 1]303 names.append(rule_name(birth, survive))304 print(f" degree <= 1 rules: {len(names)}")305 for name in names:306 print(f" {name}")307 assert len(names) == 8308 expected = {"B/S", "B012345678/S012345678", "B/S012345678", "B012345678/S",309 "B1357/S1357", "B1357/S02468", "B02468/S1357", "B02468/S02468"}310 assert set(names) == expected311 print(" the four degenerate rules and the four parity rules, no others")312 sym_hist = np.bincount(sym_deg + 1, minlength=10)313 print(" degree histogram of the 512 count sets " + " ".join(f"{k - 1}:{int(v)}" for k, v in enumerate(sym_hist) if v))314 for d in range(9):315 assert int((sym_deg <= d).sum()) == 1 << (d + 1), d316 for d in range(1, 9):317 assert int((grid <= d).sum()) == 2 * 4 ** d, d318 print(" count sets of degree at most d number 2^(d+1), so rules of degree at most d number 2 * 4^d for d = 1..8")319320def t5_fredkin():321 print("T5 fredkin = rule 150 tensor rule 150")322 mask = moore(2)323 side = 2 * HORIZON + 3324 mid = side // 2325 for birth, survive, label, centred in (326 ((1, 3, 5, 7), (0, 2, 4, 6, 8), "B1357/S02468", True),327 ((1, 3, 5, 7), (1, 3, 5, 7), "B1357/S1357", False),328 ):329 grid = np.zeros((side, side), dtype=np.uint8)330 grid[mid, mid] = 1331 row = np.zeros(side, dtype=np.uint8)332 row[mid] = 1333 pops = []334 for t in range(HORIZON + 1):335 outer = np.outer(row, row)336 if centred:337 assert np.array_equal(grid, outer), t338 pops.append(int(grid.sum()))339 if t in POWERS or t == 0:340 copies = np.zeros((side, side), dtype=np.uint8)341 for di in (-t, 0, t) if t else (0,):342 for dj in (-t, 0, t) if t else (0,):343 copies[mid + di, mid + dj] = 1344 if not centred and t:345 copies[mid, mid] = 0346 assert np.array_equal(grid, copies), (label, t)347 if not centred:348 shifted = outer.copy()349 shifted[mid, mid] ^= 1350 if t in POWERS:351 assert np.array_equal(grid, shifted), (label, t)352 grid = step(grid, birth, survive, mask)353 row = step150(row)354 squares = []355 line = np.zeros(side, dtype=np.uint8)356 line[mid] = 1357 for t in range(HORIZON + 1):358 squares.append(int(line.sum()) ** 2)359 line = step150(line)360 if centred:361 assert pops == squares362 print(f" {label} slices equal the outer product cell for cell to t = {HORIZON}")363 print(f" {label} population = rule 150 population squared, first 16: {pops[:16]}")364 print(f" {label} nine copies at t = {POWERS}")365 assert tuple(pops[:20]) == A246035366 else:367 for t in POWERS:368 assert pops[t] == squares[t] - 1, (label, t)369 print(f" {label} equals the outer product with the centre copy removed at every t = 2^j")370 print(f" {label} eight copies at t = {POWERS}, first 16 populations: {pops[:16]}")371 assert tuple(pops[:20]) == A160239372 line = np.zeros(side, dtype=np.uint8)373 line[mid] = 1374 row150 = []375 for _ in range(20):376 row150.append(int(line.sum()))377 line = step150(line)378 assert tuple(row150) == A071053379 print(f" rule 150 population, first 20: {row150}")380 print(" 20 terms each against OEIS A071053, its square A246035 for B1357/S02468, A160239 for B1357/S1357")381382def main():383 start = time.time()384 t1_moore()385 t2_life()386 pins = pin_grid()387 assert np.array_equal(pins, pin_brute())388 assert bool(pins[subset_code(LIFE[0]), subset_code(LIFE[1])]) == is_pin(rule_table(*LIFE))389 levels = level_set_grid()390 grid = degree_grid()391 hist = t3_census(pins, levels)392 t3_rules(grid, pins, levels, hist)393 t3_genus(pins, levels, grid)394 t4_affine(grid)395 t5_fredkin()396 print(f"domain: every one of the 2^18 life-like rules, all 2^18 axial checks against the closed form, {time.time() - start:.1f} s")397398if __name__ == "__main__":399 main()