carpet_tube.py

1import math2import time3from fractions import Fraction45import numpy as np6from mpmath import mp, mpf7from scipy.ndimage import distance_transform_edt89mp.dps = 5010BASE = 311FILL = 812DIM = mp.log(FILL) / mp.log(BASE)13DIM_F = float(DIM)14LOW = (Fraction(1), Fraction(4, 5), Fraction(-4, 7))15HIGH = (Fraction(9, 8), Fraction(3, 10), Fraction(-1, 14))16SEAM = Fraction(379, 280)171819def hole_tube(side, eps):20    return side * side if side <= 2 * eps else 4 * eps * side - 4 * eps * eps212223def split_index(eps):24    m = 125    while Fraction(1, BASE**m) > 2 * eps:26        m += 127    return m282930def direct_sum(eps):31    m0 = split_index(eps)32    ring = sum(FILL ** (m - 1) * hole_tube(Fraction(1, BASE**m), eps) for m in range(1, m0))33    return ring + Fraction(FILL, BASE**2) ** (m0 - 1)343536def closed_v(eps):37    m0 = split_index(eps)38    a = Fraction(4, 5) * eps * (Fraction(FILL, BASE) ** (m0 - 1) - 1)39    b = Fraction(4, 7) * eps * eps * (FILL ** (m0 - 1) - 1)40    return a - b + Fraction(FILL, BASE**2) ** (m0 - 1)414243def branch(t):44    return LOW if t < Fraction(1, 2) else HIGH454647def to_mpf(x):48    return mpf(x.numerator) / x.denominator if isinstance(x, Fraction) else mpf(x)495051def profile(t):52    x = to_mpf(t)53    c0, c1, c2 = branch(t if isinstance(t, Fraction) else Fraction(float(x)))54    return x ** (DIM - 2) * (to_mpf(c0) + to_mpf(c1) * x + to_mpf(c2) * x * x)555657def profile_float(t):58    c0, c1, c2 = branch(Fraction(t))59    return t ** (DIM_F - 2) * (float(c0) + float(c1) * t + float(c2) * t * t)606162def stationary(coeffs):63    c0, c1, c2 = (to_mpf(c) for c in coeffs)64    a = -c2 * DIM65    b = -c1 * (DIM - 1)66    c = -c0 * (DIM - 2)67    disc = mp.sqrt(b * b - 4 * a * c)68    return sorted([(-b - disc) / (2 * a), (-b + disc) / (2 * a)])697071def filled_cells(level):72    cells = np.ones((1, 1), dtype=bool)73    for _ in range(level):74        block = np.ones((BASE, BASE), dtype=bool)75        block[1, 1] = False76        cells = np.kron(cells, block)77    return cells787980def hole_census(level):81    cells = filled_cells(level)82    side = BASE**level83    padded = np.pad(cells, 1, constant_values=True)84    counted = 085    for m in range(1, level + 1):86        hole = BASE ** (level - m)87        parents = filled_cells(m - 1)88        origins = np.argwhere(parents) * (BASE * hole) + hole89        assert len(origins) == FILL ** (m - 1)90        for i, j in origins:91            assert not cells[i : i + hole, j : j + hole].any()92            ring = padded[i : i + hole + 2, j : j + hole + 2].copy()93            ring[1:-1, 1:-1] = True94            assert ring.all()95        counted += len(origins) * hole * hole96    assert counted == side * side - FILL**level97    return level, FILL**level, side * side - FILL**level9899100def raster_tube(level, cells_in):101    cells = filled_cells(level)102    side = BASE**level103    dist = distance_transform_edt(~cells)104    count = cells.sum() + ((dist <= cells_in + 1e-9) & ~cells).sum()105    return Fraction(int(count), side * side)106107108def window(ulo, uhi, samples):109    worst = 0.0110    lo, hi = math.inf, -math.inf111    for u in np.linspace(ulo, uhi, samples):112        eps = math.exp(-u)113        n = math.floor(-math.log(eps) / math.log(BASE))114        t = eps * BASE**n115        if t >= 1:116            n -= 1117            t = eps * BASE**n118        if t < 1 / BASE:119            n += 1120            t = eps * BASE**n121        m0 = split_index(Fraction(eps))122        ring = sum(FILL ** (m - 1) * float(hole_tube(Fraction(1, BASE**m), Fraction(eps))) for m in range(1, m0))123        v = ring + (FILL / BASE**2) ** (m0 - 1)124        m = eps ** (DIM_F - 2) * v125        worst = max(worst, abs(m - profile_float(t)))126        lo, hi = min(lo, m), max(hi, m)127    return lo, hi, worst128129130def main():131    start = time.time()132    print(f"CARPET: base {BASE}, fill {FILL}, d = log {FILL} / log {BASE} = {mp.nstr(DIM, 12)}")133    print("HOLES: every hole is an open square whose boundary lies in the carpet, and the holes exhaust the complement")134    for level in (3, 4, 5):135        lv, filled, empty = hole_census(level)136        print(f"  level {lv}: {filled} filled cells, {empty} empty cells in 8^(m-1) holes of side 3^-m, m = 1..{lv}, every hole ringed by filled cells")137    print()138    print("TUBE: V(eps) = sum_m 8^(m-1) h(3^-m, eps), h(s, eps) = s^2 if s <= 2 eps else 4 eps s - 4 eps^2, against the closed form, exact rationals")139    grid = [Fraction(1, 3), Fraction(2, 5), Fraction(1, 2), Fraction(3, 5), Fraction(2, 3), Fraction(4, 5), Fraction(9, 10), Fraction(999, 1000)]140    checks = 0141    for n in range(1, 13):142        for t in grid:143            eps = t / BASE**n144            assert direct_sum(eps) == closed_v(eps)145            checks += 1146    print(f"  {checks} pairs (n, t), n = 1..12, t in {{1/3, 2/5, 1/2, 3/5, 2/3, 4/5, 9/10, 999/1000}}: direct sum equals closed form exactly")147    print()148    print("PROFILE: eps = 3^-n t, t in [1/3, 1), M(eps) = eps^(d-2) V(eps) -> G(t)")149    print("  G(t) = t^(d-2) (1 + 4t/5 - 4t^2/7) on [1/3, 1/2),  G(t) = t^(d-2) (9/8 + 3t/10 - t^2/14) on [1/2, 1)")150    half = Fraction(1, 2)151    left = sum(c * half**i for i, c in enumerate(LOW))152    right = sum(c * half**i for i, c in enumerate(HIGH))153    assert left == right == Fraction(44, 35)154    seam_low = Fraction(1, FILL) * sum(c * Fraction(1, 3) ** i for i, c in enumerate(LOW)) * 9155    seam_high = sum(HIGH)156    assert seam_low == seam_high == SEAM157    slope_low = sum(i * c * half ** (i - 1) for i, c in enumerate(LOW) if i)158    slope_high = sum(i * c * half ** (i - 1) for i, c in enumerate(HIGH) if i)159    assert slope_low == slope_high == Fraction(8, 35)160    print(f"  continuous at t = 1/2: both branches give 44/35 * 2^(2-d) = {mp.nstr(profile(half), 12)}, and C^1 there: both polynomial factors have slope 8/35")161    print(f"  periodic: G(1/3) = G(1) = 379/280 = {mp.nstr(mpf(379) / 280, 12)}")162    tmax = [r for r in stationary(LOW) if Fraction(1, 3) <= Fraction(float(r)) < half]163    tmin = [r for r in stationary(HIGH) if half <= Fraction(float(r)) < 1]164    assert len(tmax) == 1 and len(tmin) == 1165    tmax, tmin = tmax[0], tmin[0]166    gmax, gmin = profile(tmax), profile(tmin)167    scan = [profile_float(t) for t in np.linspace(1 / 3, 1, 200001)]168    assert min(scan) >= float(gmin) - 1e-12 and max(scan) <= float(gmax) + 1e-12169    swing = (gmax - gmin) / gmin170    scale = mpf(10) ** 11171    print(f"  maximum {mp.nstr(gmax, 15)}, safe {mp.nstr(mp.ceil(gmax * scale) / scale, 12)} up, at t = {mp.nstr(tmax, 9)} (root of (4/7) d t^2 - (4/5)(d-1) t - (d-2) = 0 in [1/3, 1/2))")172    print(f"  minimum {mp.nstr(gmin, 15)}, safe {mp.nstr(mp.floor(gmin * scale) / scale, 12)} down, at t = {mp.nstr(tmin, 9)} (root of (d/14) t^2 - (3/10)(d-1) t - (9/8)(d-2) = 0 in [1/2, 1))")173    print(f"  swing (max - min) / min = {mp.nstr(swing * 100, 8)} %   scan of 200001 points stays inside [min, max]")174    print(f"  Cantor design for scale: swing 3.53 %; the carpet swings {mp.nstr(swing * 100, 4)} %")175    print()176    print("MEASURED: the level sum at eps = e^-u against G(t)")177    lo, hi, worst = window(50.0, 60.0, 2001)178    print(f"  u in [50, 60]: min {lo:.12f}   max {hi:.12f}   max |M(eps) - G(t)| = {worst:.1e}")179    assert worst < 1e-9180    eps = math.exp(-60)181    print(f"  outer collar at u = 60: eps^(D-2) (4 eps + pi eps^2) = {eps ** (DIM_F - 2) * (4 * eps + math.pi * eps * eps):.1e}")182    print()183    print("RASTER: level-6 cells, Euclidean distance transform to the nearest filled cell centre, no hole lemma used")184    level, k = 6, 21185    side = BASE**level186    raster = raster_tube(level, k)187    exact = closed_v(Fraction(k, side))188    print(f"  cells within {k} cell widths of a filled centre, plus the filled cells: {raster} = {float(raster):.12f}")189    print(f"  exact V({k}/{side}) from the closed form: {exact} = {float(exact):.12f}")190    assert raster == exact191    print("  equal as rationals")192    print()193    print(f"wall {time.time() - start:.1f} s")194195196if __name__ == "__main__":197    main()