smear.rs

1use crate::lattice::{Family, Rule, FAMILIES};2use num_rational::Ratio;34type Q = Ratio<i128>;56fn whole(value: Q) -> i128 {7    value.floor().to_integer()8}910fn wrap(value: Q) -> Q {11    value - (value / 2).floor() * 212}1314fn side(step: i128, drift: i128, point: Q, grain: i128) -> Vec<([i128; 4], Q)> {15    let mut cuts: Vec<Q> = (0..=2 * grain).map(|j| Ratio::new(j, grain)).collect();16    if step != 0 {17        let span = 4 * grain * (step.abs() + drift.abs() + 2);18        for k in -span..=span {19            let v = (Ratio::new(k, grain) - point * drift) / step;20            if v >= Q::from_integer(0) && v <= Q::from_integer(2) {21                cuts.push(v);22            }23        }24    }25    cuts.sort();26    cuts.dedup();27    let mut out: Vec<([i128; 4], Q)> = Vec::new();28    for pair in cuts.windows(2) {29        if pair[0] == pair[1] {30            continue;31        }32        let mid = (pair[0] + pair[1]) / 2;33        let far = wrap(mid * step + point * drift);34        let key = [35            whole(mid).rem_euclid(2),36            whole(mid * grain) - whole(mid) * grain,37            whole(far).rem_euclid(2),38            whole(far * grain) - whole(far) * grain,39        ];40        let weight = (pair[1] - pair[0]) / 2;41        match out.iter_mut().find(|(seen, _)| *seen == key) {42            Some((_, total)) => *total += weight,43            None => out.push((key, weight)),44        }45    }46    out47}4849fn overtone(residue: i128, p: i128, q: i128) -> i128 {50    let seed = i128::from((p, q) == (0, 0) || (p, q) == (3, 1));51    if residue == 1 {52        seed53    } else {54        1 - seed55    }56}5758fn ink(rule: &Rule, sigma: i128, tau: i128, tone: i128) -> Q {59    let middle = (sigma + tau + tone).rem_euclid(2);60    Q::from_integer(i128::from(rule.filled(61        4 * sigma as i64,62        4 * middle as i64,63        4 * tau as i64,64    )))65}6667fn field(rule: &Rule, x: Q, z: Q, near: i128, far: i128, step: i128, drift: i128) -> [Q; 3] {68    let across = side(step, drift, x, 4);69    let down = side(step, drift, z, 2);70    let mut out = [Q::from_integer(0); 3];71    for (a, wide) in &across {72        for (c, tall) in &down {73            let weight = wide * tall;74            let one = ink(rule, a[0], c[0], overtone(near, a[1], c[1]));75            let two = ink(rule, a[2], c[2], overtone(far, a[3], c[3]));76            out[0] += weight * one;77            out[1] += weight * two;78            out[2] += weight * one * two;79        }80    }81    out82}8384fn clip(poly: &[(Q, Q)], keep: impl Fn(Q, Q) -> Q) -> Vec<(Q, Q)> {85    let zero = Q::from_integer(0);86    let mut out = Vec::new();87    for index in 0..poly.len() {88        let (p, q) = (poly[index], poly[(index + 1) % poly.len()]);89        let (fp, fq) = (keep(p.0, p.1), keep(q.0, q.1));90        if fp >= zero {91            out.push(p);92        }93        if (fp > zero && fq < zero) || (fp < zero && fq > zero) {94            let t = fp / (fp - fq);95            out.push((p.0 + (q.0 - p.0) * t, p.1 + (q.1 - p.1) * t));96        }97    }98    out99}100101fn moments(poly: &[(Q, Q)]) -> [Q; 4] {102    let mut out = [Q::from_integer(0); 4];103    for index in 0..poly.len() {104        let (x0, y0) = poly[index];105        let (x1, y1) = poly[(index + 1) % poly.len()];106        let cross = x0 * y1 - x1 * y0;107        out[0] += cross;108        out[1] += (x0 + x1) * cross;109        out[2] += (y0 + y1) * cross;110        out[3] += (x0 * y1 + x0 * y0 * 2 + x1 * y1 * 2 + x1 * y0) * cross;111    }112    [out[0] / 2, out[1] / 6, out[2] / 6, out[3] / 24]113}114115fn knots(drift: i128, grain: i128) -> Vec<Q> {116    if drift == 0 {117        return vec![Q::from_integer(0), Ratio::new(1, 2), Q::from_integer(1)];118    }119    let cells = grain * drift.abs();120    (0..=cells).map(|j| Ratio::new(j, cells)).collect()121}122123fn limit(rule: &Rule, near: i128, far: i128, step: i128, drift: i128) -> [Q; 3] {124    let across = knots(drift, 4);125    let down = knots(drift, 2);126    let mut totals = [Q::from_integer(0); 3];127    for xi in 0..across.len() - 1 {128        for zi in 0..down.len() - 1 {129            let (x0, x1) = (across[xi], across[xi + 1]);130            let (z0, z1) = (down[zi], down[zi + 1]);131            let cell = vec![(x0, z0), (x1, z0), (x1, z1), (x0, z1)];132            let cell = clip(&cell, |x, z| x + z - Ratio::new(1, 2));133            if cell.len() < 3 {134                continue;135            }136            let cell = clip(&cell, |x, z| Ratio::new(3, 2) - x - z);137            if cell.len() < 3 {138                continue;139            }140            let [area, mx, mz, mxz] = moments(&cell);141            let corner = |x: Q, z: Q| field(rule, x, z, near, far, step, drift);142            let (v00, v01) = (corner(x0, z0), corner(x0, z1));143            let (v10, v11) = (corner(x1, z0), corner(x1, z1));144            let (dx, dz) = (x1 - x0, z1 - z0);145            for slot in 0..3 {146                let c11 = (v11[slot] - v10[slot] - v01[slot] + v00[slot]) / (dx * dz);147                let c10 = (v10[slot] - v00[slot]) / dx - c11 * z0;148                let c01 = (v01[slot] - v00[slot]) / dz - c11 * x0;149                let c00 = v00[slot] - c10 * x0 - c01 * z0 - c11 * x0 * z0;150                totals[slot] += c00 * area + c10 * mx + c01 * mz + c11 * mxz;151            }152        }153    }154    let hexagon = Ratio::new(3, 4);155    [156        totals[0] / hexagon,157        totals[1] / hexagon,158        totals[2] / hexagon,159    ]160}161162const KINDS: [(&str, i128, i128); 4] = [163    ("doubling n = 2m+1", 2, 1),164    ("doubling n = 2m-1", 2, -1),165    ("adjacent n = m+2", 1, 2),166    ("gcd echo n = 3m", 3, 0),167];168169fn far_residue(near: i128, step: i128, drift: i128) -> i128 {170    (near * step + drift).rem_euclid(4)171}172173pub fn run() {174    println!("the pair limit as an exact rational: layer m of residue r mod 4 against layer n = step*m + drift, the phase (mX mod 2, mZ mod 2) equidistributed on the hexagon");175    println!("  the cut cell law s_y = s_x + s_z + w with w = 1 exactly at (p, q) = (0,0) and (3,1) when n = 1 mod 4 and its complement when n = 3 mod 4, integrated over 1/2 <= X + Z <= 3/2 in exact rationals");176    for family in FAMILIES {177        let rule = Rule::new(family);178        for (label, step, drift) in KINDS {179            let mut line = format!("  {} {label}:", family.name());180            for near in [1i128, 3] {181                let far = far_residue(near, step, drift);182                let [em, en, both] = limit(&rule, near, far, step, drift);183                let cov = both - em * en;184                let spread = em * (Q::from_integer(1) - em);185                let same = spread == en * (Q::from_integer(1) - en);186                let r = cov / spread;187                line.push_str(&format!(188                    "  m={near} mod 4: E {em} {en} cov {cov} r {r} = {:+.10}{}",189                    *r.numer() as f64 / *r.denom() as f64,190                    if same { "" } else { " SPREADS DIFFER" }191                ));192            }193            println!("{line}");194        }195    }196    println!("  the four doubling branches read one magnitude, 253/2160, with the sign law -chi4(m) chi4(n); the adjacent and echo carpet rows are -11/135 and 29/135; the tree and void doubling covariances are exactly 0");197    println!("  the exact full-hexagon Pearson against the limit, carpet: the gap times m is bounded, so the finite-layer reading is the limit plus O(1/m)");198    let carpet = Rule::new(Family::Carpet);199    for (m, n, step, drift) in [200        (301usize, 601usize, 2i128, -1i128),201        (601, 1201, 2, -1),202        (103, 205, 2, -1),203        (203, 405, 2, -1),204        (249, 251, 1, 2),205        (99, 297, 3, 0),206    ] {207        let (measured, _) = crate::pairs::pearson(m, n, &carpet);208        let near = (m as i128).rem_euclid(4);209        let [em, en, both] = limit(&carpet, near, far_residue(near, step, drift), step, drift);210        let cov = both - em * en;211        let value = cov / (em * (Q::from_integer(1) - em));212        let target = *value.numer() as f64 / *value.denom() as f64;213        println!(214            "  ({m},{n}): measured {measured:+.8}  limit {value} = {target:+.8}  gap {:.2e}  gap * m {:+.5}",215            (measured - target).abs(),216            (measured - target) * m as f64217        );218    }219}