paper-sparse-mertens-under-grh.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, save, Board, Ramp};34const BASE: i64 = 7;5const DROP: i64 = 3;6const LEVEL: u32 = 3;78fn kept() -> Vec<i64> {9    (0..BASE).filter(|d| *d != DROP).collect()10}1112fn symbol(kept: &[i64], t: f64) -> f64 {13    let (mut re, mut im) = (0.0, 0.0);14    for d in kept {15        let angle = std::f64::consts::TAU * (*d as f64) * t;16        re += angle.cos();17        im += angle.sin();18    }19    re.hypot(im)20}2122fn transform(kept: &[i64], a: i64, span: i64) -> f64 {23    let mut mass = 1.0;24    let mut scale = 1i64;25    for _ in 0..LEVEL {26        let t = (a * scale).rem_euclid(span) as f64 / span as f64;27        mass *= symbol(kept, t);28        scale *= BASE;29    }30    mass31}3233fn proved() -> f64 {34    let q = BASE as f64;35    let n = ((BASE - 2) as f64 / 2.0).ceil();36    let harmonic = n.ln() + 0.577_215_664_901_532_9 + 1.0 / (2.0 * n);37    let pi = std::f64::consts::PI;38    let phi = (4.0 / pi) * q + (2.0 * q / pi) * harmonic + (1.0 - 2.0 / pi) * (q - 2.0) + 0.727;39    1.0 + phi / q40}4142fn main() -> Result<()> {43    let kept = kept();44    let span = BASE.pow(LEVEL);45    let peak = (kept.len() as f64).powi(LEVEL as i32);46    let half = (span - 1) / 2;47    let mass: Vec<f64> = (-half..=half).map(|a| transform(&kept, a, span)).collect();4849    assert_eq!(mass.len(), span as usize);50    assert_eq!(kept.len() as i64, BASE - 1);51    assert!((mass[half as usize] - peak).abs() < 1e-9);52    let energy: f64 = mass.iter().map(|v| v * v).sum();53    assert!((energy - span as f64 * peak).abs() < 1e-6 * span as f64 * peak);54    let total: f64 = mass.iter().sum();55    assert!(total < (BASE as f64 * proved()).powi(LEVEL as i32));56    assert!(total >= span as f64);5758    let mut board = Board::square();59    let frame = board.frame(0.07);60    let ramp = Ramp::tone(ink::blue(), ink::yellow());61    let slot = frame.w / mass.len() as f64;62    let pad = slot * 0.18;63    let axis = frame.y + frame.h / 2.0;64    for (i, value) in mass.iter().enumerate() {65        let reach = frame.h / 2.0 * (value / peak).powf(1.0 / LEVEL as f64);66        board.rect(67            frame.x + i as f64 * slot + pad,68            axis - reach,69            slot - 2.0 * pad,70            2.0 * reach,71            ramp.at(2.0 * reach / frame.h),72        );73    }74    save("paper-sparse-mertens-under-grh", &board)?;75    Ok(())76}