sponge.rs
20.2 kB · rust · 546 lines
1use std::f64::consts::{PI, SQRT_2};2use std::sync::OnceLock;34/// The largest radius the closed tube formula reaches, `1/6`: past it the walls across an arm and the centre cube's edge cylinders start to meet.5pub const EDGE: f64 = 1.0 / 6.0;67/// The covering radius of the plus, `sqrt(2)/6`, the distance from the centre of the cube to the sponge: every radius from here on swallows the plus whole.8pub const COVER: f64 = SQRT_2 / 6.0;910const PLUS: f64 = 7.0 / 27.0;11const DEPTH: usize = 40;12const DIGITS: usize = 40;13const FLOOR: i64 = 32;14const REACH: i64 = 14;15const NEGLIGIBLE: f64 = 1e-17;16const NODES: usize = 16;17const WEIGHT: [f64; 3] = [3.0, 2.0, 3.0];18const RUNNING: [f64; 3] = [0.0, 3.0, 5.0];1920// DISTANCE2122fn split(x: f64) -> (usize, f64) {23 let digit = (3.0 * x).floor().clamp(0.0, 2.0);24 (digit as usize, (3.0 * x - digit).clamp(0.0, 1.0))25}2627fn carpet(mut across: f64, mut along: f64) -> f64 {28 let mut scale = 1.0 / 3.0;29 for _ in 0..DIGITS {30 let ((a, fa), (b, fb)) = (split(across), split(along));31 if a == 1 && b == 1 {32 return scale * fa.min(1.0 - fa).min(fb).min(1.0 - fb);33 }34 across = fa;35 along = fb;36 scale /= 3.0;37 }38 0.039}4041fn gap(v: f64) -> f64 {42 (v - 1.0 / 3.0).abs().min((v - 2.0 / 3.0).abs())43}4445fn plus(point: [f64; 3]) -> f64 {46 let middle = point.map(|v| (1.0 / 3.0..=2.0 / 3.0).contains(&v));47 let gaps = point.map(gap);48 match (0..3).find(|&k| !middle[k]) {49 None => {50 let mut sorted = gaps;51 sorted.sort_by(f64::total_cmp);52 sorted[0].hypot(sorted[1])53 }54 Some(axis) => {55 let v = point[axis];56 let along = 3.0 * if v < 1.0 / 3.0 { v } else { v - 2.0 / 3.0 };57 let (i, j) = ((axis + 1) % 3, (axis + 2) % 3);58 let wall_i = carpet(3.0 * point[j] - 1.0, along) / 3.0;59 let wall_j = carpet(3.0 * point[i] - 1.0, along) / 3.0;60 gaps[i].hypot(wall_i).min(gaps[j].hypot(wall_j))61 }62 }63}6465// HOLES6667fn sweep(t: f64, radius: f64) -> f64 {68 if t >= radius {69 return PI * radius * radius / 4.0;70 }71 t / 2.0 * (radius * radius - t * t).sqrt() + radius * radius / 2.0 * (t / radius).asin()72}7374fn moment(t: f64, radius: f64) -> f64 {75 let r2 = radius * radius;76 if t >= radius {77 return r2 * radius / 3.0;78 }79 let rest = r2 - t * t;80 t * t * (3.0 * r2 * r2 - 3.0 * r2 * t * t + t.powi(4))81 / (3.0 * (r2 * radius + rest * rest.sqrt()))82}8384fn piece(from: f64, to: f64, width: f64, slope: f64, radius: f64) -> f64 {85 if to <= from {86 return 0.0;87 }88 width * (sweep(to, radius) - sweep(from, radius))89 - slope * (moment(to, radius) - moment(from, radius))90}9192fn hole(side: f64, radius: f64) -> f64 {93 4.0 * piece(0.0, side / 2.0, side, 2.0, radius)94}9596fn partial(side: f64, cut: f64, radius: f64) -> f64 {97 let left = piece(0.0, (side / 2.0).min(cut), side, 2.0, radius);98 let right = piece((side - cut).max(0.0), side / 2.0, side, 2.0, radius);99 let bottom = if cut <= side / 2.0 {100 piece(0.0, cut, cut, 1.0, radius)101 } else {102 piece(0.0, side - cut, cut, 1.0, radius) + piece(side - cut, side / 2.0, side, 2.0, radius)103 };104 left + right + 2.0 * bottom105}106107fn places(mut count: f64, digits: usize) -> Vec<usize> {108 let mut out = Vec::with_capacity(digits);109 for _ in 0..digits {110 let rest = (count / 3.0).floor();111 out.push(((count - 3.0 * rest) as usize).min(2));112 count = rest;113 }114 out115}116117fn column(index: f64, digits: usize) -> f64 {118 places(index, digits).iter().map(|&k| WEIGHT[k]).product()119}120121fn below(count: f64, digits: usize) -> f64 {122 if count >= 3f64.powi(digits as i32) {123 return 8f64.powi(digits as i32);124 }125 let (mut total, mut prefix, mut eight) = (0.0, 1.0, 8f64.powi(digits as i32));126 for &k in places(count, digits).iter().rev() {127 eight /= 8.0;128 total += prefix * RUNNING[k] * eight;129 prefix *= WEIGHT[k];130 }131 total132}133134fn levels(radius: f64) -> usize {135 ((1.0 / radius).log(3.0).ceil() as i64 + REACH).max(FLOOR) as usize136}137138fn width(level: usize) -> f64 {139 3f64.powi(-(level as i32) - 1)140}141142fn wall(radius: f64, top: usize) -> f64 {143 let mut total = 0.0;144 for level in 1..=top {145 total += 8f64.powi(level as i32 - 1) * hole(width(level), radius);146 }147 total / 2.0 + radius * (8.0f64 / 9.0).powi(top as i32) / 18.0148}149150fn strip(radius: f64, top: usize) -> f64 {151 let (mut total, mut share) = (0.0, 0.0);152 for level in 1..=top {153 let side = width(level);154 let ratio = radius / side;155 let columns = 3f64.powi(level as i32 - 1);156 let mut here = 0.0;157 let full = ((ratio - 2.0) / 3.0).floor() + 1.0;158 if full > 0.0 {159 here += below(full.min(columns), level - 1) * hole(side, radius);160 }161 let cut = ((ratio - 1.0) / 3.0).floor();162 let start = (3.0 * cut + 1.0) * side;163 if cut >= 0.0 && cut < columns && start < radius && radius < start + side {164 here += column(cut, level - 1) * partial(side, radius - start, radius);165 }166 total += here;167 share = here / (8f64.powi(level as i32 - 1) * radius * side * side);168 }169 total + share * radius * (8.0f64 / 9.0).powi(top as i32) / 9.0170}171172fn inside(radius: f64) -> f64 {173 let top = levels(radius);174 (PI + 8.0) * radius * radius - 8.0 * SQRT_2 * radius.powi(3)175 + 48.0 * (wall(radius, top) - strip(radius, top))176}177178fn crown(radius: f64) -> f64 {179 let lap = (radius * radius - EDGE * EDGE).max(0.0).sqrt();180 PI * radius * radius181 - 8.0 * SQRT_2 * radius.powi(3)182 - 4.0 * radius * radius * (EDGE / radius).min(1.0).acos()183 + (2.0 / 3.0 + 16.0 * lap * lap) * lap184 + 2.0 / 9.0185}186187// CORNER188189fn legendre() -> &'static [(f64, f64)] {190 static RULE: OnceLock<Vec<(f64, f64)>> = OnceLock::new();191 RULE.get_or_init(|| {192 let n = NODES as f64;193 (1..=NODES)194 .map(|i| {195 let mut x = (PI * (i as f64 - 0.25) / (n + 0.5)).cos();196 let mut slope = 1.0;197 for _ in 0..64 {198 let (mut low, mut high) = (1.0, x);199 for k in 2..=NODES {200 let k = k as f64;201 (low, high) = (high, ((2.0 * k - 1.0) * x * high - (k - 1.0) * low) / k);202 }203 slope = n * (x * high - low) / (x * x - 1.0);204 let step = high / slope;205 x -= step;206 if step.abs() < 1e-17 {207 break;208 }209 }210 (x, 2.0 / ((1.0 - x * x) * slope * slope))211 })212 .collect()213 })214}215216fn pair(radius: f64, start: f64, side: f64, line: f64) -> f64 {217 let (half, reach) = (side / 2.0, line - start);218 let angle = |depth: f64| (depth / radius).clamp(0.0, 1.0).asin();219 let edge = |y: f64| ((radius - y) * (radius + y)).max(0.0).sqrt().atan2(y);220 let (low, high) = (edge(line), edge(start));221 let mut cuts = vec![low, high];222 for depth in [half, reach, side - reach] {223 cuts.push(angle(depth));224 }225 for shift in [start, start + side] {226 let spread = 2.0 * radius * radius - shift * shift;227 if spread >= 0.0 {228 cuts.push(angle(229 ((radius - shift) * (radius + shift) / (spread.sqrt() + shift)).abs(),230 ));231 }232 }233 cuts.retain(|phi| (low..=high).contains(phi));234 cuts.sort_by(f64::total_cmp);235 let mut total = 0.0;236 for window in cuts.windows(2) {237 let (middle, span) = ((window[0] + window[1]) / 2.0, (window[1] - window[0]) / 2.0);238 for &(x, weight) in legendre() {239 let phi = middle + span * x;240 let depth = radius * phi.sin();241 let across = radius - start - 2.0 * radius * (phi / 2.0).sin().powi(2);242 let width = reach.min(side - depth) - across.max(depth);243 if depth < half && width > 0.0 {244 total += weight * span * (half - depth) * width * depth;245 }246 }247 }248 4.0 * total249}250251// EXPORTS252253/// The Minkowski dimension of the Menger sponge, `log(20)/log(3)`, the similarity dimension of its 20 maps of ratio `1/3`.254pub fn dimension() -> f64 {255 20f64.ln() / 3f64.ln()256}257258/// The Euclidean distance from a point of the unit cube to the Menger sponge, exact to the last binary place.259///260/// The point descends the base-3 digits while its digit triple holds at most one `1`, each step261/// dividing the distance by 3, since inside a retained subcube the distance to the sponge is the262/// distance to that subcube's copy of it. Once it lands in the plus of seven removed cubes, the263/// distance is the one to the twelve edges of the centre cube or to the four carpets walling the264/// arm it sits in. Coordinates outside `[0, 1]` are clamped; the window centre `(1/2, 1/2, 0)`265/// reads `1/6` and the centre of the cube `sqrt(2)/6`.266pub fn distance(point: [f64; 3]) -> f64 {267 let mut local = point.map(|v| v.clamp(0.0, 1.0));268 let mut scale = 1.0;269 for _ in 0..DIGITS {270 let digits = local.map(split);271 if digits.iter().filter(|(d, _)| *d == 1).count() >= 2 {272 return scale * plus(local);273 }274 local = digits.map(|(_, rest)| rest);275 scale /= 3.0;276 }277 0.0278}279280/// The volume `T(radius)` of the points of the plus of seven removed level-1 cubes within `radius` of the sponge, or `None` on `(1/6, sqrt(2)/6)`, where no closed form is known.281///282/// On `(0, 1/6]` it is `(pi + 8) r^2 - 8 sqrt(2) r^3 + 48 (V1 - A1)`, the centre cube's edge283/// cylinders plus the arcsine hole sums of the 24 wall carpets, `V1` over half a wall and `A1`284/// over the strip along a shared edge, both summed level by level and closed by their geometric285/// tails. The one term left out is `24 Deep`, `Deep` the corner volume beyond both wall tubes,286/// at most `24 x 3.842e-5 < 9.23e-4` at `1/6` and `5.81e-9` at `1/12` by the paper's bound on287/// `Deep`, so this is the upper edge of the certified enclosure:288/// `T(1/6) = (pi + 8)/36 - sqrt(2)/27`. From `sqrt(2)/6` on the plus is swallowed, `7/27`.289pub fn tube(radius: f64) -> Option<f64> {290 if radius.is_nan() || (radius > EDGE && radius < COVER) {291 return None;292 }293 if radius <= 0.0 {294 return Some(0.0);295 }296 if radius >= COVER {297 return Some(PLUS);298 }299 Some(inside(radius))300}301302/// The volume `Deep(radius)` of the points of one arm's quarter beyond the tubes of both walls it touches, at every radius, to double precision.303///304/// With `line = min(radius, 1/6)`, the quarter is `[0, line]^2` across the arm times its length,305/// and every point of `Deep` sits, through both of its wall coordinates, in one and the same hole306/// of the wall carpet that the line `across = line` cuts: a crossing hole. The volume is the sum307/// over crossing holes of one arcsine integral, `2^(m-1)` holes of side `3^(-m-1)` at every level308/// `m` for `line = 1/6`, none at all for `1/12`, whose triple `1/4 = 0.0202...` has no digit `1`.309/// It vanishes from `sqrt(10)/18` on, where the arms are swallowed. Each integral is elementary, but its310/// expanded arcsine form cancels at deep levels, so it is summed by 16-point Gauss-Legendre on each smooth311/// piece in the angle `asin(R/radius)`, `R` the depth that a wall point needs.312pub fn deep(radius: f64) -> f64 {313 if !(radius > 0.0 && radius.is_finite()) {314 return 0.0;315 }316 let line = radius.min(EDGE);317 let square = radius * radius;318 let (mut local, mut start, mut count, mut total) = (3.0 * line, 0.0, 1.0, 0.0);319 for level in 1..=levels(radius) {320 let side = width(level);321 let digit = (3.0 * local).floor().clamp(0.0, 2.0);322 local = 3.0 * local - digit;323 if digit == 1.0 {324 if count * 2.0 * (side / 2.0).powi(5) < NEGLIGIBLE * square * square {325 break;326 }327 total += count * pair(radius, start + side, side, line);328 }329 start += digit * side;330 count *= if digit == 1.0 { 2.0 } else { 3.0 };331 }332 total333}334335/// The volume `T(radius)` of the points of the plus within `radius` of the sponge at every radius, `Deep` included, to double precision.336///337/// On `(0, 1/6]` it is the closed form of `tube` less `24 Deep`. On `[1/6, sqrt(2)/6]` the arms338/// hold `2/9 - 24 Deep` and the centre cube `pi r^2 - 8 sqrt(2) r^3 - 4 r^2 acos(1/(6r)) +339/// (2/3 + 16 c^2) c` with `c = sqrt(r^2 - 1/36)`, the edge cylinders overlapping along the faces;340/// from `sqrt(2)/6` on the plus is swallowed, `7/27`.341pub fn exact(radius: f64) -> f64 {342 if radius <= 0.0 {343 return 0.0;344 }345 if radius >= COVER {346 return PLUS;347 }348 let body = if radius <= EDGE {349 inside(radius)350 } else {351 crown(radius)352 };353 body - 24.0 * deep(radius)354}355356/// The volume of the points of the unit cube within `radius` of the sponge, `sum_k (20/27)^k T(3^k radius)`, or `None` when some `3^k radius` falls where `T` has no closed form.357///358/// The cube is the plus and 20 retained subcubes, each holding a copy of the whole picture at a359/// third of the scale, so the volume is `T(radius)` plus `20/27` of the volume at three times the360/// radius, and the sum stops at the first radius past `sqrt(2)/6`, where the volume is the cube.361pub fn volume(radius: f64) -> Option<f64> {362 if radius.is_nan() {363 return None;364 }365 if radius <= 0.0 {366 return Some(0.0);367 }368 let (mut total, mut weight, mut reach) = (0.0, 1.0, radius);369 while reach < COVER {370 total += weight * tube(reach)?;371 weight *= 20.0 / 27.0;372 reach *= 3.0;373 }374 Some(total + weight)375}376377/// The Minkowski reading `radius^(D-3)` times the volume inside the cube, the number whose limit as the radius shrinks would be the sponge's Minkowski content.378pub fn reading(radius: f64) -> Option<f64> {379 Some(radius.powf(dimension() - 3.0) * volume(radius)?)380}381382/// The periodic function `p` of Kombrink, Pearse and Winter at `radius`: the reading's limit profile, unchanged when the radius is multiplied by 3, or `None` on the phases `(1/6, sqrt(2)/6]` where `T` has no closed form.383///384/// The radius is first carried by powers of 3 into `(sqrt(2)/18, sqrt(2)/6]`; on385/// `(sqrt(2)/18, 1/6]` it is `r^(D-3) (20/27 + sum_(l = 0..40) (27/20)^l T(r/3^l))`, and the386/// sponge is Minkowski measurable exactly when this is constant. It is not: the certificate puts387/// `p(1/12)` in `[2.122718, 2.122723]` and `p(1/6)` in `[2.135019, 2.136794]`.388pub fn profile(radius: f64) -> Option<f64> {389 if !(radius > 0.0 && radius.is_finite()) {390 return None;391 }392 let mut phase = radius;393 while phase > COVER {394 phase /= 3.0;395 }396 while phase <= COVER / 3.0 {397 phase *= 3.0;398 }399 if phase > EDGE {400 return None;401 }402 let (mut total, mut weight, mut reach) = (20.0 / 27.0, 1.0, phase);403 for _ in 0..=DEPTH {404 total += weight * tube(reach)?;405 weight *= 27.0 / 20.0;406 reach /= 3.0;407 }408 Some(phase.powf(dimension() - 3.0) * total)409}410411#[cfg(test)]412mod tests {413 use super::*;414415 #[test]416 fn distance_reads_the_named_points() {417 let near = |a: f64, b: f64| (a - b).abs() < 1e-15;418 assert!(near(distance([0.5, 0.5, 0.5]), SQRT_2 / 6.0));419 assert!(near(distance([1.0 / 6.0, 1.0 / 6.0, 0.5]), SQRT_2 / 18.0));420 assert!(near(distance([0.5, 0.5, 1.0]), 1.0 / 6.0));421 assert!(near(distance([2.0 / 3.0, 0.5, 5.0 / 6.0]), 1.0 / 18.0));422 assert_eq!(distance([0.0, 0.25, 1.0]), 0.0);423 }424425 #[test]426 fn tube_meets_the_identity_at_one_sixth() {427 let closed = (PI + 8.0) / 36.0 - SQRT_2 / 27.0;428 assert!((tube(EDGE).unwrap() - closed).abs() < 1e-13);429 assert_eq!(tube(0.2), None);430 assert_eq!(tube(COVER), Some(PLUS));431 }432433 #[test]434 fn tube_lies_in_its_enclosures() {435 let twelfth = tube(1.0 / 12.0).unwrap();436 let eighth = tube(1.0 / 8.0).unwrap();437 assert!((0.180947086..=0.180947093).contains(&twelfth));438 assert!((0.234186414..=0.234701259).contains(&eighth));439 }440441 #[test]442 fn tube_is_bracketed_by_the_distance_raster() {443 let n = 36;444 let cell = 1.0 / (3.0 * n as f64);445 let half = cell * 3f64.sqrt() / 2.0;446 let mut seen = Vec::with_capacity(2 * n * n * n);447 for (corner, copies) in [([1.0, 1.0, 1.0], 1.0), ([1.0, 1.0, 0.0], 6.0)] {448 for a in 0..n {449 for b in 0..n {450 for c in 0..n {451 let at = [a, b, c].map(|v| (v as f64 + 0.5) * cell);452 let point = [0, 1, 2].map(|k| corner[k] / 3.0 + at[k]);453 seen.push((distance(point), copies));454 }455 }456 }457 }458 let volume = cell.powi(3);459 for (radius, closed) in [(1.0 / 12.0, tube(1.0 / 12.0).unwrap()), (0.2, exact(0.2))] {460 let count = |edge: f64| {461 seen.iter()462 .filter(|(d, _)| *d <= edge)463 .map(|(_, c)| c)464 .sum::<f64>()465 };466 assert!(467 count(radius - half) * volume <= closed && closed <= count(radius + half) * volume468 );469 }470 }471472 #[test]473 fn deep_lives_in_the_crossing_holes() {474 assert_eq!(deep(1.0 / 12.0), 0.0);475 assert_eq!(deep(0.18), 0.0);476 assert!((deep(EDGE) - 1.8882e-6).abs() < 1.4e-8);477 assert!((deep(1.0 / 8.0) - 4.7432e-8).abs() < 2.9e-9);478 let terms: Vec<f64> = (1..=10)479 .map(|m| {480 let side = width(m);481 2f64.powi(m as i32 - 1) * pair(EDGE, EDGE - side / 2.0, side, EDGE)482 })483 .collect();484 assert!(terms485 .windows(2)486 .skip(1)487 .all(|p| (p[1] / p[0] * 243.0 / 2.0 - 1.0).abs() < 0.02));488 let radius: f64 = 0.17;489 let low = (radius * radius - 1.0 / 324.0).sqrt();490 let (n, m) = (48, 240);491 let (step, pace) = ((EDGE - low) / n as f64, (1.0 / 9.0) / m as f64);492 let half = (2.0 * step * step + pace * pace).sqrt() / 2.0;493 let (mut inner, mut outer) = (0.0, 0.0);494 for a in 0..n {495 for b in 0..n {496 for c in 0..m {497 let u = 1.0 / 3.0 + low + (a as f64 + 0.5) * step;498 let v = 1.0 / 3.0 + low + (b as f64 + 0.5) * step;499 let d = distance([1.0 / 9.0 + (c as f64 + 0.5) * pace, u, v]);500 inner += f64::from(u8::from(d > radius + half));501 outer += f64::from(u8::from(d > radius - half));502 }503 }504 }505 let volume = step * step * pace;506 let closed = deep(radius);507 assert!(inner * volume <= closed && closed <= outer * volume);508 }509510 #[test]511 fn exact_joins_at_one_sixth_and_lies_in_the_enclosures() {512 assert!((crown(EDGE) - inside(EDGE)).abs() < 1e-13);513 assert!((crown(COVER) - PLUS).abs() < 1e-15);514 assert!((0.180947086..=0.180947093).contains(&exact(1.0 / 12.0)));515 assert!((0.234186414..=0.234701259).contains(&exact(1.0 / 8.0)));516 assert!((0.256188319..=0.257110405).contains(&exact(EDGE)));517 let walk: Vec<f64> = (1..=400).map(|k| exact(k as f64 * COVER / 400.0)).collect();518 assert!(walk.windows(2).all(|pair| pair[0] < pair[1]));519 assert_eq!(exact(COVER), PLUS);520 }521522 #[test]523 fn profile_lies_in_its_bands_and_swings() {524 let twelfth = profile(1.0 / 12.0).unwrap();525 let eighth = profile(1.0 / 8.0).unwrap();526 let sixth = profile(1.0 / 6.0).unwrap();527 assert!((2.122718..=2.122723).contains(&twelfth));528 assert!((2.134668..=2.135742).contains(&eighth));529 assert!((2.135019..=2.136794).contains(&sixth));530 assert!(sixth - twelfth >= 0.012296);531 assert!((profile(1.0 / 36.0).unwrap() - twelfth).abs() < 1e-12);532 assert_eq!(profile(0.2), None);533 }534535 #[test]536 fn reading_climbs_onto_the_profile() {537 let radius = 1.0 / 12.0;538 let gaps: Vec<f64> = (0..10)539 .map(|k| radius / 3f64.powi(k))540 .map(|r| profile(r).unwrap() - reading(r).unwrap())541 .collect();542 assert!(gaps[0] > 0.1 && gaps[9] > 0.0 && gaps[9] < gaps[0] / 20.0);543 assert!(gaps.windows(2).all(|pair| pair[1] < pair[0]));544 assert_eq!(volume(COVER), Some(1.0));545 }546}