coprime.rs

1use mrlymath::bang::factory::{code_to_corners, corners_to_code, residue_corners};2use mrlymath::bang::Code;3use mrlynum::factor::{divisors, factorize, gcd, mobius, radical};4use mrlynum::series::zeta;5use std::collections::BTreeSet;67const BUDGET: usize = 200_000;8const TOLERANCE: f64 = 0.06;9const EXACT_LEVEL: u32 = 4;10const ZETA_TERMS: usize = 10_000;11const HEAD: usize = 6;1213pub struct Line {14    pub code: Code,15    pub k: usize,16    pub index: usize,17    pub spanning: bool,18    pub bracket: f64,19    pub predicted: f64,20    pub measured: f64,21    pub level: u32,22    pub exact: bool,23    pub within: bool,24    pub monotone: bool,25    pub head: Vec<u64>,26}2728pub struct Family {29    pub label: String,30    pub base: usize,31    pub lines: Vec<Line>,32}3334impl Family {35    fn spanning(&self) -> usize {36        self.lines.iter().filter(|line| line.spanning).count()37    }3839    fn distinct(&self, spanning_only: bool) -> usize {40        let heads: BTreeSet<&Vec<u64>> = self41            .lines42            .iter()43            .filter(|line| line.spanning || !spanning_only)44            .map(|line| &line.head)45            .collect();46        heads.len()47    }4849    fn flagged(&self) -> usize {50        self.lines51            .iter()52            .filter(|line| line.spanning && (!line.within || !line.exact))53            .count()54    }5556    fn exact_failures(&self) -> usize {57        self.lines.iter().filter(|line| !line.exact).count()58    }59}6061fn squarefree_divisors(base: usize) -> Vec<usize> {62    divisors(radical(base))63}6465fn hits(corners: &[Vec<u8>], divisor: usize) -> usize {66    corners67        .iter()68        .filter(|corner| corner.iter().all(|&r| (r as usize).is_multiple_of(divisor)))69        .count()70}7172fn bracket(corners: &[Vec<u8>], base: usize) -> f64 {73    let sum: f64 = squarefree_divisors(base)74        .iter()75        .map(|&e| f64::from(mobius(e)) * hits(corners, e) as f64)76        .sum();77    sum / corners.len() as f6478}7980fn foreign(base: usize, dimension: usize) -> f64 {81    let mut out = 1.0 / zeta(dimension as f64, ZETA_TERMS);82    for (prime, _) in factorize(base) {83        out /= 1.0 - (prime as f64).powi(-(dimension as i32));84    }85    out86}8788fn determinant(rows: &[Vec<i64>]) -> i64 {89    if rows.len() == 1 {90        return rows[0][0];91    }92    let mut total = 0;93    for (column, head) in rows[0].iter().enumerate() {94        let minor: Vec<Vec<i64>> = rows[1..]95            .iter()96            .map(|row| {97                row.iter()98                    .enumerate()99                    .filter(|(i, _)| *i != column)100                    .map(|(_, v)| *v)101                    .collect()102            })103            .collect();104        let sign = if column % 2 == 0 { 1 } else { -1 };105        total += sign * head * determinant(&minor);106    }107    total108}109110fn choose(start: usize, count: usize, size: usize) -> Vec<Vec<usize>> {111    if size == 0 {112        return vec![Vec::new()];113    }114    let mut out = Vec::new();115    for first in start..count {116        for mut rest in choose(first + 1, count, size - 1) {117            rest.insert(0, first);118            out.push(rest);119        }120    }121    out122}123124fn index(corners: &[Vec<u8>], dimension: usize) -> usize {125    let first = &corners[0];126    let diffs: Vec<Vec<i64>> = corners[1..]127        .iter()128        .map(|corner| {129            corner130                .iter()131                .zip(first)132                .map(|(a, b)| i64::from(*a) - i64::from(*b))133                .collect()134        })135        .collect();136    let mut out = 0usize;137    for rows in choose(0, diffs.len(), dimension) {138        let minor: Vec<Vec<i64>> = rows.iter().map(|&r| diffs[r].clone()).collect();139        out = gcd(out, determinant(&minor).unsigned_abs() as usize);140        if out == 1 {141            return out;142        }143    }144    out145}146147fn expand(columns: &mut [Vec<usize>], corners: &[Vec<u8>], base: usize) {148    for (axis, column) in columns.iter_mut().enumerate() {149        let mut next = Vec::with_capacity(column.len() * corners.len());150        for corner in corners {151            next.extend(column.iter().map(|&x| base * x + corner[axis] as usize));152        }153        *column = next;154    }155}156157fn common(columns: &[Vec<usize>]) -> Vec<usize> {158    (0..columns[0].len())159        .map(|i| columns.iter().fold(0, |g, column| gcd(g, column[i])))160        .collect()161}162163fn measure(corners: &[Vec<u8>], base: usize, dimension: usize) -> (Vec<u64>, bool) {164    let k = corners.len();165    let mut columns = vec![vec![0usize]; dimension];166    let mut terms = Vec::new();167    let mut exact = None;168    let mut level = 0u32;169    while k.pow(level + 1) <= BUDGET {170        level += 1;171        expand(&mut columns, corners, base);172        let shared = common(&columns);173        terms.push(shared.iter().filter(|&&g| g == 1).count() as u64);174        if level == EXACT_LEVEL {175            exact = Some(176                squarefree_divisors(base)177                    .iter()178                    .filter(|&&e| e > 1)179                    .all(|&e| {180                        let seen = shared.iter().filter(|&&g| g % e == 0).count();181                        seen == hits(corners, e) * k.pow(level - 1)182                    }),183            );184        }185    }186    (187        terms,188        exact.expect("the point budget reaches the exact level"),189    )190}191192pub fn survey(193    base: usize,194    dimension: usize,195    codes: impl Iterator<Item = Code>,196    label: &str,197) -> Family {198    let mut lines = Vec::new();199    for code in codes {200        let corners = code_to_corners(code, dimension, base).expect("the code fits its cells");201        let k = corners.len();202        if k < 2 {203            continue;204        }205        let (terms, exact) = measure(&corners, base, dimension);206        let level = terms.len() as u32;207        let measured = terms[terms.len() - 1] as f64 / k.pow(level) as f64;208        let bracket = bracket(&corners, base);209        let predicted = bracket * foreign(base, dimension);210        let gaps: Vec<f64> = terms211            .iter()212            .enumerate()213            .map(|(step, &count)| (count as f64 / k.pow(step as u32 + 1) as f64 - predicted).abs())214            .collect();215        let monotone = gaps.windows(2).all(|pair| pair[1] <= pair[0]);216        let index = index(&corners, dimension);217        lines.push(Line {218            code,219            k,220            index,221            spanning: index == 1,222            bracket,223            predicted,224            measured,225            level,226            exact,227            within: (measured - predicted).abs() < TOLERANCE,228            monotone,229            head: terms.iter().take(HEAD).copied().collect(),230        });231    }232    Family {233        label: label.to_string(),234        base,235        lines,236    }237}238239fn sponge() -> Code {240    let filled: Vec<Vec<u8>> = residue_corners(3, 3)241        .into_iter()242        .filter(|corner| corner.iter().filter(|&&r| r == 1).count() <= 1)243        .collect();244    corners_to_code(&filled, 3, 3)245}246247fn base6() -> Vec<Code> {248    let samples: [[[u8; 2]; 8]; 2] = [249        [250            [0, 0],251            [1, 1],252            [2, 3],253            [3, 2],254            [4, 5],255            [5, 4],256            [1, 0],257            [0, 1],258        ],259        [260            [0, 0],261            [2, 0],262            [4, 0],263            [0, 3],264            [1, 1],265            [5, 5],266            [1, 2],267            [2, 1],268        ],269    ];270    samples271        .iter()272        .map(|sample| {273            let filled: Vec<Vec<u8>> = sample.iter().map(|corner| corner.to_vec()).collect();274            corners_to_code(&filled, 2, 6)275        })276        .collect()277}278279pub fn report() {280    let families = [281        survey(2, 2, 1..16, "base 2 D 2"),282        survey(2, 3, 1..256, "base 2 D 3"),283        survey(3, 2, 1..512, "base 3 D 2"),284        survey(3, 3, [sponge()].into_iter(), "menger sponge"),285        survey(6, 2, base6().into_iter(), "base 6 samples"),286    ];287    println!("coprime census: designs with k >= 2, point budget {BUDGET}, tolerance {TOLERANCE}, exact identity at n = {EXACT_LEVEL}");288    println!(289        "{:<16}{:>9}{:>10}{:>10}{:>19}{:>9}{:>13}",290        "family", "designs", "spanning", "distinct", "distinct spanning", "flagged", "exact fails"291    );292    for family in &families {293        println!(294            "{:<16}{:>9}{:>10}{:>10}{:>19}{:>9}{:>13}",295            family.label,296            family.lines.len(),297            family.spanning(),298            family.distinct(false),299            family.distinct(true),300            family.flagged(),301            family.exact_failures()302        );303    }304    println!(305        "{:<16}{:>9}{:>10}{:>10}{:>19}{:>9}{:>13}",306        "total",307        families.iter().map(|f| f.lines.len()).sum::<usize>(),308        families.iter().map(Family::spanning).sum::<usize>(),309        "-",310        "-",311        families.iter().map(Family::flagged).sum::<usize>(),312        families.iter().map(Family::exact_failures).sum::<usize>()313    );314    let by = |keep: &dyn Fn(&Line, usize) -> bool| {315        families316            .iter()317            .flat_map(|f| f.lines.iter().map(move |line| (line, f.base)))318            .filter(|(line, base)| keep(line, *base))319            .count()320    };321    println!(322        "spanning by dimension log_q(k): above {}, exactly one {}, below {}; index 3 at k = q {}",323        by(&|line, base| line.spanning && line.k > base),324        by(&|line, base| line.spanning && line.k == base),325        by(&|line, base| line.spanning && line.k < base),326        by(&|line, base| line.index == 3 && line.k == base)327    );328    println!(329        "spanning lines whose gap to the predicted density widens at some level: {}",330        by(&|line, _| line.spanning && !line.monotone)331    );332    for family in &families[3..] {333        for line in &family.lines {334            let corners =335                code_to_corners(line.code, 2 + usize::from(family.base == 3), family.base)336                    .expect("the named code fits its cells");337            let counts: Vec<String> = squarefree_divisors(family.base)338                .iter()339                .filter(|&&e| e > 1)340                .map(|&e| format!("k_{e} {}", hits(&corners, e)))341                .collect();342            println!(343                "{}: code {} k {} {} bracket {:.4} pred {:.6} meas {:.6} n {} terms {:?}",344                family.label,345                line.code,346                line.k,347                counts.join(" "),348                line.bracket,349                line.predicted,350                line.measured,351                line.level,352                line.head353            );354        }355    }356    for k in [2usize, 3] {357        let level = families358            .iter()359            .flat_map(|f| f.lines.iter())360            .find(|line| line.k == k)361            .map(|line| line.level)362            .expect("a design with that k exists");363        println!("deepest level at k = {k}: n = {level}");364    }365}