radix.py

1import itertools2from math import log3import numpy as np45N = 126SMALL = (2, 3, 5, 7, 11, 13)7SCHEDULES = [8    ("alternating base 2 {0,1} / base 3 {0,1}", ((2, (0, 1)), (3, (0, 1)))),9    ("alternating base 2 {0,1} / base 3 {0,2}", ((2, (0, 1)), (3, (0, 2)))),10    ("pure base 2 {0,1}", ((2, (0, 1)),)),11    ("pure base 3 {0,1}", ((3, (0, 1)),)),12    ("pure base 3 {0,2}", ((3, (0, 2)),)),13]1415def build(schedule, n):16    vals = [0]17    place = 118    for j in range(n):19        base, digits = schedule[j % len(schedule)]20        vals = [v + d * place for v in vals for d in digits]21        place *= base22    return sorted(vals)2324def moebius_table(n):25    mu = np.ones(n + 1, dtype=np.int64)26    prime = np.ones(n + 1, dtype=bool)27    prime[:2] = False28    for p in range(2, int(n ** 0.5) + 1):29        if prime[p]:30            prime[p * p::p] = False31    for p in np.nonzero(prime)[0]:32        mu[p::p] = -mu[p::p]33        mu[p * p::p * p] = 034    mu[0] = 035    return mu3637def present(vals, top):38    a = np.zeros(top + 1, dtype=np.int64)39    a[np.array(vals, dtype=np.int64)] = 140    a[0] = 041    return a4243def coprime_pairs(vals, mu, top):44    a = present(vals, top)45    ordered = 046    for m in range(1, top + 1):47        if mu[m]:48            c = int(a[m::m].sum())49            if c:50                ordered += int(mu[m]) * c * c51    unit = 1 if 1 in vals else 052    diagonal = unit53    with_zero = unit54    return (ordered - diagonal) // 2 + with_zero5556def divisible_fraction(vals, m):57    return sum(1 for v in vals if v % m == 0) / len(vals)5859def joint_small(vals):60    arr = np.array(vals, dtype=np.int64)61    n = len(vals)62    total = 063    for k in range(len(SMALL) + 1):64        for sub in itertools.combinations(SMALL, k):65            d = 166            for p in sub:67                d *= p68            c = int((arr % d == 0).sum())69            total += (-1 if k % 2 else 1) * c * (c - 1) // 270    return total / (n * (n - 1) // 2)7172top = max(max(build(s, N)) for _, s in SCHEDULES)73mu = moebius_table(top)7475print("DOMAIN %d digit positions, unordered distinct pairs including the zero point" % N)76print("points per schedule %d, pairs per schedule %d" % (1 << N, (1 << N) * ((1 << N) - 1) // 2))7778rows = {}79for name, sched in SCHEDULES:80    vals = build(sched, N)81    pairs = len(vals) * (len(vals) - 1) // 282    c = coprime_pairs(vals, mu, max(vals))83    rows[name] = (vals, c, c / pairs)84    print("%-40s coprime %7d  density %.6f" % (name, c, c / pairs))8586a, b = SCHEDULES[0][0], SCHEDULES[1][0]87print("alternating gap = %.6f" % (rows[b][2] - rows[a][2]))88print("pure base 2 against 1/zeta(2) = 0.607927, difference %+.6f" % (rows[SCHEDULES[2][0]][2] - 0.6079271018540267))8990vals = rows[a][0]91res = [sum(1 for v in vals if v % 6 == r) for r in range(6)]92print("residues mod 6 of the first alternating schedule: %s" % (", ".join(str(x) for x in res)))93for name in (a, b):94    vs = rows[name][0]95    print("%-40s even fraction %.6f  divisible by 3 fraction %.6f"96          % (name, divisible_fraction(vs, 2), divisible_fraction(vs, 3)))9798for name in (a, b):99    vs = rows[name][0]100    prod = 1.0101    for p in SMALL:102        prod *= 1 - divisible_fraction(vs, p) ** 2103    exact = joint_small(vs)104    print("%-40s Euler product through 13 = %.6f  exact joint = %.6f  product minus exact = %+.6f"105          % (name, prod, exact, prod - exact))106107for p in SMALL:108    fa = divisible_fraction(rows[a][0], p)109    fb = divisible_fraction(rows[b][0], p)110    ra = 1 - fa * fa111    rb = 1 - fb * fb112    print("prime %2d  factor %.6f -> %.6f  log advantage %+.6f" % (p, ra, rb, log(rb / ra)))