parity_fill.py

1import time23import numpy as np4from fractions import Fraction5from math import gcd, sqrt67N_SPIN = 558ODDS_SPIN = list(range(1, N_SPIN + 1, 2))9L_SPIN = len(ODDS_SPIN)10GOLDEN = 137.50776411EXACT_TOP = 2112SPIN_R = 102413SWEEP_R = 51214CHECK_R = 409615LAYER_COUNTS = [4, 8, 14, 28]1617def lcm(a, b):18    return a // gcd(a, b) * b1920def odd_scales(top):21    return list(range(3, top + 1, 2))2223def mass_table(scales):24    period = 125    for n in scales:26        period = lcm(period, n)27    steps = [period // n for n in scales]28    counts = np.zeros(1 << len(scales), dtype=np.int64)29    block = 1 << 2030    for lo in range(0, period, block):31        j = np.arange(lo, min(lo + block, period), dtype=np.int64)32        code = np.zeros(j.shape, dtype=np.int64)33        for i, m in enumerate(steps):34            code |= ((j // m) & 1) << i35        counts += np.bincount(code, minlength=1 << len(scales))36    for i in range(len(scales)):37        bit = 1 << i38        for mask in range(1 << len(scales)):39            if not mask & bit:40                counts[mask] += counts[mask | bit]41    return period, counts4243def fill_exact(period, counts, scales, top):44    bits = [i for i, n in enumerate(scales) if n <= top]45    total = 046    for sub in range(1 << len(bits)):47        mask = 048        size = 049        for i, b in enumerate(bits):50            if sub >> i & 1:51                mask |= 1 << b52                size += 153        total += (-2) ** size * int(counts[mask]) ** 254    return Fraction(period * period - total, 2 * period * period)5556def indicator(n, period):57    j = np.arange(period, dtype=np.int64)58    return ((j // (period // n)) & 1).astype(bool)5960def fill_literal(top):61    scales = odd_scales(top)62    period = 163    for n in scales:64        period = lcm(period, n)65    acc = np.zeros((period, period), dtype=bool)66    for n in scales:67        c = indicator(n, period)68        acc ^= c[:, None] & c[None, :]69    return Fraction(int(acc.sum()), period * period)7071def fill_independent(top):72    prod = Fraction(1)73    for n in odd_scales(top):74        p = Fraction((n - 1) ** 2, 4 * n * n)75        prod *= 1 - 2 * p76    return (1 - prod) / 27778def fill_raster(top, r):79    u = (np.arange(r, dtype=np.float64) + 0.5) / r80    acc = np.zeros((r, r), dtype=bool)81    for n in odd_scales(top):82        c = parity(0.5 * n * u)83        acc ^= c[:, None] & c[None, :]84    return float(acc.mean())8586def primes_up_to_count(k):87    out = []88    c = 289    while len(out) < k:90        if all(c % d for d in range(2, int(sqrt(c)) + 1)):91            out.append(c)92        c += 193    return out9495def schedules():96    ps = primes_up_to_count(L_SPIN - 1)97    sched = {98        "unspun": [0.0] * L_SPIN,99        "degrees": [float(k + 1) for k in range(L_SPIN)],100        "primes": [0.0] + [float(p) for p in ps],101        "golden": [GOLDEN * (k + 1) for k in range(L_SPIN)],102    }103    for q in (2, 3, 4, 5):104        sched["eye 90/%d" % q] = [90.0 / q * k for k in range(L_SPIN)]105    return sched106107def parity(z):108    return z - np.floor(z) >= 0.5109110def spun_field(angles, r, layers=L_SPIN, dtype=np.float32):111    step = (np.arange(r, dtype=dtype) + dtype(0.5)) * dtype(1.0 / r)112    d = step - dtype(0.5)113    acc = np.zeros((r, r), dtype=bool)114    for k in range(layers):115        n = ODDS_SPIN[k]116        h = 0.5 * n117        t = np.float64(angles[k]) * np.pi / 180.0118        cx = dtype(h * np.cos(t))119        sx = dtype(h * np.sin(t))120        c0 = dtype(h * 0.5)121        a = c0 + cx * d[None, :] + sx * d[:, None]122        b = c0 - sx * d[None, :] + cx * d[:, None]123        acc ^= parity(a) & parity(b)124    return acc125126def disc_mask(r):127    d = (np.arange(r, dtype=np.float64) + 0.5) / r - 0.5128    return d[:, None] ** 2 + d[None, :] ** 2 <= 0.25129130def spun_fill(angles, r, layers=L_SPIN, disc=True):131    f = spun_field(angles, r, layers)132    return float(f[disc_mask(r)].mean()) if disc else float(f.mean())133134def raster_row(top, r):135    u = (np.arange(r, dtype=np.float64) + 0.5) / r136    acc = np.zeros((r, r), dtype=bool)137    out = {}138    for n in odd_scales(top):139        c = parity(0.5 * n * u)140        acc ^= c[:, None] & c[None, :]141        out[n] = float(acc.mean())142    return out143144def section_exact():145    scales = odd_scales(EXACT_TOP)146    period, counts = mass_table(scales)147    tops = list(range(3, EXACT_TOP + 1, 2))148    ras = raster_row(EXACT_TOP, CHECK_R)149    print("exact fill of the parity fold, odd scales 3..N, 1D grid period %d, %d subsets" % (period, 1 << len(scales)))150    exact = {}151    for top in tops:152        f = fill_exact(period, counts, scales, top)153        exact[top] = f154        lit = fill_literal(top) if top <= 9 else None155        print(156            " N %2d L %2d fill %s = %.9f raster %.9f gap %.2e" % (top, top // 2 + 1, f, float(f), ras[top], abs(float(f) - ras[top])),157            "" if lit is None else ("literal 2D XOR %s %s" % (lit, "match" if lit == f else "MISMATCH")),158        )159    zero = sum(1 for k in range(1 << len(scales)) if counts[k] == 0)160    print(" joint 1D masses: %d of the %d subsets have m_S = 0, the smallest being {3,5,7}" % (zero, 1 << len(scales)))161    print(" joint 1D masses m_S of the first scales", " ".join("m{%s} = %s" % (",".join(str(scales[i]) for i in range(len(scales)) if k >> i & 1), Fraction(int(counts[k]), period)) for k in (1, 2, 3, 5, 7)))162    return exact163164def section_independent(exact):165    print("exact against independent Bernoulli layers, p_n = ((n-1)/(2n))^2")166    for top in sorted(exact):167        f = exact[top]168        g = fill_independent(top)169        de = Fraction(1, 2) - f170        di = Fraction(1, 2) - g171        print(172            " N %2d exact %.9f independent %.9f difference %+.3e deviation %.6e against %.6e ratio %.6f"173            % (top, float(f), float(g), float(f - g), float(de), float(di), float(de / di))174        )175    print(" decay of the exact deviation, successive ratios", " ".join("%.6f" % float((Fraction(1, 2) - exact[t + 2]) / (Fraction(1, 2) - exact[t])) for t in sorted(exact)[:-1]))176    print(" decay of the independent deviation, successive ratios", " ".join("%.6f" % float((Fraction(1, 2) - fill_independent(t + 2)) / (Fraction(1, 2) - fill_independent(t))) for t in sorted(exact)[:-1]))177178def section_spun(exact):179    sched = schedules()180    print("spun parity fill at N %d, %d layers, R %d, inscribed disc" % (N_SPIN, L_SPIN, SPIN_R))181    rows = []182    for name, angles in sched.items():183        v = spun_fill(angles, SPIN_R)184        rows.append((name, v))185        print(" %-9s fill %.6f distance to 1/2 %.6f" % (name, v, abs(v - 0.5)))186    for name in ("unspun", "degrees", "eye 90/3"):187        print(" %-9s against resolution" % name, " ".join("R %d %.6f" % (r, spun_fill(sched[name], r)) for r in (256, 512, 1024, 2048)))188    best = min(rows, key=lambda r: abs(r[1] - 0.5))189    worst = max(rows, key=lambda r: abs(r[1] - 0.5))190    print(" closest to 1/2 %s at %.6f, furthest %s at %.6f" % (best[0], best[1], worst[0], worst[1]))191    print("fill against layer count at R %d, exact and square raster on the unit square, spun readings on the disc" % SPIN_R)192    for lc in LAYER_COUNTS:193        top = 2 * lc - 1194        ex = "%.6f" % float(exact[top]) if top in exact else "-"195        print(196            " L %2d N %2d exact %10s square %.6f disc" % (lc, top, ex, spun_fill(sched["unspun"], SPIN_R, lc, disc=False)),197            " ".join("%s %.6f" % (name, spun_fill(angles, SPIN_R, lc)) for name, angles in sched.items() if name in ("unspun", "degrees", "primes", "golden")),198        )199200def section_sweep():201    print("fixed increment sweep at N %d, R %d, disc fill against increment in whole degrees" % (N_SPIN, SWEEP_R))202    rows = []203    for deg in range(0, 91):204        angles = [float(deg) * k for k in range(L_SPIN)]205        rows.append((deg, spun_fill(angles, SWEEP_R)))206    for deg, v in rows:207        print(" increment %2d fill %.6f" % (deg, v))208    order = sorted(rows, key=lambda r: r[1])209    print(" lowest", " ".join("%d:%.6f" % r for r in order[:5]))210    print(" highest", " ".join("%d:%.6f" % r for r in order[-5:]))211    near = sorted(rows, key=lambda r: abs(r[1] - 0.5))212    print(" closest to 1/2", " ".join("%d:%.6f" % r for r in near[:5]))213    print(" eyes q = 1, 5, 3, 2 at increments 0, 18, 30, 45 and the quarter turn", " ".join("%d:%.6f" % (d, dict(rows)[d]) for d in (0, 18, 30, 45, 90)))214    v = dict(rows)215    print(" mirror check fill(d) = fill(90 - d): %d of 46 pairs equal to the bit, largest gap %.3e, quarter-turn check fill(90) = fill(0): %s" % (sum(1 for d in range(46) if v[d] == v[90 - d]), max(abs(v[d] - v[90 - d]) for d in range(46)), v[0] == v[90]))216    print(" spread of the 89 nonzero whole increments: min %.6f max %.6f, unspun %.6f" % (min(v[d] for d in range(1, 90)), max(v[d] for d in range(1, 90)), v[0]))217218def main():219    t0 = time.perf_counter()220    exact = section_exact()221    section_independent(exact)222    section_spun(exact)223    section_sweep()224    print("seconds", round(time.perf_counter() - t0, 2))225226main()