verify.py

1# MENGER23import cmath4from decimal import Decimal, getcontext5from fractions import Fraction6from itertools import product7from math import comb, factorial, gcd, log89getcontext().prec = 1001011DESIGN = [v for v in product(range(3), repeat=3) if sum(d == 1 for d in v) <= 1]12CENSUS = [0, 60, 1434, 32268, 721524, 15141288]13MENGER = 0.25162086845125514FULL = 0.28674742843447915FULL_DECIMALS = "0.286747428434478734107892712789"16MENGER_DECIMALS = "0.251620868451255089179675855473"1718def fibre(coords):19    best = 020    for w in product(range(3), repeat=len(coords)):21        best = max(best, sum(1 for v in DESIGN if tuple(v[i] for i in coords) == w))22    return best2324def points(level):25    pts = [(0, 0, 0)]26    for _ in range(level):27        pts = [(3 * x + a, 3 * y + b, 3 * z + c) for (x, y, z) in pts for (a, b, c) in DESIGN]28    return pts2930def split(level):31    low = level // 232    return points(low), points(level - low), 3 ** low3334def census(level):35    lows, highs, shift = split(level)36    total = 037    for hx, hy, hz in highs:38        bx, by, bz = hx * shift, hy * shift, hz * shift39        for lx, ly, lz in lows:40            x = bx + lx41            y = by + ly42            if gcd(x, y) != 1:43                continue44            z = bz + lz45            if gcd(x, z) == 1 and gcd(y, z) == 1:46                total += 147    return total4849def base_census(level):50    lows, highs, shift = split(level)51    total = 052    for hx, hy, hz in highs:53        bx, by, bz = hx * shift, hy * shift, hz * shift54        for lx, ly, lz in lows:55            zeros = ((bx + lx) % 3 == 0) + ((by + ly) % 3 == 0) + ((bz + lz) % 3 == 0)56            if zeros <= 1:57                total += 158    return total5960def parity_census(level):61    lows, highs, shift = split(level)62    total = 063    for hx, hy, hz in highs:64        bx, by, bz = hx * shift, hy * shift, hz * shift65        for lx, ly, lz in lows:66            odd = ((bx + lx) & 1) + ((by + ly) & 1) + ((bz + lz) & 1)67            if odd >= 2:68                total += 169    return total7071def parity_transfer(level):72    state = {(0, 0, 0): Fraction(1)}73    for _ in range(level):74        nxt = {}75        for s, w in state.items():76            for v in DESIGN:77                t = tuple((s[i] + (v[i] == 1)) % 2 for i in range(3))78                nxt[t] = nxt.get(t, Fraction(0)) + w * Fraction(1, 20)79        state = nxt80    return sum(w for s, w in state.items() if sum(s) >= 2)8182def parity_closed(level):83    return Fraction(1, 2) - Fraction(3, 4) * Fraction(3, 5) ** level + Fraction(1, 4) * Fraction(-1, 5) ** level8485def primes_to(limit):86    sieve = bytearray([1]) * (limit + 1)87    sieve[0:2] = b"\x00\x00"88    for i in range(2, int(limit ** 0.5) + 1):89        if sieve[i]:90            sieve[i * i::i] = bytearray(len(range(i * i, limit + 1, i)))91    return [i for i in range(2, limit + 1) if sieve[i]]9293def euler(primes, skip):94    out = 1.095    for p in primes:96        if p != skip:97            out *= 1.0 - 3.0 / p ** 2 + 2.0 / p ** 398    return out99100def bernoulli(kmax):101    out = [Fraction(0)] * (kmax + 1)102    out[0] = Fraction(1)103    for m in range(1, kmax + 1):104        out[m] = -sum(comb(m + 1, j) * out[j] for j in range(m)) / (m + 1)105    return out106107BERNOULLI = bernoulli(32)108109def decimalise(fr):110    return Decimal(fr.numerator) / Decimal(fr.denominator)111112def zeta(n, cut=100, terms=15):113    total = sum(Decimal(1) / Decimal(j) ** n for j in range(1, cut))114    total += Decimal(cut) ** (1 - n) / Decimal(n - 1)115    total += Decimal(1) / (Decimal(2) * Decimal(cut) ** n)116    for k in range(1, terms + 1):117        rising = 1118        for i in range(2 * k - 1):119            rising *= n + i120        total += decimalise(BERNOULLI[2 * k] * Fraction(rising, factorial(2 * k))) * Decimal(cut) ** (1 - n - 2 * k)121    return total122123ZETA = {}124125def zeta_at(n):126    if n not in ZETA:127        ZETA[n] = zeta(n)128    return ZETA[n]129130def mobius(limit):131    mu = [1] * (limit + 1)132    for p in primes_to(limit):133        for m in range(p, limit + 1, p):134            mu[m] = -mu[m]135        for m in range(p * p, limit + 1, p * p):136            mu[m] = 0137    return mu138139def prime_zeta(n, top, mu):140    total = Decimal(0)141    k = 1142    while k * n <= top:143        if mu[k]:144            total += Decimal(mu[k]) / Decimal(k) * zeta_at(k * n).ln()145        k += 1146    return total147148def constant(bound, nmax, top):149    small = primes_to(bound)150    mu = mobius(top)151    total = Decimal(0)152    for p in small:153        x = Decimal(1) / Decimal(p)154        total += ((1 - x) ** 2 * (1 + 2 * x)).ln()155    for n in range(2, nmax + 1):156        weight = Fraction(-2 + (-1) ** (n + 1) * 2 ** n, n)157        tail = prime_zeta(n, top, mu) - sum(Decimal(1) / Decimal(p) ** n for p in small)158        total += decimalise(weight) * tail159    return total.exp()160161def multiplier(p, t):162    total = sum(cmath.exp(2j * cmath.pi * sum(t[i] * v[i] for i in range(3)) / p) for v in DESIGN)163    return abs(total) / 20.0164165def order_of_three(p):166    k, r = 1, 3 % p167    while r != 1:168        r = (r * 3) % p169        k += 1170    return k171172def contraction(p):173    d = order_of_three(p)174    best = 0.0175    for t in product(range(p), repeat=3):176        if not any(t):177            continue178        run = 1.0179        u = t180        for _ in range(d):181            run *= multiplier(p, u)182            u = tuple((3 * x) % p for x in u)183        best = max(best, run ** (1.0 / d))184    return best185186def main():187    got = len(DESIGN)188    assert got == 20, f"design size: got {got}, want 20"189    got = sum(1 for v in DESIGN if sum(d == 0 for d in v) <= 1)190    assert got == 13, f"vectors with at most one zero: got {got}, want 13"191    print(f"design: {len(DESIGN)} vectors, {got} of them with at most one zero")192    for pair in [(0, 1), (0, 2), (1, 2)]:193        got = fibre(pair)194        assert got == 3, f"pair fibre {pair}: got {got}, want 3"195    for i in range(3):196        got = fibre((i,))197        assert got == 8, f"single fibre {i}: got {got}, want 8"198    alpha = log(20 / 3) / log(3)199    single = log(20 / 8) / log(3)200    assert alpha > 1.0, f"pair exponent: got {alpha}, want above 1"201    assert single < 1.0, f"single exponent: got {single}, want below 1"202    print(f"fibres: pair 3 giving exponent {alpha:.6f}, single 8 giving exponent {single:.6f}")203204    for level in range(1, 7):205        got = census(level)206        want = CENSUS[level - 1]207        assert got == want, f"census level {level}: got {got}, want {want}"208        print(f"level {level}: {got} pairwise-coprime points of {20 ** level}, density {got / 20 ** level:.9f}")209210    for level in range(1, 6):211        got = Fraction(base_census(level), 20 ** level)212        want = Fraction(13, 20)213        assert got == want, f"base factor at level {level}: got {got}, want {want}"214    print("base prime: no coordinate pair divisible by 3 for exactly 13/20 of every level up to 5")215216    for level in range(5):217        got = Fraction(parity_census(level), 20 ** level)218        want = parity_closed(level)219        assert got == want, f"parity census level {level}: got {got}, want {want}"220    for level in range(41):221        got = parity_transfer(level)222        want = parity_closed(level)223        assert got == want, f"parity transfer level {level}: got {got}, want {want}"224    print(f"bias at two: closed form exact to level 40, level 5 value {float(parity_closed(5)):.6f}")225226    high = constant(100, 60, 340)227    got = str(high)[:32]228    assert got == FULL_DECIMALS, f"full constant to thirty decimals: got {got}, want {FULL_DECIMALS}"229    sponge = high * Decimal(351) / Decimal(400)230    got = str(sponge)[:32]231    assert got == MENGER_DECIMALS, f"sponge constant to thirty decimals: got {got}, want {MENGER_DECIMALS}"232    again = constant(200, 45, 300)233    gap = abs(high - again)234    assert gap < Decimal(10) ** -60, f"two truncation settings: got gap {gap}, want below 1e-60"235    print(f"high precision: full {str(high)[:32]}, sponge {str(sponge)[:32]}, settings agree to {gap:.3e}")236237    primes = primes_to(10 ** 6)238    got = euler(primes, None)239    assert abs(got - FULL) < 1e-6, f"truncated full product: got {got}, want {FULL} within 1e-6"240    got = 0.65 * euler(primes, 3)241    assert abs(got - MENGER) < 1e-6, f"truncated sponge product: got {got}, want {MENGER} within 1e-6"242    print(f"products to 1e6: full {euler(primes, None):.9f}, sponge {got:.9f}")243244    local = Fraction(1) - 3 * Fraction(1, 9) + 2 * Fraction(1, 27)245    assert local == Fraction(20, 27), f"local factor at 3: got {local}, want 20/27"246    ratio = Fraction(13, 20) / local247    assert ratio == Fraction(351, 400), f"ratio: got {ratio}, want 351/400"248    loss = Fraction(1) - ratio249    assert loss == Fraction(49, 400), f"loss: got {loss}, want 49/400"250    print(f"exact rationals: 20/27, {ratio}, {loss} = {float(loss) * 100:.2f} percent")251252    for p in [2, 5, 7, 11, 13, 17, 19]:253        got = contraction(p)254        assert got <= 0.6 + 1e-9, f"contraction at {p}: got {got}, want at most 3/5"255        if p == 2:256            assert abs(got - 0.6) < 1e-9, f"contraction at 2: got {got}, want 3/5"257        print(f"contraction at {p}: {got:.6f}")258259    print("all green")260261if __name__ == "__main__":262    main()