classes.py

1from fractions import Fraction2from itertools import product3from math import pi, prod45import numpy as np6from mpmath import zeta78CORNERS = [tuple(v) for v in product((0, 1), repeat=3)]9Z3 = float(zeta(3))10Z2 = float(zeta(2))1112def code_set(mask):13    return frozenset(v for i, v in enumerate(CORNERS) if mask >> i & 1)1415def code_mask(P):16    return sum(1 << i for i, v in enumerate(CORNERS) if v in P)1718def origin_filled(P):19    return 1 if (0, 0, 0) in P else 02021def prime_factors(n):22    ps, d = [], 223    while d * d <= n:24        if n % d == 0:25            ps.append(d)26            while n % d == 0:27                n //= d28        d += 129    if n > 1:30        ps.append(n)31    return ps3233def radical(n):34    r = 135    for p in prime_factors(n):36        r *= p37    return r3839def mobius(n):40    ps = prime_factors(n)41    r = 142    for p in ps:43        r *= p44    return 0 if r != n else (-1) ** len(ps)4546def divisors(n):47    return [d for d in range(1, n + 1) if n % d == 0]4849def digits_with(q, step, parity):50    return sum(1 for d in range(0, q, step) if d % 2 == parity)5152def corner_count(P, q, step=1):53    return sum(prod(digits_with(q, step, c) for c in v) for v in P)5455def bracket(P, q):56    k = corner_count(P, q)57    return sum(Fraction(mobius(e) * corner_count(P, q, e), k) for e in divisors(radical(q)))5859def rational_part(P, q):60    d = len(next(iter(P)))61    r = bracket(P, q)62    for p in prime_factors(q):63        r *= Fraction(p ** d, p ** d - 1)64    return r6566def bits(P):67    return [v[0] * 4 + v[1] * 2 + v[2] for v in sorted(P)]6869def difference_span(P):70    L = bits(P)71    basis = []72    for w in sorted({a ^ b for a in L for b in L}, reverse=True):73        v = w74        for b in basis:75            v = min(v, v ^ b)76        if v:77            basis.append(v)78    H = {0}79    for b in basis:80        H |= {h ^ b for h in H}81    return len(basis), H8283def classify(P):84    s2, H = difference_span(P)85    if s2 == 3:86        return "span", s287    if all(x in H for x in bits(P)):88        return "inside", s289    return "offset", s29091def weight(P):92    return tuple(sum(1 for v in P if sum(v) == j) for j in range(4))9394def is_subgroup(P):95    L = set(bits(P))96    return 0 in L and all((a ^ b) in L for a in L for b in L)9798def pinned_odd(P):99    return any(all(v[i] == 1 for v in P) for i in range(3))100101def even_band(P):102    return Fraction(8, 7) * (1 - Fraction(origin_filled(P), len(P)))103104def odd_limits(P, q):105    kind, s2 = classify(P)106    naive = rational_part(P, q)107    if kind == "span":108        return [naive]109    corrected = Fraction(8, 7) * naive * (1 - Fraction(1, 2 ** s2))110    if kind == "inside":111        return [corrected]112    return [Fraction(8, 7) * naive, corrected]113114def odd_limit_over_bases(P):115    kind, s2 = classify(P)116    if kind == "span":117        return Fraction(1)118    if kind == "inside":119        return Fraction(8, 7) * (1 - Fraction(1, 2 ** s2))120    return None121122def even_point_fraction(P, q, n):123    e, o = Fraction(q + 1, 2), Fraction(q - 1, 2)124    w = {v: e ** (3 - sum(v)) * o ** sum(v) for v in P}125    k = sum(w.values())126    total = Fraction(0)127    for t in CORNERS:128        lam = sum(w[v] * (-1) ** (t[0] * v[0] + t[1] * v[1] + t[2] * v[2]) for v in P) / k129        total += lam ** n130    return total / 8131132def predicted_level(P, q, n):133    if q % 2 == 0:134        return rational_part(P, q)135    return Fraction(8, 7) * rational_part(P, q) * (1 - even_point_fraction(P, q, n))136137def design(P, q):138    return np.array([v for v in product(range(q), repeat=3) if (v[0] % 2, v[1] % 2, v[2] % 2) in P], dtype=np.int64)139140def triangular(rows):141    M = [list(r) for r in rows]142    piv = 0143    for col in range(3):144        r = next((i for i in range(piv, len(M)) if M[i][col]), None)145        if r is None:146            continue147        M[piv], M[r] = M[r], M[piv]148        for i in range(piv + 1, len(M)):149            while M[i][col]:150                f = M[piv][col] // M[i][col]151                M[piv] = [a - f * b for a, b in zip(M[piv], M[i])]152                M[piv], M[i] = M[i], M[piv]153        piv += 1154    return piv, M155156def lattice(P, q):157    F = design(P, q).tolist()158    rows = [[a - b for a, b in zip(v, F[0])] for v in F[1:]]159    if not rows:160        return 0, 0161    rank, M = triangular(rows)162    if rank < 3:163        return rank, 0164    return 3, abs(M[0][0] * M[1][1] * M[2][2])165166def stack(F, q, n):167    pts = np.zeros((1, 3), dtype=np.int64)168    for _ in range(n):169        pts = (pts[:, None, :] * q + F[None, :, :]).reshape(-1, 3)170    return pts171172def visible(pts):173    g = np.gcd(np.gcd(pts[:, 0], pts[:, 1]), pts[:, 2])174    return int(np.count_nonzero(g == 1))175176def visible_fraction(P, q, n):177    F = design(P, q)178    if n == 1:179        return visible(F) / len(F)180    base = stack(F, q, n - 1)181    hits = 0182    for v in F:183        hits += visible(base * q + v)184    return hits / len(F) ** n185186CARPET = frozenset([(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)])187NET = frozenset([(1, 1, 0), (1, 0, 1), (0, 1, 1), (1, 1, 1)])188TREE = frozenset([(0, 0, 0), (0, 0, 1)])189VOID = frozenset([(0, 0, 0), (1, 1, 1)])190AXES = frozenset([(1, 0, 0), (0, 1, 0), (0, 0, 1)])191TOP = frozenset([(1, 1, 1)])192PINNED = frozenset([(1, 1, 0), (1, 0, 0)])193TWIN_SUB = frozenset([(0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 0)])194TWIN_SPAN = frozenset([(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 1, 1)])195FLAT_CARPET = frozenset([(0, 0), (1, 0), (0, 1)])196197FAMILIES = [("carpet", CARPET), ("net", NET), ("tree", TREE), ("void", VOID)]198199CASES = [200    ("carpet", CARPET, 2, 12),201    ("carpet", CARPET, 4, 5),202    ("net", NET, 4, 5),203    ("tree", TREE, 3, 7),204    ("tree", TREE, 5, 5),205    ("void", VOID, 3, 8),206    ("void", VOID, 5, 5),207    ("twin subgroup", TWIN_SUB, 3, 6),208    ("twin spanning", TWIN_SPAN, 3, 6),209    ("twin subgroup", TWIN_SUB, 4, 4),210    ("twin spanning", TWIN_SPAN, 4, 4),211]212213def band_report(codes):214    bad = [(code_mask(P), q) for P in codes for q in range(2, 41, 2) if rational_part(P, q) != even_band(P)]215    print("even bases: 255 codes x even q<=40, exceptions to (8/7)(1-t/|P|):", len(bad))216    band = sorted({even_band(P) for P in codes})217    print("even band values of delta*zeta(3):", len(band), " ".join(str(x) for x in band))218    flat = {rational_part(FLAT_CARPET, q) for q in range(2, 41, 2)}219    v = float(min(flat)) / Z2220    print("D=2 parity carpet at even q<=40: delta*zeta(2) =", " ".join(str(x) for x in flat),221          "delta =", repr(v), "16/(3 pi^2) =", repr(16 / (3 * pi ** 2)))222223def similarity_report(codes):224    bad = 0225    for P in codes:226        for q in range(3, 76, 2):227            for e in divisors(radical(q)):228                if corner_count(P, q, e) != corner_count(P, q // e):229                    bad += 1230    print("odd bases: 255 codes x odd q<=75, exceptions to k_e(q)=k(q/e):", bad)231232def class_report(codes):233    kinds = {}234    for P in codes:235        kind, s2 = classify(P)236        kinds.setdefault(kind, []).append(s2)237    print("class counts: span", len(kinds["span"]), "inside", len(kinds["inside"]), "offset", len(kinds["offset"]),238          "total", sum(len(v) for v in kinds.values()))239    print("inside by s2:", {s: kinds["inside"].count(s) for s in sorted(set(kinds["inside"]))})240    classes = {}241    for P in codes:242        classes.setdefault(weight(P), set()).add(classify(P)[0])243    mixed = [k for k in classes.values() if len(k) > 1]244    print("weight classes:", len(classes), "mixing regimes:", len(mixed),245          "spanning against not:", sum(1 for k in mixed if "span" in k),246          "corrected against no-limit:", sum(1 for k in mixed if k == {"inside", "offset"}))247    subgroups = sorted(code_mask(P) for P in codes if is_subgroup(P))248    print("subgroup codes:", len(subgroups), " ".join(str(m) for m in subgroups))249    stable = sorted(code_mask(P) for P in codes if odd_limit_over_bases(P) == even_band(P))250    print("parity-stable codes:", len(stable), "equal to subgroup codes:", stable == subgroups)251    print("codes with two subsequential limits:", 255 - len(stable))252253def lattice_report(codes):254    full = [P for P in codes if lattice(P, 2)[0] == 3]255    print("q=2: full-rank codes:", len(full), "smallest |P|:", min(len(P) for P in full),256          "indices:", sorted({lattice(P, 2)[1] for P in full}))257    print("even q=4,6,8,10: lattice indices over all codes:",258          sorted({lattice(P, q)[1] for P in codes for q in (4, 6, 8, 10)}))259    print("q=3: codes of lattice rank <=1:", sum(1 for P in codes if lattice(P, 3)[0] <= 1),260          "codes with a coordinate pinned odd:", sum(1 for P in codes if pinned_odd(P)))261262def family_report():263    print("family table, delta*zeta(3):", " ".join(name for name, _ in FAMILIES))264    print("  every even:", " ".join(str(even_band(P)) for _, P in FAMILIES))265    for q in (3, 5, 7, 9, 11):266        print("  q=%-2d       :" % q, " ".join(str(odd_limits(P, q)[-1]) for _, P in FAMILIES))267    print("  odd limit  :", " ".join(str(odd_limit_over_bases(P)) for _, P in FAMILIES))268    print("carpet rational parts at odd q=3,5,7,9,11:",269          " ".join("%.5f" % float(odd_limits(CARPET, q)[-1]) for q in (3, 5, 7, 9, 11)),270          "k(q):", " ".join(str(corner_count(CARPET, q)) for q in (3, 5, 7, 9, 11)),271          "dips 1/k(7) = 1/%d, k(3)/k(9) = %d/%d"272          % (corner_count(CARPET, 7), corner_count(CARPET, 3), corner_count(CARPET, 9)))273    print("net: odd q<=81 with rational part <= 1:", [q for q in range(3, 82, 2) if rational_part(NET, q) <= 1],274          "B(9) =", bracket(NET, 9), "= 1 - %d/%d" % (corner_count(NET, 3), corner_count(NET, 9)),275          "k(1) =", corner_count(NET, 1))276277def exact_report():278    codes = [code_set(m) for m in range(1, 256)]279    band_report(codes)280    similarity_report(codes)281    class_report(codes)282    lattice_report(codes)283    family_report()284285def measured_report():286    for name, P, q, n in CASES:287        k = len(design(P, q))288        got = visible_fraction(P, q, n)289        level = float(predicted_level(P, q, n)) / Z3290        lims = [float(x) / Z3 for x in odd_limits(P, q)] if q % 2 else [float(rational_part(P, q)) / Z3]291        print("%-14s q=%d n=%2d k=%3d measured=%.6f level=%.6f limits=%s"292              % (name, q, n, k, got, level, " ".join("%.6f" % x for x in lims)))293    print("pinned {110,100} q=3 odd levels:", " ".join("%.6f" % visible_fraction(PINNED, 3, n) for n in (3, 5, 7)),294          "outside the trichotomy, whose value is %.6f" % (float(odd_limits(PINNED, 3)[0]) / Z3))295296def oscillation_report(name, P, q, levels):297    lims = [float(x) / Z3 for x in odd_limits(P, q)]298    print("%s q=%d limits: odd %.6f even %.6f" % (name, q, lims[0], lims[1]))299    for n in range(1, levels + 1):300        print("  n=%d measured=%.6f level=%.6f" % (n, visible_fraction(P, q, n), float(predicted_level(P, q, n)) / Z3))301302def main():303    exact_report()304    print()305    measured_report()306    print()307    oscillation_report("top corner", TOP, 5, 8)308    print()309    oscillation_report("axes", AXES, 3, 7)310311if __name__ == "__main__":312    main()