checks.py
17.7 kB · python · 451 lines
1import os2import sys3import time4from fractions import Fraction as Fr56import numpy as np7from mpmath import asin, exp, log, mp, mpf, pi, quad, sqrt8from scipy.spatial import cKDTree910HERE = os.path.dirname(os.path.abspath(__file__))11sys.path.insert(0, HERE)12import sponge_tube as st1314LEVEL = 415SIDE = 3**LEVEL16REMOVED = np.array([(1, 1, 1), (0, 1, 1), (2, 1, 1), (1, 0, 1), (1, 2, 1), (1, 1, 0), (1, 1, 2)])17PLUS = 7 / 2718SQRT2 = np.sqrt(2)1920# THE LEVEL-4 PREFRACTAL212223def prefractal():24 idx = np.arange(SIDE)25 digits = np.stack([(idx // 3**k) % 3 for k in range(LEVEL)], 1)26 one = digits == 127 i, j, k = np.meshgrid(idx, idx, idx, indexing="ij")28 count = one[i].astype(int) + one[j].astype(int) + one[k].astype(int)29 keep = (count <= 1).all(axis=-1)30 return np.stack([i[keep], j[keep], k[keep]], 1)313233class Oracle:34 def __init__(self, cubes):35 h = 1.0 / SIDE36 self.lo = cubes * h37 self.hi = self.lo + h38 self.tree = cKDTree(self.lo + h / 2)39 self.halfdiag = np.sqrt(3) * h / 24041 def dist(self, points, k=300):42 d, ii = self.tree.query(points, k=k)43 q = points[:, None, :]44 g = np.maximum(np.maximum(self.lo[ii] - q, q - self.hi[ii]), 0)45 best = np.sqrt((g * g).sum(-1)).min(1)46 complete = d[:, -1] - self.halfdiag > best47 return best, complete484950def wall_faces(cubes):51 n = SIDE52 faces = []53 for a in range(3):54 for arm in (0, 2):55 for b in ((a + 1) % 3, (a + 2) % 3):56 c = 3 - a - b57 for side in (0, 1):58 plane = 1 / 3 if side == 0 else 2 / 359 sel = (cubes[:, b] == n // 3 - 1) if side == 0 else (cubes[:, b] == 2 * n // 3)60 sel &= (cubes[:, a] >= arm * n // 3) & (cubes[:, a] < (arm + 1) * n // 3)61 sel &= (cubes[:, c] >= n // 3) & (cubes[:, c] < 2 * n // 3)62 faces.append((b, plane, sel))63 return faces646566def dist_to_faces(points, oracle, faces, own):67 c = np.floor(points * 3).astype(int)68 centre = (c == 1).all(1)69 best = np.full(len(points), np.inf)70 for b, plane, sel in faces:71 lo, hi = oracle.lo[sel], oracle.hi[sel]72 others = [k for k in range(3) if k != b]73 q = points[:, None, :]74 g = np.maximum(np.maximum(lo[None, :, others] - q[:, :, others], q[:, :, others] - hi[None, :, others]), 0)75 dd = np.sqrt((g * g).sum(-1) + (points[:, None, b] - plane) ** 2).min(1)76 if own:77 inside = np.ones(len(points), bool)78 for k in others:79 inside &= (points[:, k] >= lo[:, k].min() - 1e-12) & (points[:, k] <= hi[:, k].max() + 1e-12)80 inside &= np.abs(points[:, b] - plane) <= 1 / 3 + 1e-1281 best = np.where(inside & ~centre, np.minimum(best, dd), best)82 else:83 best = np.minimum(best, dd)84 if own:85 m = np.minimum(points[centre] - 1 / 3, 2 / 3 - points[centre])86 best[centre] = np.minimum.reduce([np.hypot(m[:, 0], m[:, 1]), np.hypot(m[:, 0], m[:, 2]), np.hypot(m[:, 1], m[:, 2])])87 return best888990# THE REDUCED DISTANCE OF THE LEMMAS919293def carpet_dist(v, z, side, depth=40):94 v, z = v / side, z / side95 s = np.full(v.shape, side)96 out = np.zeros(v.shape)97 done = np.zeros(v.shape, dtype=bool)98 for _ in range(depth):99 s = s / 3100 dv, dz = np.floor(v * 3).astype(int), np.floor(z * 3).astype(int)101 hole = (dv == 1) & (dz == 1) & ~done102 rel = np.minimum.reduce([v * 3 - 1, 2 - v * 3, z * 3 - 1, 2 - z * 3])103 out[hole] = s[hole] * rel[hole]104 done |= hole105 v, z = v * 3 - dv, z * 3 - dz106 return out107108109def reduced(points):110 c = np.floor(points * 3).astype(int)111 centre = (c == 1).all(1)112 out = np.full(len(points), np.inf)113 for a in range(3):114 arm = (c[:, a] != 1) & (c[:, (a + 1) % 3] == 1) & (c[:, (a + 2) % 3] == 1)115 along = points[arm, a] - c[arm, a] / 3116 for b in ((a + 1) % 3, (a + 2) % 3):117 other = (a + 2) % 3 if b == (a + 1) % 3 else (a + 1) % 3118 u = points[arm, b] - 1 / 3119 dd = carpet_dist(points[arm, other] - 1 / 3, along, 1 / 3)120 out[arm] = np.minimum(out[arm], np.sqrt(u * u + dd * dd))121 out[arm] = np.minimum(out[arm], np.sqrt((1 / 3 - u) ** 2 + dd * dd))122 m = np.minimum(points[centre] - 1 / 3, 2 / 3 - points[centre])123 out[centre] = np.minimum.reduce([np.hypot(m[:, 0], m[:, 1]), np.hypot(m[:, 0], m[:, 2]), np.hypot(m[:, 1], m[:, 2])])124 return out125126127def sample_plus(rng, n):128 base = REMOVED[rng.integers(0, 7, n)] / 3.0129 return base + rng.random((n, 3)) / 3.0130131132def three_sigma(hits, n, scale):133 p = hits / n134 return p * scale, 3 * np.sqrt(p * (1 - p) / n) * scale135136137# THE ORACLE VERB138139140def verb_oracle():141 t0 = time.time()142 rng = np.random.default_rng(7)143 cubes = prefractal()144 oracle = Oracle(cubes)145 faces = wall_faces(cubes)146 print(f"ORACLE: level-{LEVEL} prefractal, {len(cubes)} cubes, {sum(int(f[2].sum()) for f in faces)} cube faces on the 24 walls")147 points = sample_plus(rng, 20000)148 near = sample_plus(rng, 4000)149 axis = rng.integers(0, 3, 4000)150 near[np.arange(4000), axis] = np.where(rng.random(4000) < 0.5, 1 / 3, 2 / 3) + rng.normal(0, 1e-3, 4000)151 near = np.clip(near, 0, 1)152 c = np.minimum(np.floor(near * 3).astype(int), 2)153 near = near[(c == 1).sum(1) >= 2]154 points = np.vstack([points, near])155 d_f4, complete = oracle.dist(points)156 d_walls = dist_to_faces(points, oracle, faces, own=False)157 d_own = dist_to_faces(points, oracle, faces, own=True)158 print(f" {len(points)} points of the plus, {len(near)} of them within about 1e-3 of a face; candidate set complete for every point: {bool(complete.all())}")159 print(f" max |dist(x, F_4) - dist(x, 24 level-4 wall carpets)| = {np.abs(d_f4 - d_walls).max():.1e}")160 print(f" max |dist(x, F_4) - own four walls or centre edges at level 4| = {np.abs(d_f4 - d_own).max():.1e}")161 d_red = reduced(points)162 print(f" reduced distance of the lemmas (infinite carpets) minus dist(x, F_4): min {(d_red - d_f4).min():.1e}, max {(d_red - d_f4).max():.2e}, against sqrt(2)/{3**LEVEL} = {SQRT2 / 3**LEVEL:.2e}")163 centre = oracle.dist(np.array([[0.5, 0.5, 0.5]]))[0][0]164 print(f" dist(centre, F_4) = {centre:.17f}, sqrt(2)/6 = {SQRT2 / 6:.17f}, difference {abs(centre - SQRT2 / 6):.1e}")165 print(f" wall {time.time() - t0:.1f} s")166167168# THE MONTE CARLO VERB169170171def verb_montecarlo():172 t0 = time.time()173 rng = np.random.default_rng(7)174 cubes = prefractal()175 oracle = Oracle(cubes)176 print("MONTECARLO: T(delta) from the reduced distance of the lemmas, and a lemma-free bracket from the level-4 prefractal")177 n, chunk = 1_000_000, 100_000178 margin = SQRT2 / (6 * SIDE)179 low = {Fr(1, 8): 0, Fr(1, 12): 0}180 high = {Fr(1, 8): 0, Fr(1, 12): 0}181 for _ in range(n // chunk):182 d_f4, complete = oracle.dist(sample_plus(rng, chunk), k=200)183 assert complete.all()184 for d in low:185 high[d] += int((d_f4 <= float(d)).sum())186 low[d] += int((d_f4 <= float(d) - margin).sum())187 for d in low:188 print(f" lemma-free, {n} points against F_4: T({d}) in about [{low[d] / n * PLUS:.5f}, {high[d] / n * PLUS:.5f}], the count within delta an upper bound on T since F lies in F_4, the count within delta - sqrt(2)/486 a lower bound since every point of a kept cube is within sqrt(2)/486 of F, Monte Carlo 3 sigma {3 * np.sqrt(0.25 / n) * PLUS:.5f}")189 n, chunk = 40_000_000, 2_000_000190 hits = {Fr(1, 12): 0, Fr(1, 8): 0, Fr(1, 6): 0}191 for _ in range(n // chunk):192 d_red = reduced(sample_plus(rng, chunk))193 for d in hits:194 hits[d] += int((d_red <= float(d)).sum())195 for d in hits:196 value, err = three_sigma(hits[d], n, PLUS)197 print(f" reduced distance, {n} points: T({d}) = {value:.6f} +- {err:.6f} (3 sigma)")198 delta = Fr(1, 12)199 d = float(delta)200 a1 = st.lower(st.strip(delta)[0])201 n, chunk = 20_000_000, 5_000_000202 hits = 0203 for _ in range(n // chunk):204 u, v, z = rng.random(chunk) * d, rng.random(chunk) * d, rng.random(chunk) / 3205 du, dv = carpet_dist(v, z, 1 / 3), carpet_dist(u, z, 1 / 3)206 hits += int(((u * u + du * du <= d * d) & (v * v + dv * dv <= d * d)).sum())207 value, err = three_sigma(hits, n, d * d / 3)208 print(f" overlap of the two wall tubes, {n} points of the column: V2({delta}) = {value:.8f} +- {err:.1e}, against 2 A1 - delta^2/3 = {float(2 * a1 - mpf(d) ** 2 / 3):.8f} with A1 from the generator")209 print(f" wall {time.time() - t0:.1f} s")210211212# THE SEEDED VERB213214215def verb_seeded():216 t0 = time.time()217 rng = np.random.default_rng(20260921)218 print("SEEDED: Deep sampled on its confinement box [delta - w, delta]^2 x [0, 1/3], w = (1/81)/(4 delta), and V2 on the column, seed 20260921")219 for delta in (Fr(1, 6), Fr(1, 8)):220 d = float(delta)221 w = (1 / 81) / (4 * d)222 bound = float(st.deep_bound(delta))223 n, chunk = 10_000_000, 1_000_000224 deep_hits = 0225 v2_hits = 0226 for _ in range(n // chunk):227 u, v, z = d - w * rng.random(chunk), d - w * rng.random(chunk), rng.random(chunk) / 3228 du, dv = carpet_dist(v, z, 1 / 3), carpet_dist(u, z, 1 / 3)229 deep_hits += int(((u * u + du * du > d * d) & (v * v + dv * dv > d * d)).sum())230 u, v, z = d * rng.random(chunk), d * rng.random(chunk), rng.random(chunk) / 3231 du, dv = carpet_dist(v, z, 1 / 3), carpet_dist(u, z, 1 / 3)232 v2_hits += int(((u * u + du * du <= d * d) & (v * v + dv * dv <= d * d)).sum())233 deep, deep_err = three_sigma(deep_hits, n, w * w / 3)234 v2, v2_err = three_sigma(v2_hits, n, d * d / 3)235 a1 = st.lower(st.strip(delta)[0])236 ident = float(2 * a1 - mpf(d) ** 2 / 3) + deep237 print(f" delta = {delta}: Deep = {deep:.4e} +- {deep_err:.1e} (3 sigma, {deep_hits} hits of {n}), bound {bound:.4e}, ratio {bound / deep:.1f}")238 print(f" delta = {delta}: V2 = {v2:.7f} +- {v2_err:.1e} (3 sigma), against 2 A1 - delta^2/3 + Deep = {ident:.7f}, difference {v2 - ident:.1e}")239 print(f" wall {time.time() - t0:.1f} s")240241242# THE RECOMPUTE VERB, OWN CLOSED FORMS AT 60 DIGITS243244245def R(x):246 return mpf(x.numerator) / x.denominator if isinstance(x, Fr) else mpf(x)247248249def a0(tau, d):250 b = sqrt(d * d - tau * tau)251 return (tau * b + d * d * asin(tau / d)) / 2252253254def a1(tau, d):255 b = sqrt(d * d - tau * tau)256 return tau * tau * (d * d + d * b + b * b) / (3 * (d + b))257258259def a1_naive(tau, d):260 return (d**3 - (d * d - tau * tau) ** mpf(1.5)) / 3261262263def hole(s, d, prim=a1):264 tau = R(min(d, s / 2))265 s, d = R(s), R(d)266 return 4 * s * a0(tau, d) - 8 * prim(tau, d)267268269def cut_hole(s, c, d, prim=a1):270 if c <= s / 2:271 tau = R(min(c, d))272 s, c, d = R(s), R(c), R(d)273 return (s + 2 * c) * a0(tau, d) - 4 * prim(tau, d)274 t1, t2 = R(min(s - c, d)), R(min(s / 2, d))275 s, c, d = R(s), R(c), R(d)276 return (s + 2 * c) * a0(t1, d) - 4 * prim(t1, d) + 4 * s * (a0(t2, d) - a0(t1, d)) - 8 * (prim(t2, d) - prim(t1, d))277278279def column_holes(i, m):280 n = 1281 for _ in range(m - 1):282 n *= 2 if i % 3 == 1 else 3283 i //= 3284 return n285286287def columns_below(count, digits):288 if digits == 0:289 return 1 if count >= 1 else 0290 q, rem = divmod(count, 3 ** (digits - 1))291 if q >= 3:292 return 8**digits293 w = (3, 2, 3)294 return sum(w[:q]) * 8 ** (digits - 1) + w[q] * columns_below(rem, digits - 1)295296297def wall_half(d, levels, prim=a1):298 return sum(8 ** (m - 1) * hole(Fr(1, 3 ** (m + 1)), d, prim) for m in range(1, levels + 1)) / 2299300301def strip(d, levels, prim=a1):302 total = mpf(0)303 for m in range(1, levels + 1):304 s = Fr(1, 3 ** (m + 1))305 count = (d / s - 2) // 3 + 1 if d / s >= 2 else 0306 count = min(count, 3 ** (m - 1))307 if count > 0:308 total += columns_below(count, m - 1) * hole(s, d, prim)309 i = (d / s - 1) // 3310 if 0 <= i < 3 ** (m - 1) and (3 * i + 1) * s < d < (3 * i + 2) * s:311 total += column_holes(i, m) * cut_hole(s, d - (3 * i + 1) * s, d, prim)312 return total313314315def tube_open(d, levels):316 x = R(d)317 return (pi + 8) * x**2 - 8 * sqrt(2) * x**3 + 48 * (wall_half(d, levels) - strip(d, levels))318319320def tube_upper(d):321 x = R(d)322 total = mpf(0)323 for m in range(1, 200):324 s = Fr(1, 3 ** (m + 1))325 total += 8 ** (m - 1) * R(s) * x * R(min(s, 4 * d))326 return pi * x**2 + 24 * (total + x * R(Fr(8, 9) ** 199 / 8))327328329def hole_tables(levels):330 cells = 3**levels331 tables = {}332 for m in range(1, levels + 1):333 span = 3 ** (levels - m)334 table = np.zeros((3 ** (m - 1), cells))335 for i in range(3 ** (m - 1)):336 rows = [0]337 for k in range(m - 1):338 allowed = (0, 2) if (i // 3**k) % 3 == 1 else (0, 1, 2)339 rows = [r + a * 3**k for r in rows for a in allowed]340 for r in rows:341 table[i, (3 * r + 1) * span : (3 * r + 2) * span] = 1.0342 tables[m] = table343 return tables344345346def deep_bound(delta, tables, levels=7, tail_levels=80):347 delta = float(delta)348 cells = 3**levels349 boxes = 0.0350 for m in range(1, levels + 1):351 sm = 3.0 ** (-(m + 1))352 wm = sm * sm / (4 * delta)353 for n in range(1, levels + 1):354 sn = 3.0 ** (-(n + 1))355 wn = sn * sn / (4 * delta)356 i = np.arange(3 ** (m - 1))357 vlen = np.clip(np.minimum((3 * i + 2) * sm, delta) - np.maximum((3 * i + 1) * sm, delta - wn), 0, None)358 k = np.arange(3 ** (n - 1))359 ulen = np.clip(np.minimum((3 * k + 2) * sn, delta) - np.maximum((3 * k + 1) * sn, delta - wm), 0, None)360 boxes += float((vlen @ tables[m]) @ (ulen @ tables[n])) / cells / 3361 tail = 0.0362 for m in range(1, tail_levels + 1):363 for n in range(1, tail_levels + 1):364 if max(m, n) <= levels:365 continue366 wm = min(3.0 ** (-(2 * m + 2)) / (4 * delta), delta)367 wn = min(3.0 ** (-(2 * n + 2)) / (4 * delta), delta)368 count = min(np.ceil(wn * 3.0**m) + 2, np.ceil(wm * 3.0**n) + 2)369 tail += wm * wn * min(count / 9, 1 / 3)370 tail += 3.0 ** (-(2 * tail_levels + 4)) / (1536 * delta * delta)371 return boxes + tail372373374def verb_recompute():375 t0 = time.time()376 mp.dps = 60377 levels = 200378 print("RECOMPUTE: the tube and the bands from the closed forms of the paper alone, 60 decimal digits, 200 hole levels")379 for m in range(1, 7):380 digits = m - 1381 brute = [sum(1 for r in range(3**digits) if all(not ((i // 3**k) % 3 == 1 and (r // 3**k) % 3 == 1) for k in range(digits))) for i in range(3**digits)]382 assert brute == [column_holes(i, m) for i in range(3**digits)], m383 assert all(columns_below(count, digits) == sum(brute[:count]) for count in range(3**digits + 1)), m384 print(" column counts of the strip validated by brute force to level 6")385 for s, d in ((Fr(1, 27), Fr(1, 12)), (Fr(1, 9), Fr(1, 6))):386 q = quad(lambda t: 4 * (R(s) - 2 * t) * sqrt(R(d) ** 2 - t**2), [0, min(R(d), R(s) / 2)])387 print(f" J({s}, {d}) closed form minus quadrature: {mp.nstr(hole(s, d) - q, 3)}")388 mp.dps = 20389 for s, c, d in ((Fr(1, 9), Fr(1, 72), Fr(1, 8)), (Fr(1, 9), Fr(7, 108), Fr(1, 6)), (Fr(1, 81), Fr(1, 200), Fr(1, 12))):390 s_, c_, d_ = R(s), R(c), R(d)391392 def inner(v):393 def f(z):394 r = min(v, s_ - v, z, s_ - z)395 return sqrt(d_**2 - r**2) if r < d_ else mpf(0)396397 return quad(f, [0, min(v, s_ - v), s_ / 2, s_ - min(v, s_ - v), s_])398399 q = quad(inner, [0, min(c_, s_ / 2), c_] if c_ > s_ / 2 else [0, c_])400 print(f" cut hole (s, c, delta) = ({s}, {c}, {d}) closed form minus 2D quadrature: {mp.nstr(cut_hole(s, c, d) - q, 3)}")401 mp.dps = 60402 tail_v1 = lambda d: R(d) * R(Fr(8, 9) ** levels) / 18403 tail_a1 = lambda d: R(d) * R(Fr(8, 9) ** levels) / 9404 tables = hole_tables(7)405 for d in (Fr(1, 12), Fr(1, 8), Fr(1, 6)):406 v1, s1 = wall_half(d, levels), strip(d, levels)407 t = tube_open(d, levels)408 deep = deep_bound(d, tables)409 print(f" delta = {d}: V1 = {mp.nstr(v1, 15)}, A1 = {mp.nstr(s1, 15)}, Deep bound {deep:.4e}, T in [{mp.nstr(t - 48 * tail_a1(d) - 24 * mpf(deep), 12)}, {mp.nstr(t + 48 * tail_v1(d), 12)}]")410 weight = mpf(27) / 20411 dim = log(20) / log(3)412 for eps in (Fr(1, 12), Fr(1, 8), Fr(1, 6)):413 lo = hi = mpf(20) / 27414 for l in range(41):415 d = eps / 3**l416 t = tube_open(d, levels)417 deep = deep_bound(d, tables)418 lo += weight**l * (t - 48 * tail_a1(d) - 24 * mpf(deep))419 hi += weight**l * (t + 48 * tail_v1(d))420 series_tail = weight**41 * tube_upper(eps / 3**41) * mpf(20) / 9421 hi += series_tail422 scale = exp((dim - 3) * log(R(eps)))423 print(f" eps = {eps}: p in [{mp.nstr(lo * scale, 12)}, {mp.nstr(hi * scale, 12)}], series tail {mp.nstr(series_tail, 3)}")424 mp.dps = 200425 d = Fr(1, 12)426 v1, s1 = wall_half(d, levels, a1_naive), strip(d, levels, a1_naive)427 mp.dps = 60428 print(f" delta = 1/12 at 200 digits with a1 = (delta^3 - (delta^2 - t^2)^(3/2))/3: V1 = {mp.nstr(v1, 15)}, A1 = {mp.nstr(s1, 15)}")429 n = 3000430 g = (np.arange(n) + 0.5) / (3 * n)431 v, z = np.meshgrid(g, g, indexing="ij")432 f = np.sqrt(np.maximum(float(d) ** 2 - carpet_dist(v, z, 1 / 3) ** 2, 0))433 cell = (1 / (3 * n)) ** 2434 print(f" delta = 1/12 on a {n}^2 midpoint grid of the carpet distance: V1 = {f.sum() * cell / 2:.9f}, A1 = {f[v <= float(d)].sum() * cell:.9f}")435 print(f" wall {time.time() - t0:.1f} s")436437438def main():439 verbs = {440 "oracle": verb_oracle,441 "montecarlo": verb_montecarlo,442 "seeded": verb_seeded,443 "recompute": verb_recompute,444 }445 want = sys.argv[1:] or list(verbs)446 for v in want:447 verbs[v]()448449450if __name__ == "__main__":451 main()