sibling.py
29.0 kB · python · 707 lines
1import sys2import time3from itertools import combinations, product4from math import gcd56import numpy as np78SEED = 17299OUT = []1011def say(line):12 print(line)13 OUT.append(line)1415def check(cond, what):16 if not cond:17 print("FAIL " + what)18 sys.exit(1)1920def corner_filled(code, coords):21 idx = 022 for c in coords:23 idx = 2 * idx + (c & 1)24 return (code >> idx) & 12526def tile(code, dim, side):27 grid = np.zeros((side,) * dim, dtype=np.uint8)28 for coords in product(range(side), repeat=dim):29 grid[coords] = corner_filled(code, coords)30 return grid3132def fractal(code, dim, side, level):33 base = tile(code, dim, side)34 out = base35 for _ in range(level - 1):36 out = np.kron(out, base)37 return out3839def mask_offsets(code, dim, side, level):40 grid = fractal(code, dim, side, level)41 centre = (side ** level - 1) // 242 cells = np.argwhere(grid == 1) - centre43 keep = np.any(cells != 0, axis=1)44 return [tuple(int(v) for v in row) for row in cells[keep]]4546# DICTIONARY4748def rule_table(number):49 return np.array([(number >> i) & 1 for i in range(8)], dtype=np.uint8)5051def outer_totalistic_elementary():52 found = []53 for number in range(256):54 t = rule_table(number)55 ok = True56 for c in range(2):57 for count in range(3):58 vals = {int(t[4 * l + 2 * c + r]) for l in range(2) for r in range(2) if l + r == count}59 if len(vals) > 1:60 ok = False61 if ok:62 found.append(number)63 return found6465def dictionary():66 ot = outer_totalistic_elementary()67 say(f"dictionary: outer-totalistic elementary rules {len(ot)}")68 check(len(ot) == 64, "64 outer-totalistic elementary rules")69 check(110 not in ot and 30 not in ot and 184 not in ot, "110, 30, 184 not outer-totalistic")70 check(90 in ot and 150 in ot and 204 in ot and 232 in ot, "90, 150, 204, 232 outer-totalistic")71 t = rule_table(110)72 say(f"dictionary: rule 110 on 001 -> {t[1]} and on 100 -> {t[4]}, so it is not outer-totalistic")73 say("dictionary: the 64 are " + " ".join(str(n) for n in ot))74 m1 = mask_offsets(1, 1, 3, 1)75 check(sorted(m1) == [(-1,), (1,)], "bang dim 1, code 1 at side 3 is the elementary mask")76 m7 = mask_offsets(7, 2, 3, 1)77 check(len(m7) == 8 and (0, 0) not in m7, "bang dim 2, code 7 at side 3 popped is Moore")78 say("dictionary: bang dim 1, code 1 at side 3 is [101]; bang dim 2, code 7 at side 3 popped is the 8-cell Moore mask")7980# TOWER8182def popcount_le_one_code(dim):83 code = 084 for i in range(2 ** dim):85 if bin(i).count("1") <= 1:86 code |= 1 << i87 return code8889def not_all_odd_code(dim):90 return (1 << (2 ** dim)) - 2 + 1 if dim == 0 else (1 << (2 ** dim)) - 1 - (1 << (2 ** dim - 1))9192def tower():93 codes_a = [popcount_le_one_code(d) for d in range(1, 5)]94 codes_b = [not_all_odd_code(d) for d in range(1, 5)]95 say(f"tower: popcount <= 1 codes {codes_a}; not all odd codes {codes_b}")96 check(codes_a == [3, 7, 23, 279], "popcount <= 1 codes")97 check(codes_b == [1, 7, 127, 32767], "not all odd codes")98 for d in range(1, 5):99 ma = mask_offsets(codes_a[d - 1], d, 3, 1)100 mb = mask_offsets(codes_b[d - 1], d, 3, 1)101 fa = int(tile(codes_a[d - 1], d, 3).sum())102 fb = int(tile(codes_b[d - 1], d, 3).sum())103 check(fa == 2 ** (d - 1) * (d + 2), "popcount <= 1 tile fill 2^(D-1)(D+2)")104 check(fb == 3 ** d - 1 and len(mb) == 3 ** d - 1, "not all odd tile is Moore")105 check(all(max(abs(v) for v in o) == 1 for o in mb), "Moore offsets within Chebyshev 1")106 same = sorted(ma) == sorted(mb)107 check(same == (d <= 2), "masks coincide exactly at D <= 2")108 say(f"tower: D={d} popcount<=1 tile {fa} mask {len(ma)}, not-all-odd tile {fb} mask {len(mb)}, masks equal {same}")109 menger = mask_offsets(23, 3, 3, 1)110 check(len(menger) == 20 and all(sum(1 for v in o if v == 0) <= 1 for o in menger), "Menger mask is 20 cells with at most one zero offset")111 cantor3 = mask_offsets(1, 1, 3, 3)112 check(sorted(abs(o[0]) for o in cantor3) == [5, 5, 7, 7, 11, 11, 13, 13], "Cantor level 3 is +-5 +-7 +-11 +-13")113 say(f"tower: Cantor level 3 offsets {sorted(o[0] for o in cantor3)}, Menger level 1 mask {len(menger)} cells")114115# DECOUPLING116117def hnf2(rows):118 rows = [list(r) for r in rows if r != (0, 0) and r != [0, 0]]119 while True:120 first = [r for r in rows if r[0] != 0]121 if len(first) <= 1:122 break123 first.sort(key=lambda r: abs(r[0]))124 p = first[0]125 new = [p]126 for r in first[1:]:127 q = r[0] // p[0]128 new.append([r[0] - q * p[0], r[1] - q * p[1]])129 rows = new + [r for r in rows if r[0] == 0]130 rows = [r for r in rows if r != [0, 0]]131 pivot = [r for r in rows if r[0] != 0]132 rest = [r for r in rows if r[0] == 0]133 g = 0134 for r in rest:135 g = gcd(g, abs(r[1]))136 return pivot, g137138def lattice_index(offsets, dim):139 if dim == 1:140 g = 0141 for (o,) in offsets:142 g = gcd(g, abs(o))143 return g if g else None144 if not offsets:145 return None146 arr = np.array(offsets, dtype=np.int64)147 v1 = arr[0]148 cross = v1[0] * arr[:, 1] - v1[1] * arr[:, 0]149 nz = np.flatnonzero(cross)150 if len(nz) == 0:151 return None152 basis = [tuple(int(x) for x in v1), tuple(int(x) for x in arr[nz[0]])]153 while True:154 pivot, g = hnf2(basis)155 a, b = pivot[0]156 if a < 0:157 a, b = -a, -b158 x = arr[:, 0]159 y = arr[:, 1]160 bad = (x % a != 0)161 q = x // a162 bad |= ((y - q * b) % g != 0)163 idx = np.flatnonzero(bad)164 if len(idx) == 0:165 return a * g166 basis = [(a, b), (0, g), tuple(int(x) for x in arr[idx[0]])]167168def eca_ring_step(state, offsets, table_bits, kind):169 if kind == "general":170 idx = state.astype(np.int64)171 for k, o in enumerate(offsets):172 idx |= np.roll(state, -o).astype(np.int64) << (k + 1)173 return table_bits[idx]174 count = np.zeros_like(state, dtype=np.int64)175 for o in offsets:176 count += np.roll(state, -o)177 birth, survive = table_bits178 return np.where(state == 1, survive[count], birth[count]).astype(np.uint8)179180def decoupling():181 rng = np.random.default_rng(SEED)182 seen = {}183 rows = []184 for dim in (1, 2):185 for code in range(2 ** (2 ** dim)):186 for side in (3, 5, 7, 9):187 for level in (1, 2, 3):188 offs = mask_offsets(code, dim, side, level)189 key = (dim, frozenset(offs))190 if key in seen:191 continue192 idx = lattice_index(offs, dim)193 seen[key] = idx194 rows.append((dim, code, side, level, len(offs), idx))195 hist = {}196 for dim, code, side, level, m, idx in rows:197 k = (dim, "rank<D" if idx is None else str(idx))198 hist[k] = hist.get(k, 0) + 1199 say(f"decoupling: distinct masks {len(rows)}; index histogram {sorted(hist.items())}")200 for dim in (1, 2):201 for code in range(2 ** (2 ** dim)):202 items = [f"n{side}L{level}:{m}:{'r' if idx is None else idx}" for d, c, side, level, m, idx in rows if d == dim and c == code and (dim == 1 or idx != 1)]203 if items:204 say(f"decoupling: D={dim} code={code} side.level:cells:index " + " ".join(items))205 one_d = {(side, level): lattice_index(mask_offsets(1, 1, side, level), 1) for side in (3, 5, 7, 9) for level in (1, 2, 3)}206 check(one_d[(3, 1)] == 1 and one_d[(3, 2)] == 2 and one_d[(3, 3)] == 1, "Cantor tower index 1, 2, 1 at levels 1..3")207 check(one_d[(5, 1)] == 2 and one_d[(9, 1)] == 2 and one_d[(7, 1)] == 1, "parity tile index by side")208 diag = lattice_index(mask_offsets(9, 2, 3, 1), 2)209 check(diag == 2, "diagonal 4-mask has index 2")210 vn = lattice_index(mask_offsets(6, 2, 3, 1), 2)211 check(vn == 1, "von Neumann 4-mask has index 1")212 tested = 0213 for dim, code, side, level, m, idx in rows:214 if dim != 1 or idx is None or idx == 1:215 continue216 offs = [o[0] for o in mask_offsets(code, dim, side, level)]217 small = [o // idx for o in offs]218 n = 12 * idx * 5219 for kind in ("general", "life"):220 if kind == "general" and m > 12:221 continue222 if kind == "general":223 table = rng.integers(0, 2, size=2 ** (m + 1)).astype(np.uint8)224 bits = table225 else:226 bits = (rng.integers(0, 2, size=m + 1).astype(np.uint8), rng.integers(0, 2, size=m + 1).astype(np.uint8))227 x = rng.integers(0, 2, size=n).astype(np.uint8)228 ys = [x[i::idx].copy() for i in range(idx)]229 for _ in range(40):230 x = eca_ring_step(x, offs, bits, kind)231 ys = [eca_ring_step(y, small, bits, kind) for y in ys]232 for i in range(idx):233 check(np.array_equal(x[i::idx], ys[i]), f"interleaving D=1 code={code} side={side} level={level} kind={kind}")234 tested += 1235 say(f"decoupling: interleaving equality checked on {tested} (mask, kind) pairs at D=1, 40 steps each")236237# COMPOSITES238239BLOCKS = np.arange(512, dtype=np.int64)240ROW = [(BLOCKS >> (6 - 3 * i)) & 7 for i in range(3)]241CENTRE = (BLOCKS >> 4) & 1242POP9 = np.array([int(b).bit_count() for b in range(512)], dtype=np.int64)243OUTER = POP9 - CENTRE244RULES = np.arange(256, dtype=np.int64)245TABLE = ((RULES[:, None] >> np.arange(8)[None, :]) & 1).astype(np.int64)246247def life_like_table(birth, survive):248 b = np.isin(OUTER, list(birth))249 s = np.isin(OUTER, list(survive))250 return np.where(CENTRE == 1, s, b).astype(np.uint8)251252def composites():253 idx = 4 * TABLE[:, ROW[0]] + 2 * TABLE[:, ROW[1]] + TABLE[:, ROW[2]]254 h = ((RULES[None, :, None] >> idx[:, None, :]) & 1).astype(np.uint8)255 flat = h.reshape(65536, 512)256 packed = np.packbits(flat, axis=1).view(np.uint64)257 distinct = np.unique(packed, axis=0).shape[0]258 say(f"composites: 65536 ordered pairs give {distinct} distinct 9-input rules")259 xor9 = life_like_table((1, 3, 5, 7), (0, 2, 4, 6, 8))260 check(np.array_equal(h[150, 150], xor9), "150 then 150 is the nine-cell XOR B1357/S02468")261 check(np.array_equal(h[105, 105], xor9), "105 then 105 is the nine-cell XOR too")262 for f in range(256):263 check(np.array_equal(h[f, 204], TABLE[f, ROW[1]].astype(np.uint8)), "f then 204 is f on the middle row")264 check(np.array_equal(h[f, 170], TABLE[f, ROW[2]].astype(np.uint8)), "f then 170 is f on the row below")265 check(np.array_equal(h[f, 240], TABLE[f, ROW[0]].astype(np.uint8)), "f then 240 is f on the row above")266 check(np.array_equal(h[f, 51], 1 - h[f, 204]) and np.array_equal(h[f, 85], 1 - h[f, 170]) and np.array_equal(h[f, 15], 1 - h[f, 240]), "the complementing three negate")267 say("composites: f then 204 / 170 / 240 is f on the middle / lower / upper row for all 256 f; then 51 / 85 / 15 is its negation")268 life = life_like_table((3,), (2, 3))269 hit = np.flatnonzero(np.all(flat == life[None, :], axis=1))270 check(len(hit) == 0, "Life is not a composite")271 say(f"composites: B3/S23 occurs among the 65536 composites {len(hit)} times")272 ot = np.ones(65536, dtype=bool)273 for c in range(2):274 for n in range(9):275 cols = np.flatnonzero((CENTRE == c) & (OUTER == n))276 sub = flat[:, cols]277 ot &= sub.min(axis=1) == sub.max(axis=1)278 say(f"composites: outer-totalistic composites {int(ot.sum())} ordered pairs")279 named = {}280 for k in np.flatnonzero(ot):281 f, g = divmod(int(k), 256)282 birth = tuple(n for n in range(9) if flat[k, np.flatnonzero((CENTRE == 0) & (OUTER == n))[0]])283 survive = tuple(n for n in range(9) if flat[k, np.flatnonzero((CENTRE == 1) & (OUTER == n))[0]])284 named.setdefault((birth, survive), []).append((f, g))285 say(f"composites: distinct life-like composites {len(named)}")286 for (birth, survive), pairs in sorted(named.items(), key=lambda kv: (len(kv[1]), kv[0])):287 bs = "B" + "".join(map(str, birth)) + "/S" + "".join(map(str, survive))288 shown = " ".join(f"{f}.{g}" for f, g in pairs[:12])289 say(f"composites: {bs} from {len(pairs)} pairs f.g: {shown}" + (" ..." if len(pairs) > 12 else ""))290 cell = [(i, j) for i in range(3) for j in range(3)]291 def perm_blocks(mapping):292 out = np.zeros(512, dtype=np.int64)293 for b in range(512):294 v = 0295 for (i, j) in cell:296 si, sj = mapping(i, j)297 bit = (b >> (8 - (3 * si + sj))) & 1298 v |= bit << (8 - (3 * i + j))299 out[b] = v300 return out301 syms = [lambda i, j: (j, i), lambda i, j: (i, 2 - j), lambda i, j: (2 - i, j), lambda i, j: (j, 2 - i)]302 transpose = perm_blocks(syms[0])303 tsym = np.all(flat == flat[:, transpose], axis=1)304 full = tsym.copy()305 for m in syms[1:]:306 pb = perm_blocks(m)307 full &= np.all(flat == flat[:, pb], axis=1)308 say(f"composites: transpose-symmetric {int(tsym.sum())} pairs, full dihedral symmetry {int(full.sum())} pairs")309 swapped = h.transpose(1, 0, 2).reshape(65536, 512)[:, transpose]310 same_order = np.all(flat == swapped, axis=1)311 say(f"composites: rows-first equals columns-first for {int(same_order.sum())} of 65536 ordered pairs")312 check(bool(ot[150 * 256 + 150]) and bool(full[150 * 256 + 150]), "150.150 is life-like and dihedral")313 return h314315# CANTOR LIFE316317CANTOR = (-13, -11, -7, -5, 5, 7, 11, 13)318NAMED = (319 ("B3/S23", (3,), (2, 3)),320 ("B36/S23", (3, 6), (2, 3)),321 ("B2/S", (2,), ()),322 ("B3/S012345678", (3,), tuple(range(9))),323 ("B3678/S34678", (3, 6, 7, 8), (3, 4, 6, 7, 8)),324 ("B1357/S02468", (1, 3, 5, 7), (0, 2, 4, 6, 8)),325)326RING = 1024327SOUP_STEPS = 2000328DENSITIES = (0.1, 0.25, 0.5)329SOUP_SEEDS = 5330WIDTH = 14331HORIZON = 256332REACH = 13333334def masks_of(birth, survive, m):335 b = np.zeros(m + 1, dtype=np.uint8)336 s = np.zeros(m + 1, dtype=np.uint8)337 b[list(birth)] = 1338 s[list(survive)] = 1339 return b, s340341def count_1d(x, offsets):342 c = np.zeros(x.shape, dtype=np.uint8)343 for o in offsets:344 c += np.roll(x, -o, axis=-1)345 return c346347def life_step(x, c, b, s):348 return np.where(x == 1, s[c], b[c]).astype(np.uint8)349350def fate_name(period, disp):351 if disp != 0:352 return "mover"353 return "still" if period == 1 else "oscillator"354355def cantor_soups():356 rng = np.random.default_rng(SEED)357 for name, birth, survive in NAMED:358 b, s = masks_of(birth, survive, 8)359 for dens in DENSITIES:360 x = (rng.random((SOUP_SEEDS, RING)) < dens).astype(np.uint8)361 seen = [dict() for _ in range(SOUP_SEEDS)]362 fates = ["undecided"] * SOUP_SEEDS363 for t in range(SOUP_STEPS + 1):364 for i in range(SOUP_SEEDS):365 if fates[i] != "undecided":366 continue367 if not x[i].any():368 fates[i] = "dies"369 continue370 key = x[i].tobytes()371 if key in seen[i]:372 p = t - seen[i][key]373 fates[i] = "fixed" if p == 1 else f"period {p}"374 seen[i][key] = t375 if t < SOUP_STEPS:376 x = life_step(x, count_1d(x, CANTOR), b, s)377 final = [int(r.sum()) for r in x]378 hist = np.bincount(count_1d(x, CANTOR).ravel(), minlength=9)379 say(f"cantor {name} density {dens}: fates {fates}; final live {final} of {RING}; count histogram {hist.tolist()}")380381def cantor_xor_period():382 kernel = {0} | set(CANTOR)383 def power_two(k):384 out = {}385 for o in kernel:386 r = (o * (1 << k)) % RING387 out[r] = out.get(r, 0) ^ 1388 return {r for r, v in out.items() if v}389 check(power_two(8) == {0}, "Cantor XOR kernel to the 256th power is 1 on the ring 1024")390 check(power_two(7) != {0}, "Cantor XOR kernel to the 128th power is not 1")391 say(f"cantor B1357/S02468: kernel^128 has support {sorted(power_two(7))} and kernel^256 = 1 on the ring 1024, so the period divides 256 and is 256 for a generic soup")392393def cantor_patterns():394 seeds = [1] + [(1 << (w - 1)) | 1 | (rest << 1) for w in range(2, WIDTH + 1) for rest in range(1 << (w - 2))]395 check(len(seeds) == 8192, "8192 seeds of width at most 14")396 b, s = masks_of((3,), (2, 3), 8)397 narrow = WIDTH + 2 * REACH * 16 + 2398 wide = WIDTH + 2 * REACH * HORIZON + 2399 x = np.zeros((len(seeds), narrow), dtype=np.uint8)400 start = (narrow - WIDTH) // 2401 for i, code in enumerate(seeds):402 for k in range(WIDTH):403 x[i, start + k] = (code >> k) & 1404 ids = np.arange(len(seeds))405 seen = [dict() for _ in seeds]406 fates = {}407 witness = {}408 shift = 0409 t = 0410 while t <= HORIZON and len(ids):411 if t == 16:412 pad = (wide - narrow) // 2413 x = np.pad(x, ((0, 0), (pad, pad)))414 shift = pad415 alive = x.any(axis=1)416 left = np.argmax(x, axis=1)417 right = x.shape[1] - 1 - np.argmax(x[:, ::-1], axis=1)418 keep = []419 for row, i in enumerate(ids):420 if not alive[row]:421 fates[i] = ("death", 0, 0, 0)422 continue423 key = x[row, left[row]:right[row] + 1].tobytes()424 store = seen[i]425 if key in store:426 t0, l0 = store[key]427 p = t - t0428 v = int(left[row]) - shift - l0429 fates[i] = (fate_name(p, v), p, v, int(x[row].sum()))430 witness[i] = key431 continue432 store[key] = (t, int(left[row]) - shift)433 keep.append(row)434 if t == HORIZON:435 for row in keep:436 fates[ids[row]] = ("undecided", 0, 0, int(x[row].sum()))437 break438 x = x[keep]439 ids = ids[keep]440 x = life_step(x, count_1d(x, CANTOR), b, s)441 t += 1442 tally = {}443 for i, (kind, p, v, n) in fates.items():444 tally[kind] = tally.get(kind, 0) + 1445 say(f"cantor B3/S23 patterns: {sorted(tally.items())} over {len(seeds)} seeds, horizon {HORIZON}")446 def smallest(kind):447 cands = [(n, p, i) for i, (k, p, v, n) in fates.items() if k == kind]448 if not cands:449 return None450 n, p, i = min(cands)451 cells = [k for k, c in enumerate(witness[i]) if c]452 return n, p, i, cells453 for kind in ("still", "oscillator"):454 w = smallest(kind)455 if w:456 n, p, i, cells = w457 say(f"cantor B3/S23 smallest {kind}: {n} cells, period {p}, cells {cells}, from seed {seeds[i]:b}")458 else:459 say(f"cantor B3/S23 smallest {kind}: not found")460 movers = [(i, p, v, n) for i, (k, p, v, n) in fates.items() if k == "mover"]461 for i, p, v, n in sorted(movers, key=lambda m: (m[3], m[1]))[:10]:462 cells = [k for k, c in enumerate(witness[i]) if c]463 say(f"cantor B3/S23 mover: {n} cells, period {p}, displacement {v}, cells {cells}, from seed {seeds[i]:b}")464 if not movers:465 say("cantor B3/S23 mover: not found over width <= 14")466 periods = sorted({p for k, p, v, n in fates.values() if k == "oscillator"})467 say(f"cantor B3/S23 oscillator periods found {periods}")468469# MENGER LIFE470471MENGER = tuple(o for o in product((-1, 0, 1), repeat=3) if sum(1 for v in o if v == 0) <= 1)472TORUS = 32473MENGER_STEPS = 200474MENGER_DENSITIES = (0.15, 0.3)475MENGER_SEEDS = 2476FIELD = 24477BOX = 5478MOVER_SEEDS = 200479MOVER_STEPS = 128480481def axis_sum(x, axis):482 return x + np.roll(x, 1, axis) + np.roll(x, -1, axis)483484def count_menger(x):485 box = axis_sum(axis_sum(axis_sum(x, -1), -2), -3)486 faces = sum(np.roll(x, d, axis) for axis in (-1, -2, -3) for d in (1, -1))487 return box - x - faces488489def menger_rules():490 subsets = [c for k in (1, 2) for c in combinations((3, 4, 5, 6), k)]491 return [(b, s) for b in subsets for s in subsets] + [((3,), (2, 3))]492493def rule_name(birth, survive):494 return "B" + "".join(map(str, birth)) + "/S" + "".join(map(str, survive))495496def menger_soups():497 rng = np.random.default_rng(SEED)498 rules = menger_rules()499 check(len(rules) == 101 and len(MENGER) == 20, "100 grid rules plus B3/S23 on the 20-cell mask")500 verdicts = {}501 for birth, survive in rules:502 b, s = masks_of(birth, survive, 20)503 x = np.zeros((len(MENGER_DENSITIES) * MENGER_SEEDS, TORUS, TORUS, TORUS), dtype=np.uint8)504 for j, dens in enumerate(MENGER_DENSITIES):505 for k in range(MENGER_SEEDS):506 x[j * MENGER_SEEDS + k] = rng.random((TORUS,) * 3) < dens507 seen = [dict() for _ in range(x.shape[0])]508 fates = ["active"] * x.shape[0]509 for t in range(MENGER_STEPS + 1):510 for i in range(x.shape[0]):511 if fates[i] != "active":512 continue513 if not x[i].any():514 fates[i] = "dies"515 continue516 key = x[i].tobytes()517 if key in seen[i]:518 p = t - seen[i][key]519 fates[i] = "fixed" if p == 1 else f"period {p}"520 seen[i][key] = t521 if t < MENGER_STEPS:522 x = life_step(x, count_menger(x), b, s)523 dens = [round(float(r.mean()), 3) for r in x]524 fates = ["explodes" if f == "active" and d > 0.4 else f for f, d in zip(fates, dens)]525 verdicts[(birth, survive)] = (fates, dens)526 tally = {}527 for fates, dens in verdicts.values():528 for f in fates:529 k = f.split()[0]530 tally[k] = tally.get(k, 0) + 1531 say(f"menger soups: 101 rules x 4 runs, run fates {sorted(tally.items())}")532 quiet = []533 active = []534 for (birth, survive), (fates, dens) in verdicts.items():535 is_quiet = all(f != "active" and f != "explodes" for f in fates)536 say(f"menger {rule_name(birth, survive)}: {fates} density {dens}" + (" QUIET" if is_quiet else ""))537 if is_quiet:538 quiet.append((birth, survive))539 elif all(f == "active" for f in fates):540 active.append(max(dens))541 say(f"menger soups: {len(quiet)} rules quiet at both densities; {len(active)} rules active in all four runs with final density between {min(active)} and {max(active)}; 0 runs above 0.4")542 return rules, quiet543544def bbox(x):545 axes = [np.flatnonzero(x.any(axis=tuple(a for a in range(3) if a != k))) for k in range(3)]546 return [int(a[0]) for a in axes], [int(a[-1]) for a in axes]547548def menger_movers(rules, quiet):549 rng = np.random.default_rng(SEED + 1)550 found = {}551 tallies = {}552 for birth, survive in rules:553 b, s = masks_of(birth, survive, 20)554 tally = {}555 for seed_no in range(MOVER_SEEDS):556 box = rng.integers(0, 2, size=(BOX,) * 3).astype(np.uint8)557 if not box.any():558 tally["death"] = tally.get("death", 0) + 1559 continue560 x = np.zeros((FIELD,) * 3, dtype=np.uint8)561 lo = (FIELD - BOX) // 2562 x[lo:lo + BOX, lo:lo + BOX, lo:lo + BOX] = box563 origin = np.zeros(3, dtype=np.int64)564 seen = {}565 kind = "undecided"566 for t in range(MOVER_STEPS + 1):567 if not x.any():568 kind = "death"569 break570 mn, mx = bbox(x)571 extent = max(hi - lo_ + 1 for lo_, hi in zip(mn, mx))572 if extent > FIELD - 2:573 kind = "growing"574 break575 shift = [FIELD // 2 - (lo_ + hi + 1) // 2 for lo_, hi in zip(mn, mx)]576 x = np.roll(x, shift, axis=(0, 1, 2))577 origin -= np.array(shift)578 mn = [m + sh for m, sh in zip(mn, shift)]579 mx = [m + sh for m, sh in zip(mx, shift)]580 crop = x[mn[0]:mx[0] + 1, mn[1]:mx[1] + 1, mn[2]:mx[2] + 1]581 key = (crop.shape, crop.tobytes())582 pos = tuple(int(m) for m in np.array(mn) + origin)583 if key in seen:584 t0, pos0 = seen[key]585 p = t - t0586 v = tuple(a - c for a, c in zip(pos, pos0))587 kind = fate_name(p, 1 if any(v) else 0)588 if kind == "mover":589 cells = [tuple(int(c) for c in cell) for cell in np.argwhere(crop == 1)]590 entry = (int(crop.sum()), p, v, cells, seed_no)591 if (birth, survive) not in found or entry < found[(birth, survive)]:592 found[(birth, survive)] = entry593 break594 seen[key] = (t, pos)595 if t < MOVER_STEPS:596 x = life_step(x, count_menger(x), b, s)597 tally[kind] = tally.get(kind, 0) + 1598 tallies[(birth, survive)] = tally599 if len(tally) > 1 or (birth, survive) in quiet or (birth, survive) == ((3,), (2, 3)):600 say(f"menger movers {rule_name(birth, survive)}: {sorted(tally.items())}" + (" QUIET" if (birth, survive) in quiet else ""))601 kinds = {}602 for tally in tallies.values():603 for k, v in tally.items():604 kinds[k] = kinds.get(k, 0) + v605 say(f"menger movers: seed fates over {len(rules)} rules {sorted(kinds.items())}")606 say(f"menger movers: {len(found)} of {len(rules)} rules carry a mover from {MOVER_SEEDS} seeds in a {BOX}^3 box, {sum(1 for r in found if r in quiet)} among the {len(quiet)} quiet rules")607 for (birth, survive), (n, p, v, cells, seed_no) in sorted(found.items(), key=lambda kv: kv[1][:2]):608 say(f"menger mover {rule_name(birth, survive)}: {n} cells, period {p}, displacement {v}, seed {seed_no}, cells {cells}")609610def menger_block():611 b, s = masks_of((3,), (2, 3), 20)612 depth = 32613 x = np.zeros((8, 8, depth), dtype=np.uint8)614 x[3:5, 3:5, depth // 2] = 1615 block = x[:, :, depth // 2].copy()616 row = np.zeros(depth, dtype=np.uint8)617 row[depth // 2] = 1618 for t in range(1, 9):619 x = life_step(x, count_menger(x), b, s)620 row = np.roll(row, 1) ^ np.roll(row, -1)621 profile = x.sum(axis=(0, 1))622 check(np.array_equal((profile == 4).astype(np.uint8), row) and np.all((profile == 0) | (profile == 4)), f"block stack is rule 90 at t={t}")623 for z in np.flatnonzero(row):624 check(np.array_equal(x[:, :, z], block), f"block shape kept at t={t}")625 say("menger B3/S23 block: a 2x2 plane block becomes a rule 90 stack of blocks along the normal, checked to t=8")626627# PRODUCTS628629TORUS2 = 256630PRODUCT_STEPS = 512631PRODUCT_G = (204, 170, 110, 54, 30, 90)632FIELD2 = 128633SEED_STEPS = 96634635def eca_axis(x, table, axis):636 l = np.roll(x, 1, axis)637 r = np.roll(x, -1, axis)638 return table[4 * l.astype(np.int64) + 2 * x + r]639640def composite_step(x, f_table, g_table):641 return eca_axis(eca_axis(x, f_table, -1), g_table, -2)642643def products_torus():644 rng = np.random.default_rng(SEED + 2)645 f_table = rule_table(110)646 for g in PRODUCT_G:647 g_table = rule_table(g)648 for seed_no in range(2):649 x = (rng.random((TORUS2, TORUS2)) < 0.3).astype(np.uint8)650 dens = [float(x.mean())]651 churn = []652 seen = {x.tobytes(): 0}653 period = 0654 for t in range(1, PRODUCT_STEPS + 1):655 nxt = composite_step(x, f_table, g_table)656 churn.append(float((nxt != x).mean()))657 x = nxt658 dens.append(float(x.mean()))659 key = x.tobytes()660 if key in seen and not period:661 period = t - seen[key]662 seen[key] = t663 say(f"product 110.{g} seed {seed_no}: density start {dens[0]:.3f} min {min(dens):.3f} max {max(dens):.3f} end {dens[-1]:.3f}, churn {np.mean(churn[-64:]):.3f}, periodic {'period ' + str(period) if period else 'no'}")664665def products_seeds():666 table = rule_table(110)667 n = 511668 x = np.zeros((n, FIELD2, FIELD2), dtype=np.uint8)669 at = FIELD2 - 16670 for code in range(1, 512):671 for i in range(3):672 for j in range(3):673 x[code - 1, at + i, at + j] = (code >> (8 - (3 * i + j))) & 1674 def corner(x):675 rows = x.any(axis=2)676 cols = x.any(axis=1)677 return np.argmax(rows, axis=1), np.argmax(cols, axis=1), FIELD2 - 1 - np.argmax(rows[:, ::-1], axis=1), FIELD2 - 1 - np.argmax(cols[:, ::-1], axis=1)678 y0, x0, y1, x1 = corner(x)679 pops = [x.sum(axis=(1, 2))]680 for t in range(1, SEED_STEPS + 1):681 x = composite_step(x, table, table)682 ya, xa, yb, xb = corner(x)683 check(np.all(ya == y0 - t) and np.all(xa == x0 - t) and np.all(yb == y1) and np.all(xb == x1), f"110.110 bounding box law at t={t}")684 pops.append(x.sum(axis=(1, 2)))685 final = pops[-1]686 say(f"product 110.110 seeds: all 511 grow, upper-left corner moves (-1,-1) per step and the lower-right corner is fixed, checked to t={SEED_STEPS}")687 say(f"product 110.110 seeds: population at t={SEED_STEPS} min {int(final.min())} median {int(np.median(final))} max {int(final.max())}; single cell {int(final[(1 << 4) - 1])}")688689def main():690 t0 = time.time()691 say(f"seed {SEED}")692 dictionary()693 tower()694 decoupling()695 composites()696 cantor_soups()697 cantor_xor_period()698 cantor_patterns()699 menger_block()700 rules, quiet = menger_soups()701 menger_movers(rules, quiet)702 products_torus()703 products_seeds()704 say(f"elapsed {time.time() - t0:.1f}s")705706if __name__ == "__main__":707 main()