slices.py

1import numpy as np2from itertools import product34CARPET = frozenset([(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)])5NET = frozenset([(1, 1, 0), (1, 0, 1), (0, 1, 1), (1, 1, 1)])6TREE = frozenset([(0, 0, 0), (0, 0, 1)])78def corners(P, q):9    return np.array([v for v in product(range(q), repeat=3)10                     if (v[0] % 2, v[1] % 2, v[2] % 2) in P], dtype=np.int64)1112def level_points(F, q, n):13    pts = np.zeros((1, 3), dtype=np.int64)14    for _ in range(n):15        pts = (pts[:, None, :] * q + F[None, :, :]).reshape(-1, 3)16    return pts1718def prime_factors(m):19    ps, d = [], 220    while d * d <= m:21        if m % d == 0:22            ps.append(d)23            while m % d == 0:24                m //= d25        d += 126    if m > 1:27        ps.append(m)28    return ps2930def squarefree_upto(limit):31    out = []32    for d in range(1, limit + 1):33        ps, prod = prime_factors(d), 134        for p in ps:35            prod *= p36        if prod == d:37            out.append((d, (-1) ** len(ps)))38    return out3940def squarefree_divisors(m):41    out = [(1, 1)]42    for p in prime_factors(m):43        out = out + [(d * p, -mu) for d, mu in out]44    return out4546def height_counts(F, q, n):47    weights = {}48    for x in F.sum(axis=1).tolist():49        weights[x] = weights.get(x, 0) + 150    top, arr = max(weights), np.array([1], dtype=np.int64)51    for _ in range(n):52        L = len(arr)53        new = np.zeros(q * (L - 1) + top + 1, dtype=np.int64)54        for x, c in weights.items():55            new[x:x + q * L:q] += c * arr56        arr = new57    return arr5859def aggregate_local(F, q, n, p):60    vs, cur = [tuple(v) for v in F.tolist()], {(0, 0, 0): 1}61    for _ in range(n):62        nxt = {}63        for (a, b, c), t in cur.items():64            for v in vs:65                k = ((q * a + v[0]) % p, (q * b + v[1]) % p, (q * c + v[2]) % p)66                nxt[k] = nxt.get(k, 0) + t67        cur = nxt68    tot = sum(t for r, t in cur.items() if sum(r) % p == 0)69    return cur.get((0, 0, 0), 0), tot7071def slice_divisor_counts(F, q, n, s, ds):72    m = n // 273    lo, hi = level_points(F, q, m), level_points(F, q, n - m)74    scale = q ** m75    slo, shi = lo.sum(axis=1).tolist(), hi.sum(axis=1).tolist()76    lol, hil, out = lo.tolist(), hi.tolist(), {}77    for d in ds:78        tab = {}79        for j, v in enumerate(lol):80            k = (slo[j], v[0] % d, v[1] % d, v[2] % d)81            tab[k] = tab.get(k, 0) + 182        tot = 083        for i, v in enumerate(hil):84            need = s - shi[i] * scale85            if need >= 0:86                tot += tab.get((need, (-scale * v[0]) % d, (-scale * v[1]) % d,87                                (-scale * v[2]) % d), 0)88        out[d] = tot89    return out9091def mobius_and_prime_slices(P, q, n):92    F = corners(P, q)93    pts = level_points(F, q, n)94    s = pts.sum(axis=1)95    smax = 3 * (q ** n - 1)96    N = np.bincount(s, minlength=smax + 1)97    A = np.bincount(s[np.gcd.reduce(pts, axis=1) == 1], minlength=smax + 1)98    total = np.zeros(smax + 1, dtype=np.int64)99    for d, mu in squarefree_upto(q ** n - 1):100        total += mu * np.bincount(s[(pts % d == 0).all(axis=1)], minlength=smax + 1)101    hidden = max([int(N[h] - A[h]) for h in range(2, smax + 1)102                  if N[h] and prime_factors(h) == [h]] + [0])103    return int(np.count_nonzero(total[1:] != A[1:])), hidden104105def peel_mismatches(P, q, n):106    F = corners(P, q)107    a, b = level_points(F, q, n), level_points(F, q, n - 1)108    smax = 3 * (q ** n - 1)109    lhs = np.bincount(a.sum(axis=1)[(a % q == 0).all(axis=1)], minlength=smax + 1)110    Nb = np.bincount(b.sum(axis=1), minlength=3 * (q ** (n - 1) - 1) + 1)111    origin, bad = (0, 0, 0) in P, 0112    for h in range(smax + 1):113        r = int(Nb[h // q]) if (origin and h % q == 0 and h // q < len(Nb)) else 0114        bad += int(lhs[h]) != r115    return bad116117def parity_dichotomy(q, n):118    F = corners(TREE, q)119    base = level_points(F, q, n - 1)120    ev = od = seen = evengcd = 0121    for v in F.tolist():122        ch = base * q + np.array(v, dtype=np.int64)123        s, g = ch.sum(axis=1), np.gcd.reduce(ch, axis=1)124        pos = s > 0125        e, o = pos & (s % 2 == 0), pos & (s % 2 == 1)126        ev += int(np.count_nonzero(e))127        od += int(np.count_nonzero(o))128        seen += int(np.count_nonzero(e & (g == 1)))129        evengcd += int(np.count_nonzero(o & (g % 2 == 0)))130    return ev, od, seen, evengcd131132def walk_lambdas(P, q):133    F = corners(P, q).tolist()134    return {t: sum((-1) ** (t[0] * v[0] + t[1] * v[1] + t[2] * v[2]) for v in F) / len(F)135            for t in product((0, 1), repeat=3)}136137def central_line(P, q, n):138    F = corners(P, q)139    sc = 3 * (q ** n - 1) // 2140    ds = squarefree_divisors(sc)141    cnt = slice_divisor_counts(F, q, n, sc, [d for d, _ in ds])142    ps = prime_factors(sc)143    loc = {p: cnt[p] / cnt[1] for p in ps}144    ind = 1.0145    for p in ps:146        ind *= 1 - loc[p]147    return sc, ps, cnt[1], sum(mu * cnt[d] for d, mu in ds), loc, ind148149def order(p, q):150    o, x = 1, q % p151    while x != 1:152        x, o = x * q % p, o + 1153    return o154155print("DOMAIN carpet q=3 to n=6, net q=3 to n=7, tree q=3 to n=7, carpet q=4 to n=3, carpet q=5 to n=4")156print()157print("EXACT LAWS AT EVERY HEIGHT")158for name, P, q, nmax in [("carpet3", CARPET, 3, 4), ("net3", NET, 3, 4),159                         ("carpet4", CARPET, 4, 3), ("carpet5", CARPET, 5, 3)]:160    for n in range(1, nmax + 1):161        bad, hidden = mobius_and_prime_slices(P, q, n)162        peel = peel_mismatches(P, q, n) if q in (3, 5) else "not a prime base"163        print(f"{name} n={n} mobius mismatches={bad} peel mismatches={peel} "164              f"max hidden on a prime slice={hidden}")165pts = level_points(corners(NET, 3), 3, 7)166print(f"net3 n=7 points={len(pts)} with 3 dividing the gcd="167      f"{int(np.count_nonzero(np.gcd.reduce(pts, axis=1) % 3 == 0))}")168print()169170print("PARITY DICHOTOMY, TREE q=3")171for n in (6, 7):172    ev, od, seen, evengcd = parity_dichotomy(3, n)173    print(f"tree3 n={n} even-height points={ev} visible among them={seen} "174          f"odd-height points={od} even gcds among them={evengcd}")175print()176177print("AGGREGATED LOCAL FACTOR OVER THE HEIGHTS DIVISIBLE BY p, CARPET q=3 n=6")178F = corners(CARPET, 3)179for p in (5, 7, 11, 13):180    div, tot = aggregate_local(F, 3, 6, p)181    print(f"p={p}: {div}/{tot} = {div / tot:.6f} vs 1/p^2 = {1 / p ** 2:.6f}")182e3 = level_points(F, 3, 3)183bad = sum(aggregate_local(F, 3, 3, p) != (int(np.count_nonzero((e3 % p == 0).all(axis=1))),184                                          int(np.count_nonzero(e3.sum(axis=1) % p == 0)))185          for p in (2, 5, 7, 11, 13))186print(f"transfer against enumeration at n=3 over p=2,5,7,11,13: mismatches={bad}")187lam = walk_lambdas(CARPET, 3)188num = sum(v ** 6 for v in lam.values()) / 8 * 20 ** 6189den = (1 + lam[(1, 1, 1)] ** 6) / 2 * 20 ** 6190div, tot = aggregate_local(F, 3, 6, 2)191print(f"p=2 counted {div}/{tot} = {div / tot:.7f}")192print(f"p=2 walk formula {round(num)}/{round(den)} = {num / den:.7f}")193print()194195print("CENTRAL SLICE")196for q, nmax in ((3, 6), (5, 4)):197    for n in range(1, nmax + 1):198        sc, ps, N, A, loc, ind = central_line(CARPET, q, n)199        print(f"q={q} n={n} s*={sc}={'*'.join(map(str, ps))} N={N} A={A} delta={A / N:.5f} "200              f"independence-product={ind:.5f} locals "201              + " ".join(f"{p}:{loc[p]:.5f}" for p in ps))202print()203204print("CENTRAL COUNT AND ITS PEEL, CARPET q=3")205hc = [height_counts(corners(CARPET, 3), 3, n) for n in range(15)]206cen = [int(hc[n][3 * (3 ** n - 1) // 2]) for n in range(15)]207pel = [int(hc[n - 1][(3 ** n - 1) // 2]) for n in range(2, 15)]208print("N(s*_n) n=1..8:", ", ".join(str(cen[n]) for n in range(1, 9)))209print("N^(3)(s*_n) n=2..8:", ", ".join(str(pel[n - 2]) for n in range(2, 9)))210for n in range(2, 8):211    direct = slice_divisor_counts(corners(CARPET, 3), 3, n, 3 * (3 ** n - 1) // 2, [3])[3]212    print(f"n={n} peel {direct} vs level {n - 1} at R_n={(3 ** n - 1) // 2}: {pel[n - 2]} "213          f"equal={direct == pel[n - 2]}")214off = sum(cen[n] != 9 * cen[n - 1] - 12 * cen[n - 2] for n in range(2, 15))215off += sum(pel[i] != 9 * pel[i - 1] - 12 * pel[i - 2] for i in range(2, len(pel)))216print(f"(9,-12) recurrence exceptions on both sequences to n=14: {off}")217print(f"peel ratio n=14 = {pel[12] / cen[14]:.9f} vs (sqrt(33)-5)/8 = {(33 ** 0.5 - 5) / 8:.10f}")218print()219220print("THE CENTRAL BILL")221bad = 0222for q in (3, 5, 7):223    for n in range(1, 9):224        ps = prime_factors(3 * (q ** n - 1) // 2)225        bad += 3 not in ps226        bad += (2 in ps) != (q % 4 == 1 or n % 2 == 0)227        for p in ps:228            if p not in (2, 3) and q % p:229                bad += n % order(p, q) != 0230print(f"bill rule exceptions over q=3,5,7 and n=1..8: {bad}")231print(f"R_7(3) = {(3 ** 7 - 1) // 2}, foreign bill at q=3 n=7 = "232      f"{[p for p in prime_factors(3 * (3 ** 7 - 1) // 2) if p != 3]}")233print(f"2^1092 mod 1093^2 = {pow(2, 1092, 1093 ** 2)}")