verify.py

1import math2from fractions import Fraction34# ALPHABET56CODES = tuple(range(1, 16))7UNIT = (1, 2, 4, 8)8ROW = (3, 12)9COLUMN = (5, 10)10DIAGONAL = (6, 9)11GASKET = (7, 11, 13, 14)12FULL = 1513FREE = UNIT + DIAGONAL14DOMINO = ROW + COLUMN15HEAVY = GASKET + (FULL,)16PAIRS = tuple((a, b) for i, a in enumerate(CODES) for b in CODES[i + 1:])1718def fill(code):19    return bin(code).count("1")2021def word_fill(word):22    out = 123    for c in word:24        out *= fill(c)25    return out2627def klass(code):28    if code in UNIT:29        return "unit"30    if code in ROW:31        return "row domino"32    if code in COLUMN:33        return "column domino"34    if code in DIAGONAL:35        return "diagonal"36    if code in GASKET:37        return "gasket"38    return "full"3940# GEOMETRY4142def tile(code):43    return [(code & 1) | ((code >> 1 & 1) << 1), (code >> 2 & 1) | ((code >> 3 & 1) << 1)]4445def grid(word):46    rows, wide = [1], 147    for code in word:48        small = tile(code)49        out = []50        for upper in rows:51            for lower in small:52                value = 053                for j in range(wide):54                    if (upper >> j) & 1:55                        value |= lower << (2 * j)56                out.append(value)57        rows, wide = out, wide * 258    return rows, wide5960def segments(row):61    out = []62    rest = row63    while rest:64        low = (rest & -rest).bit_length() - 165        tail = rest >> low66        run = ((~tail) & (tail + 1)).bit_length() - 167        out.append((low, low + run - 1))68        rest &= ~(((1 << run) - 1) << low)69    return out7071def components(rows):72    parent = []7374    def find(x):75        while parent[x] != x:76            parent[x] = parent[parent[x]]77            x = parent[x]78        return x7980    prev, prev_ids = [], []81    for row in rows:82        segs = segments(row)83        ids = []84        for _ in segs:85            ids.append(len(parent))86            parent.append(len(parent))87        for k, (low, high) in enumerate(segs):88            for m, (plow, phigh) in enumerate(prev):89                if plow <= high and low <= phigh:90                    ra, rb = find(prev_ids[m]), find(ids[k])91                    if ra != rb:92                        parent[rb] = ra93        prev, prev_ids = segs, ids94    return len({find(i) for i in range(len(parent))})9596def drawn(word):97    return components(grid(word)[0])9899def contacts(rows, wide):100    top, bottom = rows[0], rows[-1]101    h = sum(1 for r in rows if (r & 1) and (r >> (wide - 1)) & 1)102    v = sum(1 for j in range(wide) if (top >> j) & 1 and (bottom >> j) & 1)103    return h, v104105# REGIMES106107def regime(a, b):108    if a in FREE or b in FREE:109        return "cut"110    if a in DOMINO and b in DOMINO:111        if (a in ROW) == (b in ROW):112            return "one"113        return "crossed"114    if a in HEAVY and b in HEAVY:115        return "one"116    heavy = a if a in HEAVY else b117    return "block" if heavy == FULL else "runs"118119def closed(a, b, word):120    kind = regime(a, b)121    if kind == "one":122        return 1123    if kind == "cut":124        spots = [i for i, c in enumerate(word) if c in FREE]125        return word_fill(word[: spots[-1] + 1]) if spots else 1126    if kind == "crossed":127        last = word[-1]128        run = 0129        while run < len(word) and word[len(word) - 1 - run] == last:130            run += 1131        return 1 << (len(word) - run)132    if kind == "block":133        n = sum(1 for c in word if c == FULL)134        run = 0135        while run < len(word) and word[len(word) - 1 - run] == FULL:136            run += 1137        return 1 << (n - run)138    heavy = a if a in GASKET else b139    light = b if a in GASKET else a140    spots = [i for i, c in enumerate(word) if c == light]141    if not spots:142        return 1143    stop = spots[-1]144    total = 1145    for i, c in enumerate(word[: stop + 1]):146        if c == heavy:147            total += word_fill(word[:i])148    return total149150# WEIGHTS151152def fill_weight(code):153    return {1: (0, 0), 2: (1, 0), 3: (0, 1), 4: (2, 0)}[fill(code)]154155def phi_weight(code):156    return (1, 0) if code in DIAGONAL else (0, 0)157158def rate_weight(a, b):159    kind = regime(a, b)160    if kind == "one":161        return (0, 0), (0, 0)162    if kind == "cut":163        return fill_weight(a), fill_weight(b)164    if kind == "crossed":165        return (1, 0), (1, 0)166    if kind == "block":167        return ((0, 0), (1, 0)) if a in DOMINO else ((1, 0), (0, 0))168    return (fill_weight(a), fill_weight(b))169170def nats(weight):171    return weight[0] * math.log(2) + weight[1] * math.log(3)172173# TRACK174175def log2_int(n):176    shift = n.bit_length() - 53177    return math.log2(n >> shift) + shift if shift > 0 else math.log2(n)178179LOG2_3 = math.log2(3)180LOG2_FILL = {1: 0.0, 2: 1.0, 3: LOG2_3, 4: 2.0}181182def free_place(word):183    for i in range(len(word) - 1, -1, -1):184        if word[i] in FREE:185            return i + 1186    return 0187188def gd_exact(word, gasket, wanted):189    out = {}190    fill_now, total, snap = 1, 1, 1191    for step, c in enumerate(word, start=1):192        if c == gasket:193            total += fill_now194            fill_now *= 3195        else:196            snap = total197            fill_now *= 2198        if step in wanted:199            out[step] = (snap, fill_now)200    return out201202def rates(a, b, word):203    kind = regime(a, b)204    gasket = a if a in GASKET else b205    logfill, top, ratio = 0.0, 0.0, 1.0206    snap_top, snap_ratio = 0.0, 1.0207    cut_top = 0.0208    run, fulls, full_run = 0, 0, 0209    last = None210    out = []211    for step, c in enumerate(word, start=1):212        weight = LOG2_FILL[fill(c)]213        run = run + 1 if c == last else 1214        last = c215        if c == FULL:216            fulls += 1217            full_run += 1218        else:219            full_run = 0220        if kind == "runs":221            if c == gasket:222                ratio = 1.0 + ratio * (2.0 ** (top - logfill))223                top = logfill224            else:225                snap_top, snap_ratio = top, ratio226            out.append((snap_top + math.log2(snap_ratio)) / step)227        elif kind == "cut":228            if c in FREE:229                cut_top = logfill + weight230            out.append(cut_top / step)231        elif kind == "one":232            out.append(0.0)233        elif kind == "crossed":234            out.append((step - run) / step)235        else:236            out.append((fulls - full_run) / step)237        logfill += weight238    return out239240# CHECKS241242def die(message):243    raise SystemExit("verify.py: " + message)244245def check(condition, message):246    if not condition:247        die(message)248249def pair_word(a, b, mask, length):250    return [b if (mask >> i) & 1 else a for i in range(length)]251252def check_contacts():253    seen = 0254    for length in range(1, 5):255        for mask in range(15 ** length):256            word, rest = [], mask257            for _ in range(length):258                word.append(CODES[rest % 15])259                rest //= 15260            rows, wide = grid(word)261            h, v = contacts(rows, wide)262            wanted_h, wanted_v = 1, 1263            for c in word:264                ch, cv = contacts(*grid([c]))265                wanted_h *= ch266                wanted_v *= cv267            check((h, v) == (wanted_h, wanted_v), f"contacts fail on {word}")268            seen += 1269    print(f"contacts multiply on all {seen} words of length at most 4 over the 15 codes")270271def check_forms(reach):272    seen = 0273    for a, b in PAIRS:274        for length in range(1, reach + 1):275            for mask in range(1 << length):276                word = pair_word(a, b, mask, length)277                want = closed(a, b, word)278                got = drawn(word)279                check(got == want, f"{(a, b)} word {word}: drew {got}, form says {want}")280                seen += 1281    print(f"closed forms match drawn cells on all {seen} words of length at most {reach} over all {len(PAIRS)} pairs")282283def check_forms_deep(pairs, reach):284    seen = 0285    for a, b in pairs:286        for length in range(1, reach + 1):287            for mask in range(1 << length):288                word = pair_word(a, b, mask, length)289                check(drawn(word) == closed(a, b, word), f"{(a, b)} deep word {word}")290                seen += 1291    print(f"closed forms match drawn cells on all {seen} words of length at most {reach} over {len(pairs)} named pairs")292293def check_census():294    tally = {}295    for a, b in PAIRS:296        tally[regime(a, b)] = tally.get(regime(a, b), 0) + 1297    check(tally == {"cut": 69, "one": 12, "crossed": 4, "block": 4, "runs": 16}, f"regime census is {tally}")298    print("regime census: cut 69, one 12, crossed 4, block 4, runs 16, total 105")299300def check_rates():301    for a, b in PAIRS:302        wa, wb = rate_weight(a, b)303        for share in (3, 5):304            length = 4000305            word = [a if (i * share) % 8 < 4 else b for i in range(length)]306            fa = sum(1 for c in word if c == a) / length307            got = rates(a, b, word)[-1] * math.log(2)308            want = fa * nats(wa) + (1 - fa) * nats(wb)309            check(abs(got - want) < 0.02, f"{(a, b)} rate {got} against {want}")310    print(f"the closed-form rate matches the weight table on all {len(PAIRS)} pairs at two frequencies, length 4000")311312def check_ledger():313    saturating = short = exact = refuted = between = 0314    smooth = 0315    for a, b in PAIRS:316        wa, wb = rate_weight(a, b)317        if (wa, wb) == (fill_weight(a), fill_weight(b)):318            saturating += 1319        else:320            short += 1321        if (wa, wb) == (phi_weight(a), phi_weight(b)):322            exact += 1323        else:324            refuted += 1325        middle = (nats(wa) + nats(wb)) / 2326        low = (nats(phi_weight(a)) + nats(phi_weight(b))) / 2327        high = (nats(fill_weight(a)) + nats(fill_weight(b))) / 2328        if low + 1e-12 < middle < high - 1e-12:329            between += 1330            check(regime(a, b) == "block", f"{(a, b)} is strictly between and is not a domino against the full tile")331        if regime(a, b) != "runs":332            smooth += 1333    check((saturating, short) == (89, 16), f"saturation ledger is {(saturating, short)}")334    check((exact, refuted) == (27, 78), f"Phi ledger is {(exact, refuted)}")335    check(between == 4, f"{between} pairs strictly between")336    check(smooth == 89, f"{smooth} pairs outside the gasket-and-domino regime")337    print("ledger: 89 saturate and 16 fall short, Phi exact on 27 and refuted on 78, 4 strictly between")338339def check_smooth(reach):340    biggest = 0341    for a, b in PAIRS:342        for length in range(1, reach + 1):343            for mask in range(1 << length):344                word = pair_word(a, b, mask, length)345                value = closed(a, b, word)346                if regime(a, b) == "runs":347                    if (a, b) == (3, 7) and length == 8:348                        biggest = max(biggest, value)349                    continue350                rest = value351                for prime in (2, 3):352                    while rest % prime == 0:353                        rest //= prime354                check(rest == 1, f"{(a, b)} word {word} has count {value}, not 3-smooth")355    check(biggest == 1094 and 1094 == 2 * 547, f"largest count at length 8 over (3,7) is {biggest}")356    print(f"every count on the 89 pairs outside the gasket-and-domino regime is 3-smooth to length {reach}; the largest count at length 8 over (3,7) is 1094 = 2 x 547")357358# MORSE359360def morse(length):361    return [bin(i).count("1") & 1 for i in range(length)]362363def gasket_domino_track(word, gasket):364    fill_log, total_log, total_ratio = 0.0, 0.0, 1.0365    snap_log, snap_ratio = 0.0, 1.0366    out = []367    for c in word:368        if c == gasket:369            shift = total_log - fill_log370            total_ratio = 1.0 + total_ratio * (2.0 ** shift)371            total_log = fill_log372            fill_log += LOG2_3373        else:374            snap_log, snap_ratio = total_log, total_ratio375            fill_log += 1.0376        out.append((snap_log + math.log2(snap_ratio), fill_log))377    return out378379def check_morse():380    reach = 1 << 14381    bits = morse(reach)382    certificate = math.log(108) + 0.5 * math.log(1.5)383    worst = 0.0384    printed = {}385    for swap in (0, 1):386        word = [7 if (bit ^ swap) == 0 else 3 for bit in bits]387        track = gasket_domino_track(word, 7)388        exact = gd_exact(word, 7, {4096, reach})389        for length in range(4, reach + 1):390            logcomp = track[length - 1][0] * math.log(2)391            worst = max(worst, abs(logcomp - (length / 2) * math.log(6)))392        printed[swap] = (track[4095][0] / 4096, track[reach - 1][0] / reach)393        for cut in (4096, reach):394            check(abs(track[cut - 1][0] - log2_int(exact[cut][0])) < 1e-6, "the float track parts from the exact count")395    check(worst < certificate, f"largest deviation {worst} against the certificate {certificate}")396    check(abs(worst - 4.273459) < 5e-7, f"largest deviation is {worst}")397    check(abs(certificate - 4.884864) < 5e-7, f"certificate constant is {certificate}")398    check(abs(printed[1][0] - 1.291967463826) < 5e-13 and abs(printed[1][1] - 1.292352803727) < 5e-13, f"reading one prints {printed[1]}")399    check(abs(printed[0][0] - 1.291291597694) < 5e-13 and abs(printed[0][1] - 1.292183837194) < 5e-13, f"reading two prints {printed[0]}")400    limit = math.log2(6) / 2401    check(abs(limit - 1.292481250360578) < 5e-16, "the limit in log 2 units")402    check(abs(math.log(6) / 2 - 0.895879734614027) < 5e-16, "the limit in nats")403    print(f"Thue-Morse: every length from 4 to 2^14 in both readings obeys the certificate, largest deviation {worst:.6f} against {certificate:.6f}")404    print(f"Thue-Morse prefix rates in log 2 units: {printed[1][0]:.12f} and {printed[1][1]:.12f} in one reading, {printed[0][0]:.12f} and {printed[0][1]:.12f} in the other, against (1/2) log_2 6 = {limit:.12f}")405406def check_saturation():407    reach = 1 << 14408    bits = morse(reach)409    for swap, want_max in ((1, Fraction(43397, 186624)), (0, Fraction(151, 648))):410        word = [7 if (bit ^ swap) == 0 else 3 for bit in bits]411        track = gasket_domino_track(word, 7)412        ratios = [logc - logf for logc, logf in track]413        low = min(ratios)414        best = max(range(4, reach), key=lambda i: ratios[i])415        exact = gd_exact(word, 7, {best + 1, 4096})416        got = Fraction(*exact[best + 1])417        check(got == want_max, f"largest saturation at length >= 5 in reading {swap} is {got}")418        check(abs(2.0 ** low - (0.0113766545 if swap == 1 else 0.0113766545)) < 5e-10, f"minimum saturation {2.0 ** low}")419        check(2.0 ** low > 1 / 108, "the proved floor 1/108 is broken")420        if swap == 1:421            at = Fraction(*exact[4096])422            check(abs(float(at) - 0.2325367033) < 5e-11, f"saturation at length 4096 is {float(at)}")423    print("Thue-Morse saturation: minimum 0.0113766545 above the floor 1/108, exact maxima 43397/186624 and 151/648, and 0.2325367033 at length 4096")424425# BOUNDARY426427def check_boundary():428    squares = [6 if int(math.isqrt(i)) ** 2 == i else 3 for i in range(1, 1 << 12)]429    for n in range(2, 40):430        length = n * n431        check(free_place(squares[:length]) == length, "the squares word misses rate 1")432        edge = (n + 1) ** 2 - 1433        check(free_place(squares[:edge]) == n * n, "the squares word misses n/(n+2)")434        check(closed(3, 6, squares[:edge]) == 1 << (n * n), "the squares word parts from the closed form")435    powers = [6 if (i & (i - 1)) == 0 else 3 for i in range(1, 1 << 15)]436    seen = []437    for k in range(2, 14):438        length = (1 << (k + 1)) - 1439        seen.append(Fraction(free_place(powers[:length]), length))440    check(seen[0] == Fraction(4, 7) and seen[1] == Fraction(8, 15) and seen[-1] == Fraction(1 << 13, (1 << 14) - 1), f"powers-of-2 rates start {seen[:2]}")441    check(all(f > Fraction(1, 2) for f in seen), "the lower rate dips below one half")442    check(all(free_place(powers[: 1 << k]) == 1 << k for k in range(1, 14)), "the powers-of-2 word misses rate 1")443    print("boundary over (3,6): the squares word has rate 1 at every square and n/(n+2) just before the next, the powers-of-2 word reads 4/7, 8/15, ..., 8192/16383")444    reach = 1 << 15445    word = [7 if (i & (i - 1)) == 0 else 3 for i in range(1, reach + 1)]446    track = gasket_domino_track(word, 7)447    at = {length: track[length - 1][0] / length for length in (2048, 2049, 32768)}448    check(abs(at[2048] - 0.502367981) < 5e-10, f"L = 2048 reads {at[2048]}")449    check(abs(at[2049] - 1.002164269) < 5e-10, f"L = 2049 reads {at[2049]}")450    check(abs(at[32768] - 0.500219405) < 5e-10, f"L = 32768 reads {at[32768]}")451    block = [track[length - 1][0] / length for length in range(4097, 8193)]452    check(abs(block[0] - 1.00123) < 5e-6 and abs(block[-1] - 0.50073) < 5e-6, f"block runs {block[0]} to {block[-1]}")453    check(abs(max(block) - block[0]) < 1e-12 and abs(min(block) - block[-1]) < 1e-12, "the block is not monotone")454    print("boundary over (3,7): the powers-of-2 word reads 0.502367981 at 2048, 1.002164269 at 2049, 0.500219405 at 32768, and sweeps 1.00123 down to 0.50073 over one block")455    triple = [7, 3]456    while len(triple) < 4096:457        n = len(triple)458        triple = triple + [7] * n + [3] * n459    track = gasket_domino_track(triple[:4096], 7)460    band = [track[length - 1][0] / length for length in range(1024, 4097)]461    check(abs(min(band) - 0.4792) < 5e-5 and abs(max(band) - 1.4379) < 5e-5, f"tripling band is [{min(band)}, {max(band)}]")462    print(f"the tripling word keeps both letters at lower density 1/4 and its prefix rate still ranges over [{min(band):.4f}, {max(band):.4f}] on 1024 <= L <= 4096")463464# COCYCLE465466def check_cocycle():467    left, right = ((0, 1), (-2, 3)), ((2, 0), (4, 0))468    matrix = {3: left, 6: right}469470    def apply(vector, m):471        return (vector[0] * m[0][0] + vector[1] * m[1][0], vector[0] * m[0][1] + vector[1] * m[1][1])472473    seen = 0474    for length in range(1, 13):475        for mask in range(1 << length):476            word = pair_word(3, 6, mask, length)477            vector = (1, 0)478            for c in word:479                vector = apply(vector, matrix[c])480            check(vector[0] + vector[1] == closed(3, 6, word), f"the cocycle misses {word}")481            seen += 1482    power = ((1, 0), (0, 1))483    for step in range(1, 25):484        power = (485            (power[0][0] * left[0][0] + power[0][1] * left[1][0], power[0][0] * left[0][1] + power[0][1] * left[1][1]),486            (power[1][0] * left[0][0] + power[1][1] * left[1][0], power[1][0] * left[0][1] + power[1][1] * left[1][1]),487        )488        biggest = max(abs(power[0][0]), abs(power[0][1]), abs(power[1][0]), abs(power[1][1]))489        check(biggest == (1 << (step + 1)) - 1, f"the largest entry at step {step} is {biggest}")490        check(closed(3, 6, [3] * step) == 1, "the constant domino word is not connected")491    print(f"the two-by-two cocycle on (3,6) reproduces the count on all {seen} words to length 12, its largest entry is 2^(L+1) - 1, and the constant word stays at one component")492493def check_byproduct():494    for k in range(1, 9):495        word = [7, 3] * k496        want = (6 ** k + 4) // 5497        check(closed(7, 3, word) == want, f"(7,3)^{k} is not (6^k + 4)/5")498        if k <= 4:499            check(drawn(word) == want, f"(7,3)^{k} drawn")500    print("the stationary control comp((7,3)^k) = (6^k + 4)/5 reads 2, 8, 44, 260, 1556, 9332 to k = 8, drawn to k = 4")501502def main():503    check_contacts()504    check_census()505    check_forms(7)506    check_forms_deep(((3, 7), (5, 7), (3, 15), (3, 6), (3, 5), (7, 15), (1, 3), (6, 9)), 8)507    check_rates()508    check_ledger()509    check_smooth(8)510    check_morse()511    check_saturation()512    check_boundary()513    check_cocycle()514    check_byproduct()515    print("all checks green")516517if __name__ == "__main__":518    main()