anf.py

1import time2from itertools import product3from math import gcd45import numpy as np67BASE2_DIMS = (2, 3, 4)8IDENTITY_CASES = ((2, 2), (2, 3), (2, 4), (3, 2), (3, 3), (5, 2), (7, 2))9COMPOSITE_BASES = (4, 6, 8, 9, 10, 12)10PRIME_BASES = (2, 3, 5, 7, 11, 13)11CLAIMED_TABLE = {2: (2, 2, 3, 4), 3: (4, 4, 6, 6)}1213def cells(q, d):14    return list(product(range(q), repeat=d))1516def strides(q, d):17    return [q ** (d - 1 - a) for a in range(d)]1819def mod_inverse(a, q):20    old_r, r = a % q, q21    old_s, s = 1, 022    while r:23        quotient = old_r // r24        old_r, r = r, old_r - quotient * r25        old_s, s = s, old_s - quotient * s26    return old_s % q2728def vandermonde(q):29    return [[pow(i, j, q) for j in range(q)] for i in range(q)]3031def vandermonde_det_mod(q):32    out = 133    for i in range(q):34        for j in range(i + 1, q):35            out = out * (j - i) % q36    return out3738def invert_mod(matrix, q):39    n = len(matrix)40    aug = [list(row) + [int(i == j) for j in range(n)] for i, row in enumerate(matrix)]41    for col in range(n):42        pivot = next((r for r in range(col, n) if gcd(aug[r][col] % q, q) == 1), None)43        if pivot is None:44            return None45        aug[col], aug[pivot] = aug[pivot], aug[col]46        inv = mod_inverse(aug[col][col], q)47        aug[col] = [x * inv % q for x in aug[col]]48        for r in range(n):49            if r != col and aug[r][col] % q:50                factor = aug[r][col]51                aug[r] = [(x - factor * y) % q for x, y in zip(aug[r], aug[col])]52    return [row[n:] for row in aug]5354def axis_transform(values, q, d, rows):55    size = q ** d56    coeff = list(values)57    for stride in strides(q, d):58        nxt = [0] * size59        for base in range(size):60            if (base // stride) % q:61                continue62            line = [coeff[base + t * stride] for t in range(q)]63            for e in range(q):64                nxt[base + e * stride] = sum(rows[e][t] * line[t] for t in range(q)) % q65        coeff = nxt66    return coeff6768def gfq_eval(coeff, q, cs):69    out = []70    for x in cs:71        acc = 072        for value, e in zip(coeff, cs):73            if value:74                term = value75                for xi, ei in zip(x, e):76                    term = term * pow(xi, ei, q) % q77                acc = (acc + term) % q78        out.append(acc)79    return out8081def int_mobius(values, q, d):82    size = q ** d83    coeff = list(values)84    for stride in strides(q, d):85        for level in range(q - 1, 0, -1):86            for index in range(size):87                if (index // stride) % q == level:88                    coeff[index] -= coeff[index - stride]89    return coeff9091def int_eval(coeff, cs):92    out = []93    for x in cs:94        acc = 095        for value, m in zip(coeff, cs):96            if value and all(xi >= mi for xi, mi in zip(x, m)):97                acc += value98        out.append(acc)99    return out100101def degree_of(coeff, cs):102    return max((sum(e) for value, e in zip(coeff, cs) if value), default=-1)103104def tensor_power(m, d, q):105    out = np.array(m, dtype=np.int64)106    for _ in range(d - 1):107        out = np.kron(out, np.array(m, dtype=np.int64)) % q108    return out109110def eval_matrix(q, d):111    cs = cells(q, d)112    return np.array([[np.prod([pow(xi, ei, q) for xi, ei in zip(x, e)]) % q for e in cs] for x in cs], dtype=np.int64)113114def int_basis_matrix(q, d):115    cs = cells(q, d)116    return np.array([[int(all(xi >= mi for xi, mi in zip(x, m))) for m in cs] for x in cs], dtype=np.int64)117118def int_mobius_matrix(q, d):119    n = q ** d120    columns = [int_mobius([int(i == j) for i in range(n)], q, d) for j in range(n)]121    return np.array(columns, dtype=np.int64).T122123def xor_anf_masks(table, d):124    out = [0] * (1 << d)125    for s in range(1 << d):126        acc = 0127        sub = s128        while True:129            acc ^= table[sub]130            if sub == 0:131                break132            sub = (sub - 1) & s133        out[s] = acc134    return out135136def signed_subset_real(table, d):137    out = [0] * (1 << d)138    for s in range(1 << d):139        acc = 0140        sub = s141        while True:142            if table[sub]:143                acc += -1 if bin(s ^ sub).count("1") & 1 else 1144            if sub == 0:145                break146            sub = (sub - 1) & s147        out[s] = acc148    return out149150def reversed_mask(s, d):151    return sum(((s >> i) & 1) << (d - 1 - i) for i in range(d))152153def code_to_table(code, n):154    return [(code >> i) & 1 for i in range(n)]155156def named(q, d):157    cs = cells(q, d)158    rules = (159        ("void", lambda c: all(v == c[0] for v in c)),160        ("tree", lambda c: all(a == 0 or c[a] == 0 for a in range(d))),161        ("carpet", lambda c: sum(c) <= 1),162        ("net", lambda c: sum(c) >= d - 1),163    )164    return [(name, [int(rule(c)) for c in cs]) for name, rule in rules]165166def fractal_shapes():167    carpet = [int(c != (1, 1)) for c in cells(3, 2)]168    sponge = [int(c.count(1) <= 1) for c in cells(3, 3)]169    return (("sierpinski carpet", 2, carpet), ("menger sponge", 3, sponge))170171def histogram(counts):172    return ", ".join(f"{k}: {counts[k]}" for k in sorted(counts))173174def verdict(ok):175    return "PASS" if ok else "FAIL"176177def block_inverse():178    print("THE INVERSE VANDERMONDE OVER GF(q)")179    for q in PRIME_BASES:180        v = vandermonde(q)181        inv = invert_mod(v, q)182        ok = inv is not None and np.array_equal(np.array(inv) @ np.array(v) % q, np.eye(q, dtype=np.int64))183        print(f"  q={q}  gcd(det, q) = {gcd(vandermonde_det_mod(q), q)}  invertible {inv is not None}  {verdict(ok)}")184    inv2 = invert_mod(vandermonde(2), 2)185    print(f"  q=2  inverse {inv2}  {verdict(inv2 == [[1, 0], [1, 1]])}")186    for q in COMPOSITE_BASES:187        inv = invert_mod(vandermonde(q), q)188        print(f"  q={q}  gcd(det, q) = {gcd(vandermonde_det_mod(q), q)}  invertible {inv is not None}  {verdict(inv is None)}")189    print()190191def block_identity():192    print("THE TENSORED TRANSFORM INVERTS EVALUATION, ALL VALUE TABLES AT ONCE")193    for q, d in IDENTITY_CASES:194        n = q ** d195        t = tensor_power(invert_mod(vandermonde(q), q), d, q)196        e = eval_matrix(q, d)197        eye = np.eye(n, dtype=np.int64)198        ok = np.array_equal(t @ e % q, eye) and np.array_equal(e @ t % q, eye)199        print(f"  q={q} D={d}  T E = E T = I on {n} coordinates  covers all {q}^{n} tables  {verdict(ok)}")200    print("THE INTEGER MOBIUS INVERTS THE DOWNWARD-CLOSED BASIS")201    for q, d in ((2, 2), (2, 3), (2, 4), (3, 2), (3, 3), (4, 2), (4, 3), (6, 2), (9, 2)):202        n = q ** d203        ok = np.array_equal(int_basis_matrix(q, d) @ int_mobius_matrix(q, d), np.eye(n, dtype=np.int64))204        print(f"  q={q} D={d}  B M = I over Z on {n} coordinates  {verdict(ok)}")205    print()206207def block_base2():208    print("BASE 2, EVERY DESIGN, AGAINST THE CLASSICAL XOR ANF")209    rows = invert_mod(vandermonde(2), 2)210    for d in BASE2_DIMS:211        cs = cells(2, d)212        n = 1 << d213        bad_rt = bad_anf = bad_naive = bad_int = 0214        for code in range(1 << n):215            table = code_to_table(code, n)216            coeff = axis_transform(table, 2, d, rows)217            if gfq_eval(coeff, 2, cs) != table:218                bad_rt += 1219            masked = xor_anf_masks([table[reversed_mask(s, d)] for s in range(n)], d)220            support = [int(v != 0) for v in coeff]221            if [masked[reversed_mask(s, d)] for s in range(n)] != support:222                bad_anf += 1223            if masked != support:224                bad_naive += 1225            integer = int_mobius(table, 2, d)226            real = signed_subset_real([table[reversed_mask(s, d)] for s in range(n)], d)227            if int_eval(integer, cs) != table or [real[reversed_mask(s, d)] for s in range(n)] != integer:228                bad_int += 1229        ok = bad_rt == 0 and bad_anf == 0 and bad_int == 0230        print(f"  D={d}  designs {1 << n:>5}  GF(2) roundtrip fails {bad_rt}  XOR ANF diffs after reversal {bad_anf}  diffs without reversal {bad_naive}  integer fails {bad_int}  {verdict(ok)}")231    print()232233def block_base3():234    print("BASE 3, D=2, ALL 512 DESIGNS")235    cs = cells(3, 2)236    rows = invert_mod(vandermonde(3), 3)237    bad_gf = bad_int = 0238    hist_gf, hist_int = {}, {}239    for code in range(512):240        table = code_to_table(code, 9)241        coeff = axis_transform(table, 3, 2, rows)242        if gfq_eval(coeff, 3, cs) != table:243            bad_gf += 1244        integer = int_mobius(table, 3, 2)245        if int_eval(integer, cs) != table:246            bad_int += 1247        dg, di = degree_of(coeff, cs), degree_of(integer, cs)248        hist_gf[dg] = hist_gf.get(dg, 0) + 1249        hist_int[di] = hist_int.get(di, 0) + 1250    print(f"  designs 512  GF(3) roundtrip fails {bad_gf}  integer roundtrip fails {bad_int}  {verdict(bad_gf == 0 and bad_int == 0)}")251    print(f"  GF(3) degree histogram    {histogram(hist_gf)}")252    print(f"  integer degree histogram  {histogram(hist_int)}")253    print(f"  designs of GF(3) degree 1: {hist_gf.get(1, 0)}  {verdict(hist_gf.get(1, 0) == 0)}")254    print()255256def block_table():257    print("THE BASE-3 DEGREE TABLE OF THE FOUR RULES")258    rows = invert_mod(vandermonde(3), 3)259    ok = True260    for d in (2, 3):261        cs = cells(3, d)262        for (name, table), claim in zip(named(3, d), CLAIMED_TABLE[d]):263            dg = degree_of(axis_transform(table, 3, d, rows), cs)264            di = degree_of(int_mobius(table, 3, d), cs)265            ok = ok and dg == claim266            print(f"  D={d} {name:<7} cells {sum(table):>2}/{3 ** d}  GF(3) deg {dg}  page {claim}  integer deg {di}  {verdict(dg == claim)}")267    print(f"  all eight entries  {verdict(ok)}")268    print()269    print("THE FRACTALS THE ROWS ARE NAMED AFTER")270    for name, d, shape in fractal_shapes():271        cs = cells(3, d)272        dg = degree_of(axis_transform(shape, 3, d, rows), cs)273        di = degree_of(int_mobius(shape, 3, d), cs)274        rule = dict(named(3, d))["carpet"]275        rdg = degree_of(axis_transform(rule, 3, d, rows), cs)276        print(f"  {name:<18} D={d}  cells {sum(shape):>2}/{3 ** d}  GF(3) deg {dg}  integer deg {di}  ceiling D(q-1) = {2 * d}")277        print(f"  carpet row         D={d}  cells {sum(rule):>2}/{3 ** d}  GF(3) deg {rdg}  same set as the fractal {rule == shape}")278    print()279280def main():281    start = time.time()282    block_inverse()283    block_identity()284    block_base2()285    block_base3()286    block_table()287    print(f"domain: base 2 at D = {', '.join(map(str, BASE2_DIMS))} exhaustive, base 3 at D = 2 exhaustive, base 3 at D = 3 by the matrix identity")288    print(f"wall {time.time() - start:.1f}s")289290if __name__ == "__main__":291    main()