main.rs

1mod design;23use std::collections::BTreeSet;45use mrlynum::factor::gcd;6use mrlynum::lattice::{farey, new_nodes, nodes, totients};78const QS: [usize; 7] = [125, 250, 500, 1000, 2000, 4000, 8000];910const SMALL: [usize; 4] = [10, 30, 60, 125];1112fn literal_stack(q: usize) -> BTreeSet<(usize, usize)> {13    let mut lit = BTreeSet::new();14    for n in 1..=q {15        for k in 1..=n {16            let g = gcd(k, n);17            lit.insert((k / g, n / g));18        }19    }20    lit21}2223fn window_pairs(q: usize) -> BTreeSet<(usize, usize)> {24    nodes(q)25        .into_iter()26        .filter(|node| node.num > 0)27        .map(|node| (node.num as usize, node.den as usize))28        .collect()29}3031fn walk_pairs(q: usize) -> BTreeSet<(usize, usize)> {32    farey(q)33        .into_iter()34        .filter(|node| node.num > 0)35        .map(|node| (node.num as usize, node.den as usize))36        .collect()37}3839fn drawn_by(q: usize, pair: (usize, usize)) -> usize {40    let (a, b) = pair;41    (1..=q).filter(|n| (a * n) % b == 0).count()42}4344fn brightness_ok(q: usize) -> bool {45    nodes(q)46        .into_iter()47        .filter(|node| node.num > 0)48        .all(|node| {49            let b = node.den as usize;50            node.brightness == (q / b) as u64 && drawn_by(q, (node.num as usize, b)) == q / b51        })52}5354fn totient_sum(q: usize, phi: &[u64]) -> u64 {55    phi[1..=q].iter().sum()56}5758fn meter_walk(q: usize, m: u64) -> (f64, f64, u64) {59    let mut s1 = 0.0f64;60    let mut s2 = 0.0f64;61    let mut j = 0u64;62    for node in farey(q) {63        if node.num == 0 {64            continue;65        }66        j += 1;67        let delta = node.num as f64 / node.den as f64 - j as f64 / m as f64;68        s1 += delta.abs();69        s2 += delta * delta;70    }71    (s1, s2, j)72}7374fn meter_window(q: usize, m: u64) -> (f64, f64, u64) {75    let mut s1 = 0.0f64;76    let mut s2 = 0.0f64;77    let mut j = 0u64;78    for node in nodes(q) {79        if node.num == 0 {80            continue;81        }82        j += 1;83        let delta = node.num as f64 / node.den as f64 - j as f64 / m as f64;84        s1 += delta.abs();85        s2 += delta * delta;86    }87    (s1, s2, j)88}8990fn verdict(ok: bool) -> &'static str {91    if ok {92        "PASS"93    } else {94        "FAIL"95    }96}9798fn stack() {99    let phi = totients(*QS.last().unwrap());100    println!("LIT NODES ARE THE FAREY NODES");101    println!("      Q   lit set   window   mediant walk   sum phi(k)   floor(Q/b)   agree");102    for q in SMALL {103        let lit = literal_stack(q);104        let win = window_pairs(q);105        let walk = walk_pairs(q);106        let control = totient_sum(q, &phi);107        let bright = brightness_ok(q);108        let ok = lit == win && win == walk && lit.len() as u64 == control && bright;109        println!(110            "  {:>5}   {:>7}   {:>6}   {:>12}   {:>10}   {:>10}   {}",111            q,112            lit.len(),113            win.len(),114            walk.len(),115            control,116            verdict(bright),117            verdict(ok)118        );119    }120121    let fresh = (2..=60).all(|n| new_nodes(n) == phi[n]);122    println!(123        "  new nodes at scale n equals phi(n), n = 2..60 : {}",124        verdict(fresh)125    );126}127128fn meter() {129    let phi = totients(*QS.last().unwrap());130    println!("THE DISCREPANCY METER");131    println!("      Q       nodes   sum phi(k)     S2*Q   S1/sqrt(Q)   exp S2   exp S1");132    let mut prev: Option<(usize, f64, f64)> = None;133    for q in QS {134        let m = totient_sum(q, &phi);135        let (s1, s2, seen) = meter_walk(q, m);136        let (mut e2, mut e1) = ("-".to_string(), "-".to_string());137        if let Some((pq, p1, p2)) = prev {138            let span = (q as f64 / pq as f64).ln();139            e2 = format!("{:+.3}", (s2 / p2).ln() / span);140            e1 = format!("{:+.3}", (s1 / p1).ln() / span);141        }142        prev = Some((q, s1, s2));143        println!(144            "  {:>5}  {:>10}  {:>11}   {:.4}       {:.4}   {:>6}   {:>6}",145            q,146            seen,147            m,148            s2 * q as f64,149            s1 / (q as f64).sqrt(),150            e2,151            e1152        );153    }154155    println!();156    println!("CROSS-CHECK  the sorted window route against the mediant walk");157    println!("      Q       nodes     S2*Q   S1/sqrt(Q)   agrees");158    for q in [125usize, 250, 500, 1000] {159        let m = totient_sum(q, &phi);160        let (w1, w2, wn) = meter_window(q, m);161        let (r1, r2, rn) = meter_walk(q, m);162        let ok = wn == rn && wn == m && (w1 - r1).abs() < 1e-9 && (w2 - r2).abs() < 1e-12;163        println!(164            "  {:>5}  {:>10}   {:.4}       {:.4}   {}",165            q,166            wn,167            w2 * q as f64,168            w1 / (q as f64).sqrt(),169            verdict(ok)170        );171    }172}173174fn main() {175    let want = std::env::args().nth(1);176    let pick = |name: &str| match want.as_deref() {177        Some(verb) => verb == name,178        None => true,179    };180    if pick("stack") {181        stack();182        println!();183    }184    if pick("meter") {185        meter();186        println!();187    }188    if pick("design") {189        design::run();190    }191}