pairs.rs

1use crate::lattice::{cell, Family, Rule};2use crate::sums::pairwise;3use mrlynum::factor::gcd;4use mrlynum::series::chi4;56fn marks(first: i64, second: i64, span: i64) -> Vec<i64> {7    let mut out = Vec::new();8    let (mut a, mut b) = (0i64, 0i64);9    while a <= span || b <= span {10        let next = a.min(b);11        out.push(next);12        if a == next {13            a += first;14        }15        if b == next {16            b += second;17        }18    }19    out20}2122pub fn pearson(m: usize, n: usize, rule: &Rule) -> (f64, f64) {23    let g = gcd(m, n) as i64;24    let (unit_m, unit_n) = (n as i64 / g, m as i64 / g);25    let (m, n) = (m as i64, n as i64);26    let (wide, tall) = (4 * m * n / g, 2 * m * n / g);27    let across = marks(unit_m, unit_n, wide);28    let down = marks(unit_m, unit_n, tall);29    let mut strips: Vec<[f64; 4]> = Vec::new();30    for pair in across.windows(2) {31        let width = (pair[1] - pair[0]) as f64 / wide as f64;32        let (xm, xn) = (pair[0] / unit_m, pair[0] / unit_n);33        let mut total = [0.0; 4];34        for strip in down.windows(2) {35            let height = (strip[1] - strip[0]) as f64 / tall as f64;36            let (zm, zn) = (2 * (strip[0] / unit_m), 2 * (strip[0] / unit_n));37            if let (Some(a), Some(b)) = (cell(rule, m, xm, zm), cell(rule, n, xn, zn)) {38                let weight = width * height;39                total[0] += weight;40                total[1] += weight * f64::from(a);41                total[2] += weight * f64::from(b);42                total[3] += weight * f64::from(a && b);43            }44        }45        strips.push(total);46    }47    let sum = |k: usize| pairwise(&strips.iter().map(|row| row[k]).collect::<Vec<f64>>());48    let area = sum(0);49    let (ea, eb, eab) = (sum(1) / area, sum(2) / area, sum(3) / area);50    let covariance = eab - ea * eb;51    (52        covariance / (ea * (1.0 - ea) * eb * (1.0 - eb)).sqrt(),53        covariance,54    )55}5657fn det(a: [[f64; 3]; 3]) -> f64 {58    a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])59        - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])60        + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])61}6263fn richardson(triple: &[(usize, f64)]) -> f64 {64    let row = |m: usize| [1.0, 1.0 / m as f64, 1.0 / (m * m) as f64];65    let a = [row(triple[0].0), row(triple[1].0), row(triple[2].0)];66    let mut b = a;67    for (slot, (_, value)) in triple.iter().enumerate() {68        b[slot][0] = *value;69    }70    det(b) / det(a)71}7273fn branch(first: (usize, f64), second: (usize, f64)) -> f64 {74    (second.0 as f64 * second.1 - first.0 as f64 * first.1) / (second.0 - first.0) as f6475}7677const SMALL: [(usize, usize); 7] = [(5, 9), (5, 7), (9, 17), (13, 25), (17, 51), (5, 15), (3, 9)];7879const DOUBLING: [(usize, usize); 18] = [80    (3, 5),81    (7, 13),82    (7, 15),83    (11, 21),84    (11, 23),85    (19, 39),86    (23, 45),87    (27, 53),88    (103, 205),89    (149, 297),90    (151, 301),91    (157, 313),92    (201, 401),93    (301, 601),94    (401, 801),95    (403, 805),96    (501, 1001),97    (601, 1201),98];99100const FITTED: [usize; 6] = [157, 201, 301, 401, 501, 601];101102type Row = (&'static str, &'static [(usize, usize)], &'static [usize]);103104pub fn small_pairs(rule: &Rule) {105    println!("full-hexagon exact Pearson, carpet, the seven blind pairs");106    for (m, n) in SMALL {107        let (r, covariance) = pearson(m, n, rule);108        println!("  ({m},{n}) r = {r:+.8}  cov = {covariance:+.8}");109    }110}111112pub fn doubling(rule: &Rule) {113    println!("doubling pairs (m, 2m+-1), carpet: sign law sign r = -chi4(m) chi4(n)");114    let mut hits = 0;115    let mut fitted: Vec<(usize, f64)> = Vec::new();116    let mut fitted_cov: Vec<(usize, f64)> = Vec::new();117    let mut branches: Vec<(usize, f64)> = Vec::new();118    for (m, n) in DOUBLING {119        let (r, covariance) = pearson(m, n, rule);120        let predicted = -f64::from(chi4(m) * chi4(n));121        let agreed = r.signum() == predicted;122        hits += usize::from(agreed);123        println!(124            "  ({m},{n}) r = {r:+.8}  cov = {covariance:+.8}  predicted sign {predicted:+}  {}",125            if agreed { "ok" } else { "MISS" }126        );127        if FITTED.contains(&m) {128            fitted.push((m, r));129            fitted_cov.push((m, covariance));130        }131        if [103, 403, 501, 601].contains(&m) {132            branches.push((m, r.abs()));133        }134    }135    println!("  sign law: {hits}/{} pairs", DOUBLING.len());136    println!("Richardson r = r_inf + a/m + b/m^2 on sliding triples of m = {FITTED:?}");137    for start in 0..FITTED.len() - 2 {138        let r = richardson(&fitted[start..start + 3]);139        let c = richardson(&fitted_cov[start..start + 3]);140        println!(141            "  m = {:?}: r_inf = {r:+.8}  cov_inf = {c:+.8}",142            &FITTED[start..start + 3]143        );144    }145    let low = branch(branches[2], branches[3]);146    let high = branch(branches[0], branches[1]);147    println!("  branch 1/m extrapolation |r|: (501,601) -> {low:.7}  (103,403) -> {high:.7}");148    println!(149        "  candidates: 253/2160 = {:.8}  19/162 = {:.8}",150        253.0 / 2160.0,151        19.0 / 162.0152    );153}154155pub fn persistence() {156    println!("other pairs and families");157    let rules: Vec<(Family, Rule)> = [Family::Carpet, Family::Tree, Family::Void]158        .into_iter()159        .map(|family| (family, Rule::new(family)))160        .collect();161    let table: [Row; 4] = [162        ("doubling", &[(201, 401)], &[1, 2]),163        ("adjacent", &[(199, 201), (249, 251)], &[0, 1, 2]),164        ("gcd echo (m,3m)", &[(67, 201), (99, 297)], &[0, 1, 2]),165        ("coprime control", &[(101, 173), (97, 251)], &[0, 1, 2]),166    ];167    for (label, pairs, chosen) in table {168        for (m, n) in pairs.iter().copied() {169            let mut line = format!("  {label} ({m},{n}):");170            for slot in chosen {171                let (family, rule) = &rules[*slot];172                let (r, _) = pearson(m, n, rule);173                line.push_str(&format!("  {} {r:+.8}", family.name()));174            }175            println!("{line}");176        }177    }178    println!(179        "  candidates: adjacent -11/135 = {:+.7}  echo 29/135 = {:+.7}",180        -11.0 / 135.0,181        29.0 / 135.0182    );183}