lattice.rs

1use crate::factor::coprime;2use crate::series;3use std::f64::consts::PI;45/// Sieves the Euler totients of zero through n.6///7/// ```8/// assert_eq!(mrlynum::lattice::totients(6), vec![0, 1, 1, 2, 2, 4, 2]);9/// ```10pub fn totients(n: usize) -> Vec<u64> {11    let mut phi: Vec<u64> = (0..=n as u64).collect();12    for p in 2..=n {13        if phi[p] == p as u64 {14            for m in (p..=n).step_by(p) {15                phi[m] -= phi[m] / p as u64;16            }17        }18    }19    phi20}2122/// Counts the ordered pairs of coprime coordinates between one and n: twice the totient sum less one.23pub fn coprime_pairs(n: usize) -> u64 {24    if n == 0 {25        return 0;26    }27    2 * totients(n)[1..].iter().sum::<u64>() - 128}2930/// Estimates pi from visibility: the density of coprime pairs in the n-by-n window tends to six over pi squared.31pub fn pi_estimate(n: usize) -> f64 {32    let density = coprime_pairs(n) as f64 / (n as f64 * n as f64);33    (6.0 / density).sqrt()34}3536/// The rational factor r with zeta of the dimension equal to r times pi to the dimension, read off the Bernoulli fraction; none at an odd dimension or past twelve.37///38/// ```39/// assert!((mrlynum::lattice::zeta_factor(4).unwrap() - 1.0 / 90.0).abs() < 1e-15);40/// ```41pub fn zeta_factor(dimension: u32) -> Option<f64> {42    if dimension == 0 || !dimension.is_multiple_of(2) || dimension > 12 {43        return None;44    }45    let d = dimension as usize;46    let (num, den) = series::bernoulli(d + 1)[d];47    let sign = if (d / 2).is_multiple_of(2) { -1.0 } else { 1.0 };48    let factorial = (1..=d).fold(1.0f64, |out, k| out * k as f64);49    Some(sign * (num as f64 / den as f64) * 2f64.powi(d as i32) / (2.0 * factorial))50}5152/// The value zeta takes at a whole argument above one: the exact Bernoulli form at an even one, the Euler-Maclaurin sum at an odd one.53///54/// Panics at a whole argument of one or below, where the sum does not converge.55pub fn zeta_whole(s: u32) -> f64 {56    assert!(s > 1, "zeta needs a whole argument above one");57    match zeta_factor(s) {58        Some(factor) => factor * PI.powi(s as i32),59        None => series::zeta(f64::from(s), 20_000),60    }61}6263/// The density the visible count of a window in the dimension walks to: one over zeta of the dimension.64pub fn visible_density(dimension: u32) -> f64 {65    1.0 / zeta_whole(dimension)66}6768/// Recovers the constant the dimension hides from the visible count of the window: pi at an even dimension, zeta of the dimension at an odd one.69pub fn recovered(n: usize, dimension: u32) -> f64 {70    let density = series::visible(n, dimension) as f64 / (n as f64).powi(dimension as i32);71    let zeta = 1.0 / density;72    match zeta_factor(dimension) {73        Some(factor) => (zeta / factor).powf(1.0 / f64::from(dimension)),74        None => zeta,75    }76}7778/// A visible node: a reduced fraction and the brightness a stack of scales one through the window gives it.79#[derive(Clone, Copy, Debug, PartialEq, Eq)]80pub struct Node {81    /// The numerator, coprime to the denominator.82    pub num: u64,83    /// The denominator.84    pub den: u64,85    /// The count of scales putting a line here: the floor of the window over the denominator.86    pub brightness: u64,87}8889/// A grid crossing of two visible nodes, its brightness the separable product.90#[derive(Clone, Copy, Debug, PartialEq, Eq)]91pub struct Node2d {92    /// The horizontal node.93    pub x: Node,94    /// The vertical node.95    pub y: Node,96    /// The product of the two axis brightnesses.97    pub brightness: u64,98}99100/// Lists the visible nodes of a window in ascending value: every reduced fraction with denominator at most n.101pub fn nodes(n: usize) -> Vec<Node> {102    let mut out = Vec::new();103    for den in 1..=n {104        for num in 0..=den {105            if coprime(num, den) {106                out.push(Node {107                    num: num as u64,108                    den: den as u64,109                    brightness: (n / den) as u64,110                });111            }112        }113    }114    out.sort_by(|a, b| (a.num as u128 * b.den as u128).cmp(&(b.num as u128 * a.den as u128)));115    out116}117118/// Walks the Farey sequence of the order by the Stern-Brocot mediant recurrence from zero over one to one over one, an independent route to the same nodes.119pub fn farey(order: usize) -> Vec<Node> {120    let mut out = Vec::new();121    if order == 0 {122        return out;123    }124    let n = order as u64;125    let (mut a, mut b, mut c, mut d) = (0u64, 1u64, 1u64, n);126    out.push(Node {127        num: a,128        den: b,129        brightness: n / b,130    });131    while c <= n {132        let k = (n + b) / d;133        (a, b, c, d) = (c, d, k * c - a, k * d - b);134        out.push(Node {135            num: a,136            den: b,137            brightness: n / b,138        });139    }140    out141}142143/// Lists the grid crossings of a window's nodes, row-major over the ascending axis nodes.144pub fn grid(n: usize) -> Vec<Node2d> {145    let axis = nodes(n);146    let mut out = Vec::with_capacity(axis.len() * axis.len());147    for &y in &axis {148        for &x in &axis {149            out.push(Node2d {150                x,151                y,152                brightness: x.brightness * y.brightness,153            });154        }155    }156    out157}158159/// Counts the nodes window n lights that window n minus one lacked: two at window one, phi of n after.160pub fn new_nodes(n: usize) -> u64 {161    if n == 0 {162        return 0;163    }164    (nodes(n).len() - nodes(n - 1).len()) as u64165}166167#[cfg(test)]168mod tests {169    use super::*;170    use std::f64::consts::PI;171172    #[test]173    fn totients_match_hand_checked_values() {174        let phi = totients(12);175        assert_eq!(phi, vec![0, 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4]);176    }177178    #[test]179    fn coprime_pairs_match_a_brute_count() {180        for n in [1usize, 2, 3, 10, 20] {181            let mut brute = 0;182            for a in 1..=n {183                for b in 1..=n {184                    if coprime(a, b) {185                        brute += 1;186                    }187                }188            }189            assert_eq!(coprime_pairs(n), brute, "window {n}");190        }191    }192193    #[test]194    fn pi_estimate_converges_by_a_hundred_thousand() {195        assert!((pi_estimate(100_000) - PI).abs() < 1e-4);196    }197198    #[test]199    fn the_zeta_factor_is_the_known_even_fraction() {200        assert!((zeta_factor(2).unwrap() - 1.0 / 6.0).abs() < 1e-15);201        assert!((zeta_factor(4).unwrap() - 1.0 / 90.0).abs() < 1e-15);202        assert!((zeta_factor(6).unwrap() - 1.0 / 945.0).abs() < 1e-15);203        assert_eq!(zeta_factor(3), None);204        assert_eq!(zeta_factor(14), None);205    }206207    #[test]208    fn the_visible_density_is_one_over_the_whole_zeta() {209        assert!((visible_density(2) - 6.0 / (PI * PI)).abs() < 1e-15);210        assert!((zeta_whole(2) - PI * PI / 6.0).abs() < 1e-15);211        assert!((zeta_whole(3) - 1.202_056_903_159_594).abs() < 1e-9);212    }213214    #[test]215    fn the_window_recovers_pi_in_the_even_dimensions() {216        assert!((recovered(1_000, 2) - pi_estimate(1_000)).abs() < 1e-12);217        assert!((recovered(1_000, 2) - PI).abs() < 2e-3);218        assert!((recovered(1_000, 4) - PI).abs() < 2e-3);219        assert!((recovered(1_000, 3) - 1.202_056_903).abs() < 2e-3);220    }221222    #[test]223    fn node_counts_follow_the_farey_sequence() {224        for (n, count) in [(1, 2), (2, 3), (3, 5), (4, 7), (5, 11), (6, 13)] {225            assert_eq!(nodes(n).len(), count, "window {n}");226        }227    }228229    #[test]230    fn nodes_are_reduced_bright_and_ascending() {231        let all = nodes(12);232        for pair in all.windows(2) {233            assert!(pair[0].num * pair[1].den < pair[1].num * pair[0].den);234        }235        for node in all {236            assert!(coprime(node.num as usize, node.den as usize));237            assert_eq!(node.brightness, 12 / node.den);238        }239    }240241    #[test]242    fn the_farey_walk_lands_on_the_window_nodes_exactly() {243        assert!(farey(0).is_empty());244        for n in 1..=50 {245            assert_eq!(farey(n), nodes(n), "order {n}");246        }247        let ends = farey(7);248        assert_eq!((ends[0].num, ends[0].den), (0, 1));249        assert_eq!((ends.last().unwrap().num, ends.last().unwrap().den), (1, 1));250    }251252    #[test]253    fn the_farey_length_is_one_past_the_totient_sum() {254        let phi = totients(50);255        for n in 1..=50usize {256            let want = 1 + phi[1..=n].iter().sum::<u64>() as usize;257            assert_eq!(farey(n).len(), want, "order {n}");258        }259    }260261    #[test]262    fn grid_brightness_is_separable() {263        let flat = grid(3);264        assert_eq!(flat.len(), 25);265        for cross in flat {266            assert_eq!(cross.brightness, cross.x.brightness * cross.y.brightness);267        }268        let corner = &grid(3)[0];269        assert_eq!((corner.x.num, corner.y.num), (0, 0));270        assert_eq!(corner.brightness, 9);271    }272273    #[test]274    fn new_nodes_is_the_totient() {275        let phi = totients(50);276        assert_eq!(new_nodes(1), 2);277        for (n, &expected) in phi.iter().enumerate().skip(2) {278            assert_eq!(new_nodes(n), expected, "window {n}");279        }280    }281}