sponge_tube.py

1import math2import time3from fractions import Fraction as Fr45import numpy as np6from mpmath import iv, mp78mp.prec = 133910iv.prec = 13311D = iv.log(20) / iv.log(3)12WEIGHT = iv.mpf(27) / 2013CELLW = (3, 2, 3)14LEVELS = 40015DEEP_LEVELS = 716TAIL_LEVELS = 80171819def num(x):20    return iv.mpf(x.numerator) / x.denominator if isinstance(x, Fr) else iv.mpf(x)212223def asin(u):24    return iv.atan2(u, iv.sqrt(1 - u * u))252627def a0(t, d):28    t = min(t, d)29    x, dd = num(t), num(d)30    if t == d:31        return iv.pi * dd * dd / 432    return x / 2 * iv.sqrt(dd * dd - x * x) + dd * dd / 2 * asin(x / dd)333435def a1(t, d):36    t = min(t, d)37    x, dd = num(t), num(d)38    if t == d:39        return dd**3 / 340    y = dd * dd - x * x41    return x * x * (3 * dd**4 - 3 * dd * dd * x * x + x**4) / (3 * (dd**3 + y * iv.sqrt(y)))424344def seg(t0, t1, alpha, beta, d):45    if t1 <= t0:46        return iv.mpf(0)47    return num(alpha) * (a0(t1, d) - a0(t0, d)) - num(beta) * (a1(t1, d) - a1(t0, d))484950def hole_j(s, d):51    return 4 * seg(Fr(0), s / 2, s, 2, d)525354def partial_j(s, c, d):55    left = seg(Fr(0), min(s / 2, c), s, 2, d)56    right = seg(max(Fr(0), s - c), s / 2, s, 2, d)57    if c <= s / 2:58        bottom = seg(Fr(0), c, c, 1, d)59    else:60        bottom = seg(Fr(0), s - c, c, 1, d) + seg(s - c, s / 2, s, 2, d)61    return left + right + 2 * bottom626364def column_weight(i, digits):65    w = 166    for _ in range(digits):67        w *= CELLW[i % 3]68        i //= 369    return w707172def columns_below(count, digits):73    if count >= 3**digits:74        return 8**digits75    total, prefix = 0, 176    for d in range(digits - 1, -1, -1):77        digit = (count // 3**d) % 378        total += prefix * sum(CELLW[:digit]) * 8**d79        prefix *= CELLW[digit]80    return total818283def wall_half(delta):84    total = iv.mpf(0)85    for m in range(1, LEVELS + 1):86        total += 8 ** (m - 1) * hole_j(Fr(1, 3 ** (m + 1)), delta)87    total /= 288    assert upper(total) <= lower(num(delta / 18))89    return total, num(delta * Fr(8, 9) ** LEVELS / 16)909192def strip(delta):93    total = iv.mpf(0)94    for m in range(1, LEVELS + 1):95        s = Fr(1, 3 ** (m + 1))96        ratio = delta / s97        full = int((ratio - 2) // 3) + 198        if full > 0:99            total += columns_below(min(full, 3 ** (m - 1)), m - 1) * hole_j(s, delta)100        cut = int((ratio - 1) // 3)101        if 0 <= cut < 3 ** (m - 1) and (3 * cut + 1) * s < delta < (3 * cut + 2) * s:102            total += column_weight(cut, m - 1) * partial_j(s, delta - (3 * cut + 1) * s, delta)103    assert upper(total) <= lower(num(delta * delta / 3))104    return total, num(delta * Fr(8, 9) ** LEVELS / 8)105106107def hole_rows(i, digits):108    rows = [0]109    for d in range(digits):110        digit = (i // 3**d) % 3111        allowed = (0, 2) if digit == 1 else (0, 1, 2)112        rows = [r + a * 3**d for r in rows for a in allowed]113    return sorted(rows)114115116def overlap(ints_a, ints_b):117    total, j = Fr(0), 0118    for lo, hi in ints_a:119        while j < len(ints_b) and ints_b[j][1] <= lo:120            j += 1121        k = j122        while k < len(ints_b) and ints_b[k][0] < hi:123            total += min(hi, ints_b[k][1]) - max(lo, ints_b[k][0])124            k += 1125    return total126127128def columns_meeting(lo, hi, s, digits):129    first = max(0, int((lo / s - 2) // 3))130    last = min(3**digits - 1, int((hi / s - 1) // 3) + 1)131    return [i for i in range(first, last + 1) if (3 * i + 2) * s > lo and (3 * i + 1) * s < hi]132133134def deep_bound(delta):135    total = Fr(0)136    cache = {}137    for m in range(1, DEEP_LEVELS + 1):138        sm = Fr(1, 3 ** (m + 1))139        wm = sm * sm / (4 * delta)140        for n in range(1, DEEP_LEVELS + 1):141            sn = Fr(1, 3 ** (n + 1))142            wn = sn * sn / (4 * delta)143            for i in columns_meeting(delta - wn, delta, sm, m - 1):144                vlen = min(delta, (3 * i + 2) * sm) - max(delta - wn, (3 * i + 1) * sm)145                if (m, i) not in cache:146                    cache[(m, i)] = [((3 * r + 1) * sm, (3 * r + 2) * sm) for r in hole_rows(i, m - 1)]147                for k in columns_meeting(delta - wm, delta, sn, n - 1):148                    ulen = min(delta, (3 * k + 2) * sn) - max(delta - wm, (3 * k + 1) * sn)149                    if (n, k) not in cache:150                        cache[(n, k)] = [((3 * r + 1) * sn, (3 * r + 2) * sn) for r in hole_rows(k, n - 1)]151                    total += vlen * ulen * overlap(cache[(m, i)], cache[(n, k)])152    for m in range(1, TAIL_LEVELS + 1):153        for n in range(1, TAIL_LEVELS + 1):154            if max(m, n) <= DEEP_LEVELS:155                continue156            wm = min(Fr(1, 3 ** (2 * m + 2)) / (4 * delta), delta)157            wn = min(Fr(1, 3 ** (2 * n + 2)) / (4 * delta), delta)158            count = min(math.ceil(wn * 3**m) + 2, math.ceil(wm * 3**n) + 2)159            total += wm * wn * min(Fr(count, 9), Fr(1, 3))160    total += 2 * Fr(1, 3 ** (2 * TAIL_LEVELS + 6)) / (36 * delta * delta) / Fr(8, 9) / 8161    return total162163164def tube(delta):165    d = num(delta)166    half, half_tail = wall_half(delta)167    strip_v, strip_tail = strip(delta)168    value = iv.pi * d * d - 8 * iv.sqrt(2) * d**3 + 8 * d * d + 48 * (half - strip_v)169    return value - 48 * strip_tail - 24 * num(deep_bound(delta)), value + 48 * half_tail170171172def tube_upper(delta):173    s_sum = Fr(0)174    for m in range(1, LEVELS + 1):175        s = Fr(1, 3 ** (m + 1))176        s_sum += 8 ** (m - 1) * s * delta * min(s, 4 * delta)177    s_sum += delta * Fr(8, 9) ** LEVELS / 8178    return iv.pi * num(delta) ** 2 + 24 * num(s_sum)179180181def periodic(eps, depth):182    lo = hi = iv.mpf(20) / 27183    for l in range(depth + 1):184        tl, th = tube(eps / 3**l)185        lo += WEIGHT**l * tl186        hi += WEIGHT**l * th187    hi += WEIGHT ** (depth + 1) * tube_upper(eps / 3 ** (depth + 1)) * iv.mpf(20) / 9188    scale = iv.exp((D - 3) * iv.log(num(eps)))189    return lower(lo * scale), upper(hi * scale)190191192def carpet_distance(v, z, side, depth=40):193    v, z = v / side, z / side194    s = np.full(v.shape, side)195    out = np.zeros(v.shape)196    done = np.zeros(v.shape, dtype=bool)197    for _ in range(depth):198        s = s / 3199        dv, dz = np.floor(v * 3).astype(int), np.floor(z * 3).astype(int)200        hole = (dv == 1) & (dz == 1) & ~done201        rel = np.minimum.reduce([v * 3 - 1, 2 - v * 3, z * 3 - 1, 2 - z * 3])202        out[hole] = s[hole] * rel[hole]203        done |= hole204        v, z = v * 3 - dv, z * 3 - dz205    return out206207208def raster(delta, n):209    h = 1 / (3 * n)210    grid = (np.arange(n) + 0.5) * h211    u, v, z = np.meshgrid(grid, grid, grid, indexing="ij")212    dist = np.full(u.shape, np.inf)213    dvz, duz = carpet_distance(v, z, 1 / 3), carpet_distance(u, z, 1 / 3)214    for t, dd in ((u, dvz), (1 / 3 - u, dvz), (v, duz), (1 / 3 - v, duz)):215        dist = np.minimum(dist, np.sqrt(t * t + dd * dd))216    arm = dist217    e = np.minimum.reduce([np.sqrt(np.minimum(u, 1 / 3 - u) ** 2 + np.minimum(v, 1 / 3 - v) ** 2), np.sqrt(np.minimum(u, 1 / 3 - u) ** 2 + np.minimum(z, 1 / 3 - z) ** 2), np.sqrt(np.minimum(v, 1 / 3 - v) ** 2 + np.minimum(z, 1 / 3 - z) ** 2)])218    r = np.sqrt(3) * h / 2219    inner = 6 * (arm <= delta - r).sum() + (e <= delta - r).sum()220    outer = 6 * (arm <= delta + r).sum() + (e <= delta + r).sum()221    return inner * h**3, outer * h**3222223224def lower(x):225    return mp.make_mpf(x._mpi_[0]) if hasattr(x, "_mpi_") else x226227228def upper(x):229    return mp.make_mpf(x._mpi_[1]) if hasattr(x, "_mpi_") else x230231232def down(x, k):233    return mp.mpf(int(mp.floor(lower(x) * 10**k))) / 10**k234235236def up(x, k):237    return mp.mpf(int(mp.ceil(upper(x) * 10**k))) / 10**k238239240def main():241    start = time.time()242    print(f"SPONGE: base 3, fill 20, D = log 20 / log 3 in [{mp.nstr(lower(D), 12)}, {mp.nstr(upper(D), 12)}], weight per level 27/20, constant term 20/27")243    print("TUBE: T(delta) = volume of F_delta inside the plus, delta in (0, 1/6], from the wall carpets and the centre cube edges, interval arithmetic at 133 bits")244    for delta in (Fr(1, 8), Fr(1, 12)):245        lo, hi = tube(delta)246        lo, hi = lower(lo), upper(hi)247        rl, rh = raster(float(delta), 120)248        print(f"  delta = {delta}: closed form in [{mp.nstr(down(lo, 9), 10)}, {mp.nstr(up(hi, 9), 10)}], raster band [{rl:.5f}, {rh:.5f}] at 120^3 cells per cube")249        assert rl <= float(lo) and float(hi) <= rh250    lo, hi = tube(Fr(1, 6))251    ident = (iv.pi + 8) / 36 - iv.sqrt(2) / 27252    assert lower(lo) <= upper(ident) and lower(ident) <= upper(hi)253    print(f"  delta = 1/6: A1 = V1 there, so T = (pi + 8)/36 - sqrt2/27 - 24 Deep = {mp.nstr(down(ident, 9), 10)} less at most {float(24 * deep_bound(Fr(1, 6))):.2e}; closed form in [{mp.nstr(down(lo, 9), 10)}, {mp.nstr(up(hi, 9), 10)}]")254    print("PERIODIC: p(eps) = eps^(D-3) (20/27 + sum_l (27/20)^l T(eps/3^l)), eps in (sqrt2/18, 1/6]")255    bands = {}256    for eps in (Fr(1, 12), Fr(1, 8), Fr(1, 6)):257        lo, hi = periodic(eps, 40)258        bands[eps] = (lo, hi)259        print(f"  eps = {eps}: p in [{mp.nstr(down(lo, 6), 7)}, {mp.nstr(up(hi, 6), 7)}]")260    gap = bands[Fr(1, 6)][0] - bands[Fr(1, 12)][1]261    gap2 = bands[Fr(1, 8)][0] - bands[Fr(1, 12)][1]262    assert gap > 0 and gap2 > 0263    print(f"  p(1/6) - p(1/12) >= {mp.nstr(down(gap, 6), 7)},  p(1/8) - p(1/12) >= {mp.nstr(down(gap2, 6), 7)}")264    swing = gap / bands[Fr(1, 12)][1]265    print(f"  relative swing (p(1/6) - p(1/12)) / p(1/12) >= {mp.nstr(down(swing * 100, 4), 5)} %")266    print(f"  doubly deep set, upper bound on its volume: {float(deep_bound(Fr(1, 6))):.3g} at delta = 1/6, {float(deep_bound(Fr(1, 8))):.3g} at delta = 1/8, {float(deep_bound(Fr(1, 12))):.3g} at delta = 1/12")267    print(f"wall {time.time() - start:.1f} s")268269270if __name__ == "__main__":271    main()