wiki-divisor-function.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board};34const X: usize = 36;5const ROOT: usize = 6;6const WINDOW: f64 = 38.0;7const SAMPLES: usize = 800;89fn main() -> Result<()> {10    let mut board = Board::square();11    let frame = board.frame(0.08);12    let unit = frame.w / WINDOW;13    let at = |a: f64, b: f64| (frame.x + a * unit, frame.y + frame.h - b * unit);1415    let divisors: Vec<usize> = (0..=X)16        .map(|n| (1..=n).filter(|d| n % d == 0).count())17        .collect();18    let total: usize = divisors[1..].iter().sum();19    let folded = 2 * (1..=ROOT).map(|a| X / a).sum::<usize>() - ROOT * ROOT;20    assert_eq!(ROOT * ROOT, X);21    assert_eq!(total, 140);22    assert_eq!(folded, total);23    assert_eq!(divisors[36], 9);24    assert_eq!(divisors[31], 2);2526    let mut inner = Vec::new();27    let mut arms = Vec::new();28    for a in 1..=X {29        for b in 1..=X / a {30            let p = at(a as f64, b as f64);31            if a <= ROOT && b <= ROOT {32                inner.push(p);33            } else {34                arms.push(p);35            }36        }37    }38    assert_eq!(inner.len(), ROOT * ROOT);39    assert_eq!(arms.len(), total - ROOT * ROOT);4041    let curve: Vec<(f64, f64)> = (0..=SAMPLES)42        .map(|i| {43            let a = (X as f64).powf(i as f64 / SAMPLES as f64);44            at(a, X as f64 / a)45        })46        .collect();47    let hair = (frame.w / 512.0).max(1.0);48    board.segment(at(0.0, 0.0), at(WINDOW, 0.0), hair, ink::line());49    board.segment(at(0.0, 0.0), at(0.0, WINDOW), hair, ink::line());50    board.polyline(&curve, 3.0, ink::dim());51    let r = unit * 0.36;52    plot::dots(&mut board, &arms, r, ink::blue());53    plot::dots(&mut board, &inner, r, ink::yellow());54    save("wiki-divisor-function", &board)?;55    Ok(())56}