wiki-rajchman-measure.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board, Frame};3use std::f64::consts::PI;45const TOP: usize = 243;6const DEPTH: u32 = 40;7const WIDTH: f64 = 0.1;8const STEPS: usize = 4096;910fn main() -> Result<()> {11    let mut board = Board::square();12    let frame = board.frame(0.08);13    let gap = frame.h * 0.08;14    let half = (frame.h - gap) / 2.0;15    let upper = Frame::new(frame.x, frame.y, frame.w, half);16    let lower = Frame::new(frame.x, frame.y + half + gap, frame.w, half);1718    let cantor: Vec<f64> = (1..=TOP).map(|n| cantor(n as f64)).collect();19    let bump: Vec<f64> = (1..=TOP).map(|n| bump(n as f64)).collect();2021    let powers: Vec<usize> = (0..)22        .map(|k| 3usize.pow(k))23        .take_while(|&p| p <= TOP)24        .collect();25    let ceiling = cantor[0];26    for &p in &powers {27        assert!((cantor[p - 1] - ceiling).abs() < 1e-9);28    }29    assert_eq!(powers.len(), 6);30    assert!(cantor.iter().all(|&v| v <= ceiling + 1e-9));31    assert!(bump[0] > 0.9);32    assert!(bump[40..].iter().all(|&v| v < 0.01));3334    comb(&mut board, upper, &cantor, ceiling, &powers);35    comb(&mut board, lower, &bump, ceiling.max(bump[0]), &[]);36    plot::axis(&mut board, upper, ink::line());37    plot::axis(&mut board, lower, ink::line());38    save("wiki-rajchman-measure", &board)?;39    Ok(())40}4142fn comb(board: &mut Board, frame: Frame, values: &[f64], peak: f64, lit: &[usize]) {43    let slot = frame.w / values.len() as f64;44    let pad = slot * 0.2;45    for (i, &v) in values.iter().enumerate() {46        let h = frame.h * v / peak;47        let color = if lit.contains(&(i + 1)) {48            ink::yellow()49        } else {50            ink::blue()51        };52        board.rect(53            frame.x + i as f64 * slot + pad,54            frame.y + frame.h - h,55            slot - 2.0 * pad,56            h,57            color,58        );59    }60}6162fn cantor(t: f64) -> f64 {63    (1..=DEPTH)64        .map(|k| (2.0 * PI * t / 3f64.powi(k as i32)).cos().abs())65        .product()66}6768fn bump(n: f64) -> f64 {69    let mut re = 0.0;70    let mut im = 0.0;71    let mut mass = 0.0;72    for i in 0..STEPS {73        let u = (i as f64 + 0.5) / STEPS as f64;74        let w = (-1.0 / (u * (1.0 - u))).exp();75        let x = 0.5 + WIDTH * (u - 0.5);76        re += w * (2.0 * PI * n * x).cos();77        im += w * (2.0 * PI * n * x).sin();78        mass += w;79    }80    re.hypot(im) / mass81}