paper-sponge-measurability.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board, Frame, Grid};3use mrlymath::three::{carpet, slice};4use std::f64::consts::PI;56const LEVEL: usize = 4;7const SIDE: usize = 81;8const EPS: f64 = 1.0 / 36.0;9const LEVELS: usize = 200;10const DEPTH: usize = 40;11const DIGITS: usize = 40;12const CELLW: [f64; 3] = [3.0, 2.0, 3.0];13const SUMW: [f64; 3] = [0.0, 3.0, 5.0];14const BANDS: [(f64, f64, f64); 3] = [15    (1.0 / 12.0, 2.122718, 2.122723),16    (1.0 / 8.0, 2.134668, 2.135742),17    (1.0 / 6.0, 2.135019, 2.136794),18];19const TUBES: [(f64, f64, f64); 2] = [20    (1.0 / 12.0, 0.180947086, 0.180947093),21    (1.0 / 8.0, 0.234186414, 0.234701259),22];2324// THE HOLE INTEGRALS2526fn a0(t: f64, d: f64) -> f64 {27    if t >= d {28        return PI * d * d / 4.0;29    }30    t / 2.0 * (d * d - t * t).sqrt() + d * d / 2.0 * (t / d).asin()31}3233fn a1(t: f64, d: f64) -> f64 {34    if t >= d {35        return d * d * d / 3.0;36    }37    let y = d * d - t * t;38    t * t * (3.0 * d.powi(4) - 3.0 * d * d * t * t + t.powi(4)) / (3.0 * (d.powi(3) + y * y.sqrt()))39}4041fn seg(t0: f64, t1: f64, alpha: f64, beta: f64, d: f64) -> f64 {42    if t1 <= t0 {43        return 0.0;44    }45    alpha * (a0(t1, d) - a0(t0, d)) - beta * (a1(t1, d) - a1(t0, d))46}4748fn hole(s: f64, d: f64) -> f64 {49    4.0 * seg(0.0, s / 2.0, s, 2.0, d)50}5152fn partial(s: f64, c: f64, d: f64) -> f64 {53    let left = seg(0.0, (s / 2.0).min(c), s, 2.0, d);54    let right = seg((s - c).max(0.0), s / 2.0, s, 2.0, d);55    let bottom = if c <= s / 2.0 {56        seg(0.0, c, c, 1.0, d)57    } else {58        seg(0.0, s - c, c, 1.0, d) + seg(s - c, s / 2.0, s, 2.0, d)59    };60    left + right + 2.0 * bottom61}6263fn digit(q: f64) -> usize {64    ((q - 3.0 * (q / 3.0).floor()) as usize).min(2)65}6667fn weight(mut i: f64, digits: usize) -> f64 {68    let mut w = 1.0;69    for _ in 0..digits {70        w *= CELLW[digit(i)];71        i = (i / 3.0).floor();72    }73    w74}7576fn below(count: f64, digits: usize) -> f64 {77    if count >= 3f64.powi(digits as i32) {78        return 8f64.powi(digits as i32);79    }80    let (mut total, mut prefix) = (0.0, 1.0);81    for d in (0..digits).rev() {82        let k = digit((count / 3f64.powi(d as i32)).floor());83        total += prefix * SUMW[k] * 8f64.powi(d as i32);84        prefix *= CELLW[k];85    }86    total87}8889fn wall_half(delta: f64) -> f64 {90    let mut total = 0.0;91    for m in 1..=LEVELS {92        total += 8f64.powi(m as i32 - 1) * hole(3f64.powi(-(m as i32) - 1), delta);93    }94    total / 2.095}9697fn strip(delta: f64) -> f64 {98    let mut total = 0.0;99    for m in 1..=LEVELS {100        let s = 3f64.powi(-(m as i32) - 1);101        let ratio = delta / s;102        let columns = 3f64.powi(m as i32 - 1);103        let full = ((ratio - 2.0) / 3.0).floor() + 1.0;104        if full > 0.0 {105            total += below(full.min(columns), m - 1) * hole(s, delta);106        }107        let cut = ((ratio - 1.0) / 3.0).floor();108        if cut >= 0.0109            && cut < columns110            && (3.0 * cut + 1.0) * s < delta111            && delta < (3.0 * cut + 2.0) * s112        {113            total += weight(cut, m - 1) * partial(s, delta - (3.0 * cut + 1.0) * s, delta);114        }115    }116    total117}118119fn tube(delta: f64) -> f64 {120    (PI + 8.0) * delta * delta - 8.0 * 2f64.sqrt() * delta.powi(3)121        + 48.0 * (wall_half(delta) - strip(delta))122}123124fn periodic(eps: f64) -> f64 {125    let dim = 20f64.ln() / 3f64.ln();126    let mut total = 20.0 / 27.0;127    for l in 0..=DEPTH {128        total += (27f64 / 20.0).powi(l as i32) * tube(eps / 3f64.powi(l as i32));129    }130    eps.powf(dim - 3.0) * total131}132133// THE DISTANCE IN THE MIDPLANE134135fn split(x: f64) -> (usize, f64) {136    let d = (3.0 * x).floor().min(2.0);137    (d as usize, 3.0 * x - d)138}139140fn carpet_dist(mut a: f64, mut b: f64) -> f64 {141    let mut s = 1.0 / 3.0;142    for _ in 0..DIGITS {143        let ((da, fa), (db, fb)) = (split(a), split(b));144        if da == 1 && db == 1 {145            return s * fa.min(1.0 - fa).min(fb).min(1.0 - fb);146        }147        a = fa;148        b = fb;149        s /= 3.0;150    }151    0.0152}153154fn wall(u: f64, along: f64, across: f64) -> f64 {155    let d = carpet_dist(3.0 * along, 3.0 * across) / 3.0;156    (u * u + d * d).sqrt()157}158159fn arm_dist(along: f64, u: f64) -> f64 {160    let side = wall(u, along, 1.0 / 6.0).min(wall(1.0 / 3.0 - u, along, 1.0 / 6.0));161    let far = wall(1.0 / 6.0, along, u);162    side.min(far)163}164165fn plus_dist(x: f64, y: f64) -> f64 {166    let (third, two) = (1.0 / 3.0, 2.0 / 3.0);167    let mid = |v: f64| (third..=two).contains(&v);168    if mid(x) && mid(y) {169        let mut best = f64::MAX;170        for cx in [third, two] {171            for cy in [third, two] {172                best = best.min(((x - cx).powi(2) + (y - cy).powi(2)).sqrt());173            }174        }175        let edge = (x - third)176            .abs()177            .min((x - two).abs())178            .min((y - third).abs())179            .min((y - two).abs());180        return best.min((1.0 / 36.0 + edge * edge).sqrt());181    }182    if mid(y) {183        let along = if x < third { x } else { 1.0 - x };184        return arm_dist(along, y - third);185    }186    let along = if y < third { y } else { 1.0 - y };187    arm_dist(along, x - third)188}189190fn dist(mut x: f64, mut y: f64) -> f64 {191    let mut scale = 1.0;192    for _ in 0..DIGITS {193        let ((dx, fx), (dy, fy)) = (split(x), split(y));194        if dx == 1 || dy == 1 {195            return scale * plus_dist(x, y);196        }197        x = fx;198        y = fy;199        scale /= 3.0;200    }201    0.0202}203204// THE MARKS205206fn bands(board: &mut Board, frame: Frame) {207    plot::axis(board, frame, ink::line());208    let area = frame.inset(24.0);209    let (lo, hi) = (2.120, 2.140);210    let y_of = |v: f64| area.y + area.h * (hi - v) / (hi - lo);211    let width = area.w * 0.3;212    let (a, b) = (BANDS[0], BANDS[2]);213    for (k, band) in [a, b].iter().enumerate() {214        let x = area.x + area.w * (0.25 + 0.5 * k as f64) - width / 2.0;215        let (top, foot) = (y_of(band.2), y_of(band.1));216        let h = (foot - top).max(4.0);217        board.rect(x, top, width, h, ink::orange());218    }219    let (cx, gap_lo, gap_hi) = (area.x + area.w / 2.0, y_of(a.2) - 2.0, y_of(b.1) + 2.0);220    board.segment((cx, gap_lo), (cx, gap_hi), 2.0, ink::dim());221    board.segment((cx - 10.0, gap_lo), (cx + 10.0, gap_lo), 2.0, ink::dim());222    board.segment((cx - 10.0, gap_hi), (cx + 10.0, gap_hi), 2.0, ink::dim());223}224225fn main() -> Result<()> {226    let sponge = carpet(3, LEVEL)?;227    assert_eq!(sponge.types().sum(), 160000);228    let dust = slice(&sponge, 2, (SIDE - 1) / 2)?;229    let cells = dust.types().bytes().to_vec();230    assert_eq!(cells.len(), SIDE * SIDE);231    assert_eq!(cells.iter().filter(|&&b| b != 0).count(), 256);232233    let ident = (PI + 8.0) / 36.0 - 2f64.sqrt() / 27.0;234    assert!((tube(1.0 / 6.0) - ident).abs() < 1e-9);235    for (delta, lo, hi) in TUBES {236        let t = tube(delta);237        assert!(lo <= t && t <= hi);238    }239    let mut values = [0.0; 3];240    for (k, (eps, lo, hi)) in BANDS.iter().enumerate() {241        values[k] = periodic(*eps);242        assert!(*lo <= values[k] && values[k] <= *hi);243    }244    assert!(values[2] - values[0] >= 0.012296);245    assert!((dist(0.5, 0.5) - 2f64.sqrt() / 6.0).abs() < 1e-15);246    assert!((dist(1.0 / 6.0, 1.0 / 6.0) - 2f64.sqrt() / 18.0).abs() < 1e-15);247248    let mut board = Board::square();249    let margin = (board.width as f64 * 0.08).round();250    let plate = 600.0;251    let sheet = Frame::new(margin, margin, plate, plate);252    board.rect(sheet.x, sheet.y, sheet.w, sheet.h, ink::panel());253    let n = plate as usize;254    let blue = ink::blue();255    for j in 0..n {256        for i in 0..n {257            let (x, y) = ((i as f64 + 0.5) / plate, (j as f64 + 0.5) / plate);258            let cover = (0.5 + (EPS - dist(x, y)) * plate).clamp(0.0, 1.0);259            if cover > 0.0 {260                board.blend(sheet.x as usize + i, sheet.y as usize + j, blue, cover);261            }262        }263    }264    let lattice = Grid::new(sheet, SIDE, SIDE, 0.0);265    for row in 0..SIDE {266        for col in 0..SIDE {267            if cells[row * SIDE + col] != 0 {268                lattice.fill(&mut board, col, row, ink::fg());269            }270        }271    }272    let panel = Frame::new(273        board.width as f64 - margin - 240.0,274        board.height as f64 - margin - 240.0,275        240.0,276        240.0,277    );278    bands(&mut board, panel);279    save("paper-sponge-measurability", &board)?;280    Ok(())281}