wiki-minkowski-content.rs

1use figures::{ink, plot, save, Board, Frame};2use mrlyrs::core::error::Result;34const ROWS: usize = 7;5const WIDEST: f64 = 0.1;6const THREAD: usize = 5;7const FROM: f64 = 2.0;8const TO: f64 = 6.0;9const SAMPLES: usize = 4000;10const LOW: f64 = 2.494975716;1112fn main() -> Result<()> {13    let dimension = 2f64.ln() / 3f64.ln();14    let mut board = Board::square();15    let frame = board.frame(0.08);16    let upper = Frame::new(frame.x, frame.y, frame.w, frame.h * 0.58);17    let lower = Frame::new(frame.x, frame.y + frame.h * 0.66, frame.w, frame.h * 0.34);1819    let span = 1.0 + 2.0 * WIDEST;20    let px = |u: f64| upper.x + upper.w * (u + WIDEST) / span;21    let pitch = upper.h / ROWS as f64;22    let thick = pitch * 0.62;23    let fine = pitch * 0.16;24    let thread = intervals(THREAD);25    for row in 0..ROWS {26        let eps = WIDEST * 3f64.powf(-(row as f64) / 2.0);27        let top = upper.y + row as f64 * pitch + (pitch - thick) / 2.0;28        let tube = tube(eps);29        let length: f64 = tube.iter().map(|(a, b)| b - a).sum();30        assert!((length - volume(eps)).abs() < 1e-12);31        for (a, b) in tube {32            board.rect(px(a), top, px(b) - px(a), thick, ink::blue());33        }34        let mid = top + (thick - fine) / 2.0;35        for &(a, b) in &thread {36            board.rect(px(a), mid, px(b) - px(a), fine, ink::yellow());37        }38    }3940    let xs: Vec<f64> = (0..=SAMPLES)41        .map(|i| FROM + (TO - FROM) * i as f64 / SAMPLES as f64)42        .collect();43    let ys: Vec<f64> = xs44        .iter()45        .map(|&x| reading(3f64.powf(-x), dimension))46        .collect();47    let high = 2f64.powf(2.0 - dimension);48    let lo = ys.iter().copied().fold(f64::MAX, f64::min);49    let hi = ys.iter().copied().fold(f64::MIN, f64::max);50    assert!((lo - LOW).abs() < 1e-5);51    assert!(hi <= high + 1e-9);52    assert!((reading(1.0 / 18.0, dimension) - high).abs() < 1e-9);53    for k in 2..=6 {54        assert!((reading(3f64.powi(-k), dimension) - 2.5).abs() < 1e-9);55    }56    for &x in xs.iter().step_by(97) {57        let eps = 3f64.powf(-x);58        assert!((reading(eps / 3.0, dimension) - reading(eps, dimension)).abs() < 1e-9);59    }60    let pad = (high - LOW) * 0.12;61    let (bottom, ceiling) = (LOW - pad, high + pad);62    let hair = (lower.h / 256.0).max(1.0);63    for k in (FROM as usize + 1)..(TO as usize) {64        let x = lower.x + lower.w * (k as f64 - FROM) / (TO - FROM);65        board.segment((x, lower.y), (x, lower.y + lower.h), hair, ink::dim());66    }67    let pts: Vec<(f64, f64)> = xs68        .iter()69        .zip(&ys)70        .map(|(&x, &y)| {71            (72                lower.x + lower.w * (x - FROM) / (TO - FROM),73                lower.y + lower.h * (1.0 - (y - bottom) / (ceiling - bottom)),74            )75        })76        .collect();77    board.polyline(&pts, (lower.h / 64.0).max(2.0), ink::yellow());78    plot::axis(&mut board, lower, ink::line());79    save("wiki-minkowski-content", &board)?;80    Ok(())81}8283fn intervals(level: usize) -> Vec<(f64, f64)> {84    let side = 3f64.powi(-(level as i32));85    (0..1usize << level)86        .map(|i| {87            let a: f64 = (0..level)88                .filter(|&l| (i >> (level - 1 - l)) & 1 == 1)89                .map(|l| 2.0 * 3f64.powi(-(l as i32 + 1)))90                .sum();91            (a, a + side)92        })93        .collect()94}9596fn tube(eps: f64) -> Vec<(f64, f64)> {97    let level = (0..40)98        .find(|&k| 3f64.powi(-(k as i32 + 1)) <= 2.0 * eps)99        .unwrap_or(40);100    let mut out: Vec<(f64, f64)> = Vec::new();101    for (a, b) in intervals(level) {102        let (a, b) = (a - eps, b + eps);103        match out.last_mut() {104            Some(last) if a <= last.1 => last.1 = last.1.max(b),105            _ => out.push((a, b)),106        }107    }108    out109}110111fn volume(eps: f64) -> f64 {112    2.0 * eps113        + (1..80)114            .map(|k| 2f64.powi(k - 1) * 3f64.powi(-k).min(2.0 * eps))115            .sum::<f64>()116}117118fn reading(eps: f64, dimension: f64) -> f64 {119    eps.powf(dimension - 1.0) * volume(eps)120}