density.py

1import itertools2from fractions import Fraction3from math import log4import numpy as np56Q = 37NMAX = 148FACTOR_DEG = 69EULER_N = 1010EULER_DEG = 511LEVELS = (10, 12, 14)1213def trim(c):14    while c and c[-1] == 0:15        del c[-1]16    return c1718def divmod_poly(a, b):19    r = list(a)20    db = len(b) - 121    inv = pow(b[-1], Q - 2, Q)22    q = [0] * max(0, len(r) - db)23    while r and len(r) - 1 >= db:24        d = len(r) - 1 - db25        f = (r[-1] * inv) % Q26        q[d] = f27        for i, c in enumerate(b):28            r[i + d] = (r[i + d] - f * c) % Q29        trim(r)30    return tuple(trim(q)), tuple(r)3132def monics(d):33    for tail in itertools.product(range(Q), repeat=d):34        yield tail + (1,)3536def irreducibles(maxdeg):37    out = []38    for d in range(1, maxdeg + 1):39        for p in monics(d):40            if all(divmod_poly(p, s)[1] != () for s in out if len(s) - 1 <= d // 2):41                out.append(p)42    return out4344def label(p):45    parts = []46    for k in range(len(p) - 1, -1, -1):47        c = p[k]48        if c == 0:49            continue50        if k == 0:51            parts.append(str(c))52        else:53            base = "t" if k == 1 else "t^%d" % k54            parts.append(base if c == 1 else "%d%s" % (c, base))55    return " + ".join(parts) if parts else "0"5657def poly_of(i, n):58    return tuple(trim([(i >> k) & 1 for k in range(n)]))5960def masks(primes, n):61    bits = ((np.arange(1 << n)[:, None] >> np.arange(n)) & 1).astype(np.int64)62    out = {}63    for p in primes:64        dp = len(p) - 165        rows = np.zeros((n, dp), dtype=np.int64)66        for k in range(n):67            r = divmod_poly(tuple([0] * k + [1]), p)[1]68            rows[k, :len(r)] = r69        out[p] = ~((bits @ rows) % Q).any(axis=1)70    return out7172def factor_sets(primes, msk, n):73    small = [[] for _ in range(1 << n)]74    for p in primes:75        for i in np.nonzero(msk[p])[0]:76            if i:77                small[int(i)].append(p)78    ids = {p: k for k, p in enumerate(primes)}79    degs = [len(p) - 1 for p in primes]80    seen = {}81    sets = [()] * (1 << n)82    for i in range(1, 1 << n):83        h = poly_of(i, n)84        keys = []85        for p in small[i]:86            keys.append(ids[p])87            while True:88                q, r = divmod_poly(h, p)89                if r != ():90                    break91                h = q92        if len(h) > 1:93            if h not in seen:94                seen[h] = len(degs)95                degs.append(len(h) - 1)96            keys.append(seen[h])97        sets[i] = tuple(sorted(keys))98    return sets, degs99100def sieve_count(sets, degs, n, degcap=None):101    cnt = {}102    for i in range(1, 1 << n):103        ps = sets[i]104        if degcap is not None:105            ps = tuple(k for k in ps if degs[k] <= degcap)106        for k in range(len(ps) + 1):107            for sub in itertools.combinations(ps, k):108                cnt[sub] = cnt.get(sub, 0) + 1109    total = 0110    for key, m in cnt.items():111        total += (-1 if len(key) % 2 else 1) * (m * m + 2 * m)112    return total113114def brute(n):115    polys = [poly_of(i, n) for i in range(1 << n)]116    c = 0117    for a in polys:118        for b in polys:119            x, y = a, b120            while y != ():121                x, y = y, divmod_poly(x, y)[1]122            c += 1 if len(x) == 1 else 0123    return c124125primes = irreducibles(FACTOR_DEG)126msk = masks(primes, NMAX)127sets, degs = factor_sets(primes, msk, NMAX)128129print("DOMAIN q = 3, S = {0, 1}, degree below n, n up to %d, ordered pairs over all %d^2 including the zero polynomial"130      % (NMAX, 1 << NMAX))131print("irreducibles of degree at most %d: %d" % (FACTOR_DEG, len(primes)))132133for n in range(1, 5):134    assert sieve_count(sets, degs, n) == brute(n), n135print("cross-check against Euclid on every ordered pair, n = 1 to 4: passed")136137for n in LEVELS:138    z = sieve_count(sets, degs, n)139    print("n = %2d  Z = %d / %d  density = %.6f  gap from 9/16 = %+.6f"140          % (n, z, 4 ** n, z / 4 ** n, z / 4 ** n - 9 / 16))141142print("prediction (2/3) * (3/4) / (8/9) = 9/16 = %.6f" % (9 / 16))143144pi = {p: Fraction(int(msk[p][:1 << EULER_N].sum()), 1 << EULER_N)145      for p in primes if len(p) - 1 <= EULER_DEG}146for p in primes:147    if len(p) - 1 == 1:148        print("pi(%s) at n = %d = %s = %.6f" % (label(p), EULER_N, pi[p], float(pi[p])))149for d in (2, 3):150    ps = [p for p in primes if len(p) - 1 == d]151    devs = [abs(pi[p] - Fraction(1, Q ** d)) for p in ps]152    print("degree %d, %d primes, mean |pi - 3^-%d| = %.6f, max = %.6f"153          % (d, len(ps), d, float(sum(devs) / len(devs)), float(max(devs))))154155prod = Fraction(1)156for d in range(1, EULER_DEG + 1):157    for p in primes:158        if len(p) - 1 == d:159            prod *= 1 - pi[p] ** 2160    print("marginal Euler product through degree %d = %.6f" % (d, float(prod)))161162shared = sieve_count(sets, degs, EULER_N, degcap=EULER_DEG)163exact = Fraction(shared, 4 ** EULER_N)164print("exact no shared prime of degree at most %d at n = %d = %d/%d = %.6f"165      % (EULER_DEG, EULER_N, shared, 4 ** EULER_N, float(exact)))166print("marginal product minus exact = %+.6f" % float(prod - exact))167168g = log(4) / log(3)169print("gamma = log_3(4) = %.6f  gamma/2 = %.6f  window (gamma/2, 1/2] empty = %s"170      % (g, g / 2, g / 2 > 0.5))