paper-franel-converse.rs

1use mrlycore::errors::Result;2use mrlycore::Color;3use mrlyfig::{ink, plot, save, Board, Frame};4use mrlynum::design::elements;5use mrlynum::factor::{gcd, mobius_sieve};6use mrlynum::lattice::farey;78const BASE: u64 = 3;9const DIGITS: [u64; 2] = [0, 1];10const ORDER: usize = 81;1112// THE NODES1314fn nodes_of(dens: &[usize]) -> Vec<(usize, usize)> {15    let mut out = Vec::new();16    for &b in dens {17        for a in 1..=b {18            if gcd(a, b) == 1 {19                out.push((a, b));20            }21        }22    }23    out.sort_by(|p, q| (p.0 * q.1).cmp(&(q.0 * p.1)));24    out25}2627fn sawtooth(nodes: &[(usize, usize)]) -> Vec<f64> {28    let m = nodes.len() as f64;29    nodes30        .iter()31        .enumerate()32        .map(|(j, &(a, b))| a as f64 / b as f64 - (j + 1) as f64 / m)33        .collect()34}3536fn kernel(dilates: &[i64]) -> f64 {37    let mut total = 0.0;38    for d in 1..dilates.len() {39        if dilates[d] == 0 {40            continue;41        }42        for e in 1..dilates.len() {43            if dilates[e] == 0 {44                continue;45            }46            let g = gcd(d, e) as f64;47            total += g * g / (d as f64 * e as f64) * (dilates[d] * dilates[e]) as f64;48        }49    }50    total51}5253fn dilated(holds: &[bool], mu: &[i8]) -> Vec<i64> {54    let q = holds.len() - 1;55    let mut out = vec![0i64; q + 1];56    for d in 1..=q {57        for c in 1..=q / d {58            if holds[d * c] {59                out[d] += mu[c] as i64;60            }61        }62    }63    out64}6566// THE MARKS6768fn ticks(board: &mut Board, area: Frame, nodes: &[(usize, usize)], thick: f64, color: Color) {69    let reach = (ORDER as f64).ln();70    let foot = area.y + area.h;71    for &(a, b) in nodes {72        let x = area.x + area.w * a as f64 / b as f64;73        let h = area.h * (1.0 - 0.82 * (b as f64).ln() / reach);74        board.segment((x, foot), (x, foot - h), thick, color);75    }76}7778fn trace(79    board: &mut Board,80    area: Frame,81    nodes: &[(usize, usize)],82    delta: &[f64],83    span: f64,84    thick: f64,85    color: Color,86) {87    let pts: Vec<(f64, f64)> = nodes88        .iter()89        .zip(delta)90        .map(|(&(a, b), &v)| {91            (92                area.x + area.w * a as f64 / b as f64,93                area.y + area.h * (0.5 - 0.5 * v / span),94            )95        })96        .collect();97    board.polyline(&pts, thick, color);98}99100fn main() -> Result<()> {101    let design: Vec<usize> = elements(BASE, &DIGITS, 5)102        .into_iter()103        .filter(|&n| n as usize <= ORDER)104        .map(|n| n as usize)105        .collect();106    assert_eq!(design.len(), 16);107    assert_eq!(design[..5], [1, 3, 4, 9, 10]);108    let mut holds = vec![false; ORDER + 1];109    for &b in &design {110        holds[b] = true;111    }112    let thin = nodes_of(&design);113    assert_eq!(thin.len(), 241);114    let full: Vec<(usize, usize)> = farey(ORDER)115        .iter()116        .filter(|n| n.num > 0)117        .map(|n| (n.num as usize, n.den as usize))118        .collect();119    assert_eq!(full.len(), 2020);120    assert_eq!(full, nodes_of(&(1..=ORDER).collect::<Vec<usize>>()));121122    let delta_thin = sawtooth(&thin);123    let delta_full = sawtooth(&full);124    let s2_thin: f64 = delta_thin.iter().map(|v| v * v).sum();125    let s2_full: f64 = delta_full.iter().map(|v| v * v).sum();126    assert!((s2_thin - 0.004_002_104_850).abs() < 1e-9);127    assert!((s2_full - 0.006_524_653_839).abs() < 1e-9);128129    let mu = mobius_sieve(ORDER);130    let x_thin = dilated(&holds, &mu);131    assert_eq!(x_thin[1], -2);132    assert_eq!(x_thin.iter().map(|v| v * v).sum::<i64>(), 19);133    let g_thin = kernel(&x_thin);134    assert!((g_thin - 12.574_087).abs() < 1e-5);135    assert!((g_thin - (12.0 * 241.0 * s2_thin + 1.0)).abs() < 1e-8);136    let x_full = dilated(&[true; ORDER + 1], &mu);137    assert_eq!(x_full[1], -4);138    let g_full = kernel(&x_full);139    assert!((g_full - 159.157_609).abs() < 1e-5);140    assert!((g_full - (12.0 * 2020.0 * s2_full + 1.0)).abs() < 1e-8);141142    let mut board = Board::square();143    let frame = board.frame(0.08);144    let gap = 26.0;145    let row = (frame.h - 2.0 * gap) * 0.21;146    let top = Frame::new(frame.x, frame.y, frame.w, row);147    let mid = Frame::new(frame.x, frame.y + row + gap, frame.w, row);148    let foot_y = frame.y + 2.0 * (row + gap);149    let foot = Frame::new(frame.x, foot_y, frame.w, frame.y + frame.h - foot_y);150151    plot::baseline(&mut board, top, ink::line());152    plot::baseline(&mut board, mid, ink::line());153    plot::axis(&mut board, foot, ink::line());154    let zero = foot.y + foot.h / 2.0;155    board.segment((foot.x, zero), (foot.x + foot.w, zero), 1.0, ink::line());156157    ticks(&mut board, top.inset(2.0), &thin, 2.2, ink::blue());158    ticks(159        &mut board,160        mid.inset(2.0),161        &full,162        1.0,163        ink::fade(ink::dim(), 0.55),164    );165166    let inner = foot.inset(16.0);167    let span = delta_thin.iter().fold(0.0f64, |a, v| a.max(v.abs())) * 1.12;168    trace(169        &mut board,170        inner,171        &full,172        &delta_full,173        span,174        1.2,175        ink::fade(ink::dim(), 0.7),176    );177    trace(178        &mut board,179        inner,180        &thin,181        &delta_thin,182        span,183        2.4,184        ink::orange(),185    );186    plot::dots(187        &mut board,188        &thin189            .iter()190            .zip(&delta_thin)191            .map(|(&(a, b), &v)| {192                (193                    inner.x + inner.w * a as f64 / b as f64,194                    inner.y + inner.h * (0.5 - 0.5 * v / span),195                )196            })197            .collect::<Vec<_>>(),198        2.6,199        ink::blue(),200    );201    save("paper-franel-converse", &board)?;202    Ok(())203}