verify.py

1import cmath2import math3import sys4import time56# DESIGNS78def strings(q, F, L):9    out = [0]10    for _ in range(L):11        out = [v * q + d for v in out for d in F]12    return sorted(set(v for v in out if v > 0))1314def indigits(n, q, F):15    if n == 0:16        return False17    while n:18        if n % q not in F:19            return False20        n //= q21    return True2223def member(q, F, N):24    return [False] + [indigits(n, q, F) for n in range(1, N + 1)]2526def report(name, got, want, tol):27    ok = abs(got - want) <= tol28    print("  %-56s %.3e <= %.1e  %s" % (name, abs(got - want), tol, "ok" if ok else "FAILED"))29    assert ok, "%s: got %r, wanted %r within %r" % (name, got, want, tol)3031def exact(name, got, want):32    print("  %-56s %s  %s" % (name, want, "ok" if got == want else "FAILED"))33    assert got == want, "%s: got %r, wanted %r" % (name, got, want)3435# INVERSE3637def inverse(q, F, N):38    inside = member(q, F, N)39    assert inside[1], "1 must lie in the design"40    nu = [0] * (N + 1)41    nu[1] = 142    for n in range(2, N + 1):43        total = 044        d = 245        while d * d <= n:46            if n % d == 0:47                if inside[d]:48                    total += nu[n // d]49                e = n // d50                if e != d and inside[e]:51                    total += nu[n // e]52            d += 153        if inside[n]:54            total += nu[1]55        nu[n] = -total56    return nu, inside5758def convolve(nu, inside, N):59    out = [0] * (N + 1)60    for a in range(1, N + 1):61        if not inside[a]:62            continue63        for b in range(1, N // a + 1):64            out[a * b] += nu[b]65    return out6667def mobius(N):68    mu = [1] * (N + 1)69    prime = [True] * (N + 1)70    for p in range(2, N + 1):71        if not prime[p]:72            continue73        for m in range(p, N + 1, p):74            if m != p:75                prime[m] = False76            mu[m] = -mu[m]77        for m in range(p * p, N + 1, p * p):78            mu[m] = 079    mu[0] = 080    return mu8182def semigroup(inside, N):83    reach = [False] * (N + 1)84    reach[1] = True85    for n in range(2, N + 1):86        d = 187        while d * d <= n:88            if n % d == 0:89                if inside[d] and reach[n // d]:90                    reach[n] = True91                    break92                if inside[n // d] and reach[d]:93                    reach[n] = True94                    break95            d += 196    return reach9798def check_inverse():99    print("THE REPLACEMENT IDENTITY zeta_F N_F = 1")100    N = 3 ** 8101    nu, inside = inverse(3, [0, 1], N)102    conv = convolve(nu, inside, N)103    exact("base 3 {0,1}: (1_S * nu_F)(1) to n = 6561", conv[1], 1)104    exact("base 3 {0,1}: (1_S * nu_F)(n) = 0 for 2 <= n <= 6561", sum(abs(v) for v in conv[2:]), 0)105    reach = semigroup(inside, N)106    for n in (9, 27, 36):107        exact("base 3 {0,1}: %d in the semigroup, nu_F = 0" % n, (reach[n], nu[n]), (True, 0))108    for n in (16, 48, 52):109        exact("base 3 {0,1}: %d in the semigroup, outside S_F" % n, (reach[n], inside[n]), (True, False))110    for n in range(1, N + 1):111        assert nu[n] == 0 or reach[n], "nu_F supported outside the semigroup at n = %d" % n112    exact("base 3 {0,1}: support inside the semigroup, strictly", True, True)113    for q, L in ((2, 13), (3, 8)):114        M = q ** L115        nuf, _ = inverse(q, list(range(q)), M)116        mu = mobius(M)117        exact("base %d full digit set: nu_F = mu term for term to n = %d" % (q, M), nuf[1:] == mu[1:], True)118    return nu, inside119120# THE WALL121122def repunit(q, a):123    return (q ** a - 1) // (q - 1)124125def witness(q, F):126    missing = sorted(set(range(q)) - set(F))127    inner = [c for c in missing if c >= 2]128    if inner:129        c = min(inner)130        return "repunit", repunit(q, c), repunit(q, c + 1)131    if q % 2:132        return "odd base", 2, (q * q + 1) // 2133    return "even base", q * q - 1, q * q + 1134135def check_wall():136    print("THE WALL: the indicator is multiplicative only at the full digit set")137    tally = {"repunit": 0, "odd base": 0, "even base": 0}138    seen = blind = full = 0139    for q in range(2, 8):140        for code in range(1, 1 << q):141            F = [d for d in range(q) if code >> d & 1]142            seen += 1143            if 1 not in F:144                blind += 1145                continue146            if len(F) == q:147                full += 1148                continue149            kind, m, n = witness(q, F)150            assert math.gcd(m, n) == 1, "witness not coprime at q = %d, F = %s" % (q, F)151            fm = indigits(m, q, F)152            fn = indigits(n, q, F)153            fmn = indigits(m * n, q, F)154            assert fmn != (fm and fn), "witness failed at q = %d, F = %s, pair (%d, %d)" % (q, F, m, n)155            tally[kind] += 1156    exact("nonempty digit sets examined, 2 <= q <= 7", seen, 246)157    exact("sets with 1 outside F, killed at f(1)", blind, 120)158    exact("full digit sets, multiplicative", full, 6)159    exact("sets witnessed, by branch", tally, {"repunit": 114, "odd base": 3, "even base": 3})160    exact("base 3, F = {0,1}: the constructed pair", witness(3, [0, 1]), ("repunit", 4, 13))161    exact("base 2, F = {1}: the constructed pair", witness(2, [1]), ("even base", 3, 5))162    exact("base 12, F = {1}: the constructed pair", witness(12, [1]), ("repunit", 13, 157))163164# THE POSITION PRODUCT165166def position(q, F, L):167    Q = q ** L168    root = [cmath.exp(-2j * math.pi * j / Q) for j in range(Q)]169    grid = []170    for a in range(Q):171        z = complex(1)172        for i in range(L):173            w = 0j174            for d in F:175                w += root[(-d * (q ** i) * a) % Q]176            z *= w177        grid.append(z)178    live = set(strings(q, F, L))179    if 0 in F:180        live.add(0)181    worst = 0.0182    for n in range(Q):183        acc = 0j184        step = (n % Q)185        idx = 0186        for a in range(Q):187            acc += grid[a] * root[idx]188            idx += step189            if idx >= Q:190                idx -= Q191        acc /= Q192        want = 1.0 if n in live else 0.0193        worst = max(worst, abs(acc - want))194    return worst195196def check_position():197    print("THE POSITION PRODUCT: int G_L(t) e(-nt) dt is the indicator of the level")198    report("base 3 {0,1}, L = 5, every n < 243", position(3, [0, 1], 5), 0.0, 1e-10)199    report("base 10 missing 9, L = 3, every n < 1000", position(10, list(range(9)), 3), 0.0, 1e-9)200201# THE DESIGN ZETA202203def design_zeta(q, F, s, J=260, cut=14.0):204    k = len(F)205    nonzero = [a for a in F if a > 0]206    gamma = [float(k)] + [float(sum(a ** l for a in F)) for l in range(1, J + 1)]207    depth = 1208    while q ** depth < 5000:209        depth += 1210    small = strings(q, F, depth)211    value = {}212    for j in range(J, -1, -1):213        w = s + j214        if w.real >= cut:215            value[j] = sum(complex(n) ** (-w) for n in small)216            continue217        acc = sum(complex(a) ** (-w) for a in nonzero)218        binom = complex(1)219        for l in range(1, J - j + 1):220            binom = binom * (-w - (l - 1)) / l221            acc += binom * q ** (-w - l) * gamma[l] * value[j + l]222        value[j] = acc / (1.0 - k * q ** (-w))223    return value[0]224225def check_zeta():226    print("THE DESIGN ZETA: the peel evaluator, then the two boxed zeros")227    for q, F, L, s in ((3, [0, 1], 15, 2.0), (10, list(range(9)), 5, 2.0)):228        k = len(F)229        brute = sum(float(n) ** (-s) for n in strings(q, F, L))230        tail = k * (k * q ** (-s)) ** L / (1.0 - k * q ** (-s))231        got = design_zeta(q, F, complex(s, 0.0))232        report("base %d, k = %d: peel against the direct sum at s = 2" % (q, k), got.real, brute, tail)233    rho3 = complex(0.720790, 28.605680)234    rho10 = complex(1.001590, 2.739200)235    at3 = abs(design_zeta(3, [0, 1], rho3))236    off3 = abs(design_zeta(3, [0, 1], rho3 + 1j))237    at10 = abs(design_zeta(10, list(range(9)), rho10))238    off10 = abs(design_zeta(10, list(range(9)), rho10 + 1j))239    print("  base 3 {0,1}:      |zeta_F| = %.3e at the box centre, %.3e one unit up" % (at3, off3))240    print("  base 10 missing 9: |zeta_F| = %.3e at the box centre, %.3e one unit up" % (at10, off10))241    assert at3 < 1e-3 < off3, "base 3 box centre is not small against its control"242    assert at10 < 1e-3 < off10, "base 10 box centre is not small against its control"243    exact("both centres below 1e-3, both controls above", True, True)244    a = 2245    base, scaled = [0, 1], [0, 2]246    left = set(strings(3, scaled, 9))247    right = set(a * n for n in strings(3, base, 9))248    exact("base 3: S_{2F} = 2 S_F below 3^9", left == right, True)249    probe = complex(0.6, 3.1)250    report("base 3: zeta_{0,2}(s) = 2^(-s) zeta_{0,1}(s)", design_zeta(3, scaled, probe), complex(a) ** (-probe) * design_zeta(3, base, probe), 1e-11)251252# THE LYNDON PRODUCT253254def lyndon(k, L):255    total = 0256    for d in range(1, L + 1):257        if L % d == 0:258            total += moebius_int(d) * k ** (L // d)259    assert total % L == 0, "Lyndon count is not an integer at k = %d, L = %d" % (k, L)260    return total // L261262def moebius_int(n):263    result, d = 1, 2264    while d * d <= n:265        if n % d == 0:266            n //= d267            if n % d == 0:268                return 0269            result = -result270        d += 1271    return -result if n > 1 else result272273def check_lyndon():274    print("THE LYNDON PRODUCT: the free monoid on F factors over its Lyndon words")275    top = 16276    for k in (2, 3, 4, 9, 10):277        series = [1] + [0] * top278        for L in range(1, top + 1):279            c = lyndon(k, L)280            factor = [0] * (top + 1)281            for j in range(0, top // L + 1):282                factor[j * L] = math.comb(c + j - 1, j)283            fresh = [0] * (top + 1)284            for i in range(top + 1):285                if not series[i]:286                    continue287                for j in range(0, top + 1 - i, L):288                    fresh[i + j] += series[i] * factor[j]289            series = fresh290        want = [k ** i for i in range(top + 1)]291        exact("k = %2d: the Lyndon product is 1/(1 - k u) through u^%d" % (k, top), series, want)292    exact("c_2(L), L = 1..10, is A001037", [lyndon(2, L) for L in range(1, 11)], [2, 1, 2, 3, 6, 9, 18, 30, 56, 99])293294# VERIFY295296def main():297    clock = time.time()298    check_wall()299    check_inverse()300    check_position()301    check_lyndon()302    check_zeta()303    print("all checks green in %.1f seconds" % (time.time() - clock))304305if __name__ == "__main__":306    main()