verify.py

1# WALSH23from fractions import Fraction4from itertools import permutations56KRAW = [[1, 1, 1, 1], [3, 1, -1, -3], [3, -1, -1, 3], [1, -1, 1, -1]]78QUASI = {9    (0, 1): lambda n: (9 * n * n + 18 * n + 21) // 8,10    (0, 3): lambda n: (3 * n * n + 6 * n + 3) // 8,11    (1, 1): lambda n: (3 * n * n + 2 * n - 5) // 8,12    (1, 3): lambda n: (9 * n * n + 6 * n - 3) // 8,13    (2, 1): lambda n: (9 * n * n - 6 * n - 3) // 8,14    (2, 3): lambda n: (3 * n * n - 2 * n - 5) // 8,15    (3, 1): lambda n: (3 * n * n - 6 * n + 3) // 8,16    (3, 3): lambda n: (9 * n * n - 18 * n + 21) // 8,17}181920def corners():21    return [(a, b, c) for a in (0, 1) for b in (0, 1) for c in (0, 1)]222324def bit(p):25    return 4 * p[0] + 2 * p[1] + p[2]262728def design(code):29    return frozenset(p for p in corners() if code >> bit(p) & 1)303132def quasi(k, n):33    return QUASI[(k, n % 4)](n)343536def brute_counts(n):37    got = {e: 0 for e in corners()}38    for x in range(4 * n):39        for y in range(4 * n):40            z = 6 * n - 2 - x - y41            if z < 0 or z >= 4 * n or z % 2:42                continue43            got[((x // 4) % 2, (y // 4) % 2, (z // 4) % 2)] += 144    return got454647def macro_counts(n):48    got = {e: 0 for e in corners()}49    top = 6 * n - 250    for x in range(4 * n):51        lo = max(0, top - x - (4 * n - 1))52        hi = min(4 * n - 1, top - x)53        if lo > hi:54            continue55        ex = (x // 4) % 256        for c in range(8):57            if (c - x) % 2:58                continue59            first = lo + (c - lo) % 860            if first > hi:61                continue62            cnt = (hi - first) // 8 + 163            z = top - x - first64            got[(ex, (first // 4) % 2, (z // 4) % 2)] += cnt65    return got666768def solid_cut(n):69    top = 6 * n - 270    out = set()71    for x in range(4 * n):72        for y in range(4 * n):73            for z in range(0, 4 * n, 2):74                if x + y + z == top:75                    out.add((x, y, z))76    return out777879def plane_points(n):80    top = 6 * n - 281    out = set()82    for x in range(4 * n):83        for y in range(4 * n):84            z = top - x - y85            if 0 <= z < 4 * n and z % 2 == 0:86                out.add((x, y, z))87    return out888990def polymul(a, b, cap):91    r = [0] * (cap + 1)92    for i, ai in enumerate(a):93        if ai == 0 or i > cap:94            continue95        for j, bj in enumerate(b):96            if i + j > cap:97                break98            r[i + j] += ai * bj99    return r100101102def gf_extract(k, n):103    h = (n - 1) // 2104    cap = 3 * h + 1105    even = [0] * (cap + 1)106    for i in range(0, 2 * h + 1, 2):107        if i <= cap:108            even[i] = 1109    odd = [0] * (cap + 1)110    for i in range(1, 2 * h, 2):111        if i <= cap:112            odd[i] = 1113    p = [0] * (cap + 1)114    p[0] = 1115    if cap >= 1:116        p[1] = 6117    if cap >= 2:118        p[2] = 1119    for _ in range(3 - k):120        p = polymul(p, even, cap)121    for _ in range(k):122        p = polymul(p, odd, cap)123    return p[cap]124125126def box(m):127    if m >= 0 and m % 2 == 0:128        return (m + 2) * (m + 4) // 8129    return 0130131132def split_sum(k, n):133    h = (n - 1) // 2134    total = 0135    for a in (0, 1, 2):136        if (a - h - k - 1) % 2:137            continue138        rho = 6 if a == 1 else 1139        total += rho * (box(3 * h + 1 - k - a)140                        - (3 - k) * box(h - k - a - 1)141                        - k * box(h + 1 - k - a))142    return total143144145def spectrum(D):146    out = []147    for j in range(4):148        acc = 0149        for S in corners():150            if sum(S) != j:151                continue152            acc += sum((-1) ** sum(p[i] * S[i] for i in range(3)) for p in D)153        out.append(Fraction(acc, 8))154    return out155156157def spectrum_kraw(D):158    m = [0, 0, 0, 0]159    for p in D:160        m[sum(p)] += 1161    return [Fraction(sum(m[w] * KRAW[j][w] for w in range(4)), 8) for j in range(4)]162163164def ink_law(D, n):165    S = spectrum_kraw(D)166    s = -1 if n % 4 == 1 else 1167    return (S[0] - Fraction(1, 2) * S[3] * s168            + (Fraction(2, 3) * S[1] - Fraction(1, 3) * S[2] * s) / n169            + (Fraction(2, 3) * S[2]170               - (Fraction(1, 3) * S[1] + Fraction(1, 2) * S[3]) * s) / n ** 2)171172173def ink_exact(D, counts, n):174    return Fraction(sum(counts[p] for p in D), 6 * n * n)175176177def family_laws(n):178    s = -1 if n % 4 == 1 else 1179    carpet = (Fraction(1, 2) + Fraction(s, 8) + Fraction(1, 2 * n)180              - Fraction(s, 8 * n * n))181    tree = (Fraction(1, 4) + (Fraction(1, 3) - Fraction(s, 12)) / n182            + Fraction(1 - s, 6 * n * n))183    void = Fraction(1, 4) - Fraction(s, 4 * n) + Fraction(1, 2 * n * n)184    return {23: carpet, 232: 1 - carpet, 3: tree, 129: void}185186187def orbits():188    seen = {}189    reps = []190    for code in range(256):191        D = design(code)192        if D in seen:193            continue194        orbit = set()195        for perm in permutations(range(3)):196            for flip in corners():197                img = frozenset(198                    tuple(p[perm[i]] ^ flip[i] for i in range(3)) for p in D)199                orbit.add(img)200        for img in orbit:201            seen[img] = len(reps)202        reps.append((code, len(orbit)))203    return reps204205206def constituent(k, hpar, arg, raw=False):207    table = {(0, 0): (18, 18, 6), (0, 1): (6, 12, 6), (1, 0): (6, 4, 0),208             (1, 1): (18, 30, 12), (2, 0): (18, 6, 0), (2, 1): (6, 8, 2),209             (3, 0): (6, 0, 0), (3, 1): (18, 18, 6)}210    a, b, c = table[(k, hpar)]211    q = arg if raw else (arg - 1 - 2 * hpar) // 4212    return a * q * q + b * q + c213214215def main():216    small = list(range(1, 56, 2))217    for n in small:218        got = brute_counts(n)219        for e in corners():220            want = quasi(sum(e), n)221            assert got[e] == want, "counts n=%d eps=%s got %d want %d" % (222                n, e, got[e], want)223        total = sum(got.values())224        assert total == 6 * n * n, "slice size n=%d got %d want %d" % (225            n, total, 6 * n * n)226        for k in range(4):227            vals = {got[e] for e in corners() if sum(e) == k}228            assert len(vals) == 1, "weight split n=%d k=%d got %s" % (229                n, k, sorted(vals))230            g = gf_extract(k, n)231            assert g == quasi(k, n), "gf n=%d k=%d got %d want %d" % (232                n, k, g, quasi(k, n))233            f = split_sum(k, n)234            assert f == quasi(k, n), "finite sum n=%d k=%d got %d want %d" % (235                n, k, f, quasi(k, n))236        laws = family_laws(n)237        for code, want in laws.items():238            gotink = ink_law(design(code), n)239            assert gotink == want, "family n=%d code=%d got %s want %s" % (240                n, code, gotink, want)241        for code in range(256):242            D = design(code)243            gotink = ink_exact(D, got, n)244            want = ink_law(D, n)245            assert gotink == want, "ink n=%d code=%d got %s want %s" % (246                n, code, gotink, want)247        print("odd n = %d: eight counts, three routes, 256 inks all agree" % n)248249    for n in (101, 555, 999, 9991):250        got = macro_counts(n)251        for e in corners():252            want = quasi(sum(e), n)253            assert got[e] == want, "cold n=%d eps=%s got %d want %d" % (254                n, e, got[e], want)255        total = sum(got.values())256        assert total == 6 * n * n, "cold size n=%d got %d want %d" % (257            n, total, 6 * n * n)258        for code in range(256):259            D = design(code)260            gotink = ink_exact(D, got, n)261            want = ink_law(D, n)262            assert gotink == want, "cold ink n=%d code=%d got %s want %s" % (263                n, code, gotink, want)264        print("cold n = %d: interval count matches, 256 inks exact" % n)265266    for n in (1, 3, 5, 7, 9):267        cut = solid_cut(n)268        pts = plane_points(n)269        assert cut == pts, "solid cut n=%d got %d cells want %d" % (270            n, len(cut), len(pts))271        assert len(cut) == 6 * n * n, "solid cut size n=%d got %d want %d" % (272            n, len(cut), 6 * n * n)273    print("solid cut: n = 1, 3, 5, 7, 9 agree with the slice cell for cell")274275    nine = tuple(quasi(k, 9) for k in range(4))276    assert nine == (114, 32, 84, 24), "n=9 got %s want (114, 32, 84, 24)" % (277        nine,)278    assert tuple(gf_extract(k, 9) for k in range(4)) == nine, \279        "n=9 gf got %s want %s" % (tuple(gf_extract(k, 9) for k in range(4)), nine)280    print("n = 9 anchor: got %s want (114, 32, 84, 24)" % (nine,))281282    one = tuple(quasi(k, 1) for k in range(4))283    assert one == (6, 0, 0, 0), "n=1 got %s want (6, 0, 0, 0)" % (one,)284    print("n = 1 boundary: got %s want (6, 0, 0, 0)" % (one,))285286    ladder = set()287    for code in range(256):288        D = design(code)289        a = spectrum(D)290        b = spectrum_kraw(D)291        assert a == b, "krawtchouk code=%d got %s want %s" % (code, b, a)292        blink = -a[3] / 2293        assert (blink == 0) == (a[3] == 0), "blink zero code=%d" % code294        ladder.add(abs(blink))295    want = {Fraction(0), Fraction(1, 16), Fraction(1, 8), Fraction(3, 16),296            Fraction(1, 4)}297    assert ladder == want, "ladder got %s want %s" % (298        sorted(ladder), sorted(want))299    print("blink ladder: got %s want the same five rungs" % (300        sorted(str(v) for v in ladder),))301302    reps = orbits()303    assert len(reps) == 22, "classes got %d want 22" % len(reps)304    assert sum(size for _, size in reps) == 256, "orbit sizes got %d want 256" % \305        sum(size for _, size in reps)306    still = sum(1 for code, _ in reps if spectrum_kraw(design(code))[3] == 0)307    assert still == 9, "stationary classes got %d want 9" % still308    print("symmetry: got 22 classes covering 256 designs, 9 stationary")309310    for q in range(0, 12):311        for k in range(4):312            assert constituent(k, 0, 4 * q + 1) == quasi(k, 4 * q + 1), \313                "constituent k=%d n=%d" % (k, 4 * q + 1)314            assert constituent(k, 1, 4 * q + 3) == quasi(k, 4 * q + 3), \315                "constituent k=%d n=%d" % (k, 4 * q + 3)316            left = constituent(k, 0, -q - 1, True)317            right = constituent(3 - k, 1, q, True)318            assert left == right, "reciprocity q=%d k=%d got %d want %d" % (319                q, k, left, right)320    print("reciprocity: q -> -q-1 sends constituent k to constituent 3-k")321322    hexa = [18 * q * q + 18 * q + 6 for q in range(6)]323    assert hexa == [6, 42, 114, 222, 366, 546], "A164016 got %s" % hexa324    assert all(quasi(0, 4 * q + 1) == hexa[q] for q in range(6)), "A164016 k=0"325    assert all(quasi(3, 4 * q + 3) == hexa[q] for q in range(6)), "A164016 k=3"326    print("A164016: got %s want 6, 42, 114, 222, 366, 546" % hexa)327328    print("all green")329330331if __name__ == "__main__":332    main()