census.py

1import sys2import time3from fractions import Fraction4from itertools import product56import numpy as np7from mpmath import mp, zeta89Q = 310D = 311PARITY = {(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)}12DIGITS = [v for v in product(range(Q), repeat=D) if tuple(c % 2 for c in v) in PARITY]13K = len(DIGITS)14WALK = 1215BASE_LOCAL = Fraction(sum(1 for v in DIGITS if any(v)), K)16BRACKET = BASE_LOCAL * Fraction(Q ** D, Q ** D - 1)17LEGAL = np.zeros((Q,) * D, dtype=bool)18for _v in DIGITS:19    LEGAL[_v] = True202122def mobius(limit):23    mu = np.ones(limit + 1, dtype=np.int64)24    mu[0] = 025    composite = np.zeros(limit + 1, dtype=bool)26    primes = []27    for i in range(2, limit + 1):28        if not composite[i]:29            primes.append(i)30            mu[i] = -131        for p in primes:32            if i * p > limit:33                break34            composite[i * p] = True35            if i % p == 0:36                mu[i * p] = 037                break38            mu[i * p] = mu[i] * mu[p]39    return mu404142def divisible_count(d, n):43    counts = np.zeros((d, d, d), dtype=np.int64)44    counts[0, 0, 0] = 145    place = 1 % d46    for _ in range(n):47        nxt = np.zeros_like(counts)48        for v in DIGITS:49            shift = (place * v[0] % d, place * v[1] % d, place * v[2] % d)50            nxt += np.roll(counts, shift, axis=(0, 1, 2))51        counts = nxt52        place = place * Q % d53    return int(counts[0, 0, 0])545556def in_design(pts, n):57    ok = np.ones(len(pts), dtype=bool)58    rest = pts59    for _ in range(n):60        dig = rest % Q61        ok &= LEGAL[dig[:, 0], dig[:, 1], dig[:, 2]]62        rest = rest // Q63    return int(np.count_nonzero(ok))646566def primitive_box(side):67    axis = np.arange(side + 1, dtype=np.int64)68    grid = np.stack(np.meshgrid(axis, axis, axis, indexing="ij"), axis=-1)69    pts = grid.reshape(-1, D)70    return pts[np.gcd.reduce(pts, axis=1) == 1]717273def visible_count(n, cutoff):74    span = Q ** n75    mu = mobius(cutoff)76    head = 077    for d in range(1, cutoff + 1):78        if mu[d]:79            head += int(mu[d]) * (divisible_count(d, n) - 1)80    partial = np.zeros(span, dtype=np.int64)81    for d in range(1, cutoff + 1):82        if mu[d]:83            partial[d::d] += mu[d]84    tail = 085    for g in range(cutoff + 1, span):86        side = (span - 1) // g87        weight = int(partial[g])88        if weight == 0:89            continue90        tail += weight * in_design(g * primitive_box(side), n)91    return head - tail929394def brute_count(n):95    digits = np.array(DIGITS, dtype=np.int64)96    total = 097    for lead in digits:98        pts = lead.reshape(1, D)99        for _ in range(n - 1):100            pts = (pts[:, None, :] * Q + digits[None, :, :]).reshape(-1, D)101        total += int(np.count_nonzero(np.gcd.reduce(pts, axis=1) == 1))102    return total103104105def ladder(n):106    return max(8, min(int(round(Q ** (n / 2))), 120))107108109def main():110    top = int(sys.argv[1]) if len(sys.argv) > 1 else 9111    mp.dps = 30112    delta = mp.mpf(BRACKET.numerator) / BRACKET.denominator / zeta(D)113    print(f"design q={Q} D={D} k={K} bracket={BRACKET} delta={mp.nstr(delta, 10)}")114    print(f"domain n=1..{top} coordinates 0..{Q ** top - 1}")115    for n in range(1, top + 1):116        cut = ladder(n)117        clock = time.time()118        a = visible_count(n, cut)119        line = f"A({n}) = {a}  G={cut}  {time.time() - clock:.1f}s"120        if n <= 6:121            line += f"  brute={brute_count(n)}"122        if n >= 7:123            gap = (delta * mp.mpf(K) ** n - a) / mp.mpf(WALK) ** n124            line += f"  gap/{WALK}^n = {mp.nstr(gap, 5)}"125        print(line, flush=True)126    if top >= 9:127        low = visible_count(9, 100)128        high = visible_count(9, 150)129        print(f"cutoff check A(9) G=100 {low} G=150 {high} agree={low == high}")130131132main()