wiki-mertens-function.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board};3use mrlynum::factor::mobius_sieve;45const TOP: usize = 1000;67fn main() -> Result<()> {8    let mut board = Board::square();9    let frame = board.frame(0.08);10    let reach = (TOP as f64).sqrt() * 1.1;11    let across = |n: f64| frame.x + frame.w * (n - 1.0) / (TOP as f64 - 1.0);12    let up = |v: f64| frame.y + frame.h * (1.0 - (v + reach) / (2.0 * reach));1314    let mu = mobius_sieve(TOP);15    let mut walk: Vec<i64> = Vec::with_capacity(TOP);16    let mut sum = 0i64;17    for n in 1..=TOP {18        sum += i64::from(mu[n]);19        walk.push(sum);20    }2122    let mut over = Vec::with_capacity(TOP);23    let mut under = Vec::with_capacity(TOP);24    for n in 1..=TOP {25        let root = (n as f64).sqrt();26        over.push((across(n as f64), up(root)));27        under.push((across(n as f64), up(-root)));28    }29    board.polyline(&over, 2.6, ink::dim());30    board.polyline(&under, 2.6, ink::dim());31    board.segment(32        (frame.x, up(0.0)),33        (frame.x + frame.w, up(0.0)),34        1.6,35        ink::line(),36    );3738    let mut steps = Vec::with_capacity(2 * TOP);39    for (i, &m) in walk.iter().enumerate() {40        let n = (i + 1) as f64;41        steps.push((across(n), up(m as f64)));42        steps.push((across(n + 1.0), up(m as f64)));43    }44    board.polyline(&steps, 2.6, ink::blue());4546    plot::axis(&mut board, frame, ink::line());47    assert_eq!(walk.len(), TOP);48    assert_eq!(walk[999], 2);49    assert_eq!(walk.iter().copied().min(), Some(-12));50    assert!(walk51        .iter()52        .enumerate()53        .all(|(i, &m)| (m as f64).abs() <= ((i + 1) as f64).sqrt()));54    save("wiki-mertens-function", &board)?;55    Ok(())56}