collapse.py
6.6 kB · python · 236 lines
1import resource2import sys3import time4from fractions import Fraction5from itertools import product6from math import log78# CELLS910CELLS = [11 ("G", ((4, (0, 1)), (8, (0, 1, 2, 3)))),12 ("H", ((4, (0, 1)), (16, (0, 1, 4, 5)))),13 ("N", ((9, (0, 1, 2)), (27, tuple(range(9))))),14 ("T", ((4, (0, 1)), (8, (0, 1, 2, 3)), (16, tuple(range(8))))),15]1617HEIGHT = 10**131819SHARP = ((2, (1,)), (4, (0, 1, 2, 3)))2021IRRATIONAL = ((4, (0, 1, 2)), (16, tuple(range(16))))2223# FRAME2425def root_of(bases):26 for cand in range(2, min(bases) + 1):27 ok = True28 for b in bases:29 x = b30 while x % cand == 0:31 x //= cand32 if x != 1:33 ok = False34 break35 if ok:36 return cand37 return None3839def exponent_of(b, r):40 e = 041 while b % r == 0:42 b //= r43 e += 144 return e4546def lcm(values):47 out = 148 for v in values:49 g, w = out, v50 while w:51 g, w = w, g % w52 out = out * v // g53 return out5455def frame(cell):56 bases = [b for b, _ in cell]57 r = root_of(bases)58 exps = [exponent_of(b, r) for b in bases]59 return r, exps, lcm(exps)6061# BLOCKS6263def block_digits(cell):64 r, exps, m = frame(cell)65 keeps = [set(a) for _, a in cell]66 out = []67 for word in product(range(r), repeat=m):68 ok = True69 for keep, e in zip(keeps, exps):70 for g in range(m // e):71 v = 072 for t in range(e - 1, -1, -1):73 v = v * r + word[g * e + t]74 if v not in keep:75 ok = False76 break77 if not ok:78 break79 if ok:80 value = 081 for t in range(m - 1, -1, -1):82 value = value * r + word[t]83 out.append(value)84 return r, exps, m, tuple(sorted(out))8586def in_design(n, b, keep):87 while n:88 if n % b not in keep:89 return False90 n //= b91 return True9293def bottom_block(cell):94 r, _, m = frame(cell)95 keeps = [(b, set(a)) for b, a in cell]96 return tuple(97 c for c in range(r**m) if all(in_design(c, b, keep) for b, keep in keeps)98 )99100# DIMENSION101102def log_exact(value, r):103 t, x = 0, value104 while x % r == 0:105 x //= r106 t += 1107 return Fraction(t) if x == 1 else None108109def design_dim(b, allowed, r):110 t = log_exact(len(allowed), r)111 return None if t is None else t / exponent_of(b, r)112113def dim_text(count, r, m):114 t = log_exact(count, r)115 if t is None:116 return "log_%d(%d) / %d" % (r, count, m), log(count) / (m * log(r))117 return str(t / m), float(t / m)118119# ENUMERATION120121def walk(base, digits, limit):122 ds = sorted(digits)123 level = [d for d in ds if d and d < limit]124 while level:125 for v in level:126 yield v127 nxt = []128 for v in level:129 head = v * base130 if head >= limit:131 break132 for d in ds:133 w = head + d134 if w < limit:135 nxt.append(w)136 level = nxt137138def thin_side(cell, r):139 best = None140 for b, allowed in cell:141 d = design_dim(b, allowed, r)142 if best is None or d < best[0]:143 best = (d, b, allowed)144 return best[1], best[2]145146# VERBS147148def verb_blocks():149 for name, cell in CELLS:150 r, exps, m, blocks = block_digits(cell)151 assert blocks == bottom_block(cell), name152 shown = list(blocks) if len(blocks) <= 16 else list(blocks[:16]) + ["..."]153 print(154 "%s bases %s root %d exponents %s M %d card A %d A %s"155 % (name, [b for b, _ in cell], r, exps, m, len(blocks), shown)156 )157 r, exps, m = frame(SHARP)158 blocks = bottom_block(SHARP)159 keep = set(SHARP[0][1])160 bad = [n for n in walk(r**m, blocks, 10**6) if not in_design(n, SHARP[0][0], keep)]161 print(162 "sharpness bases %s exponents %s M %d A %s in the collapse and not in the joint set %s"163 % ([b for b, _ in SHARP], exps, m, list(blocks), bad[:4])164 )165166def verb_dim():167 for name, cell in CELLS + [("I", IRRATIONAL)]:168 r, exps, m, blocks = block_digits(cell)169 text, value = dim_text(len(blocks), r, m)170 dims = [design_dim(b, allowed, r) for b, allowed in cell]171 if None in dims:172 print(173 "%s bases %s root %d M %d card A %d exact dim %s = %.6f card A is not a power of the root, no rational part and no budget"174 % (name, [b for b, _ in cell], r, m, len(blocks), text, value)175 )176 continue177 raw = sum(dims) - (len(cell) - 1)178 budget = max(Fraction(0), raw)179 print(180 "%s exact dim %s = %.6f parts %s raw sum %s naive budget max(0, raw) %s = %.6f gap %s"181 % (182 name,183 text,184 value,185 " + ".join(str(d) for d in dims),186 raw,187 budget,188 float(budget),189 Fraction(text) - budget,190 )191 )192193def verb_check():194 for name, cell in CELLS:195 r, exps, m, blocks = block_digits(cell)196 big = r**m197 start = time.time()198 keeps = [(b, set(a)) for b, a in cell]199 inside = 0200 for n in walk(big, blocks, HEIGHT):201 assert all(in_design(n, b, keep) for b, keep in keeps), (name, n)202 inside += 1203 base, allowed = thin_side(cell, r)204 others = [(b, set(a)) for b, a in cell if b != base]205 marks = []206 power = big207 while power <= HEIGHT:208 marks.append(power)209 power *= big210 tally = [0] * len(marks)211 joint = 0212 for n in walk(base, allowed, HEIGHT):213 if all(in_design(n, b, keep) for b, keep in others):214 joint += 1215 for i, cut in enumerate(marks):216 if n < cut:217 tally[i] += 1218 law = [len(blocks) ** (i + 1) for i in range(len(marks))]219 seen = [t + 1 for t in tally]220 assert inside == joint, (name, inside, joint)221 assert seen == law, (name, seen, law)222 print(223 "%s height %d collapse inside the joint set %d joint count %d counts %s at %s %.1f s"224 % (name, HEIGHT, inside, joint, law, marks, time.time() - start)225 )226 peak = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / (1024 * 1024)227 print("peak resident %.0f MB" % peak)228229def main():230 verbs = sys.argv[1:] or ["blocks", "dim", "check"]231 table = {"blocks": verb_blocks, "dim": verb_dim, "check": verb_check}232 for verb in verbs:233 print("-- %s" % verb)234 table[verb]()235236main()