counting.rs

1use mrlycore::errors::{value_error, Result};2use mrlynum::classics::{factorial, gcd};3use mrlynum::factor::{divisors, mobius};4use std::collections::BTreeMap;56type Cycles = BTreeMap<usize, u128>;78fn pos_block_cycles(length: usize) -> Cycles {9    let mut out = Cycles::new();10    for period in divisors(length) {11        let strings: i128 = divisors(period)12            .iter()13            .map(|&d| i128::from(mobius(period / d)) * (1i128 << d))14            .sum();15        *out.entry(period).or_insert(0) += (strings / period as i128) as u128;16    }17    out18}1920fn neg_block_cycles(length: usize) -> Cycles {21    let states = 1usize << length;22    let step = |x: usize| -> usize {23        let mut new = vec![0usize; length];24        for (i, item) in new.iter_mut().enumerate() {25            *item = (x >> ((i + length - 1) % length)) & 1;26        }27        new[0] ^= 1;28        new.iter().enumerate().map(|(i, &b)| b << i).sum()29    };30    let mut seen = vec![false; states];31    let mut out = Cycles::new();32    for start in 0..states {33        if seen[start] {34            continue;35        }36        let mut run = 0;37        let mut j = start;38        while !seen[j] {39            seen[j] = true;40            j = step(j);41            run += 1;42        }43        *out.entry(run).or_insert(0) += 1;44    }45    out46}4748fn combine(c1: &Cycles, c2: &Cycles) -> Cycles {49    let mut out = Cycles::new();50    for (&l1, &n1) in c1 {51        for (&l2, &n2) in c2 {52            let g = gcd(l1 as u128, l2 as u128) as usize;53            *out.entry(l1 * l2 / g).or_insert(0) += n1 * n2 * g as u128;54        }55    }56    out57}5859fn class_cycles(pos: &[usize], neg: &[usize]) -> u32 {60    let blocks: Vec<Cycles> = pos61        .iter()62        .map(|&l| pos_block_cycles(l))63        .chain(neg.iter().map(|&l| neg_block_cycles(l)))64        .collect();65    if blocks.is_empty() {66        return 1;67    }68    let mut acc = blocks[0].clone();69    for b in &blocks[1..] {70        acc = combine(&acc, b);71    }72    acc.values().sum::<u128>() as u3273}7475fn partitions(n: usize, m: usize) -> Vec<Vec<usize>> {76    if n == 0 {77        return vec![vec![]];78    }79    let mut out = Vec::new();80    for k in (1..=n.min(m)).rev() {81        for rest in partitions(n - k, k) {82            let mut item = vec![k];83            item.extend(rest);84            out.push(item);85        }86    }87    out88}8990fn bipartitions(dimension: usize) -> Vec<(Vec<usize>, Vec<usize>)> {91    let mut out = Vec::new();92    for s in 0..=dimension {93        for pos in partitions(s, s.max(1)) {94            for neg in partitions(dimension - s, (dimension - s).max(1)) {95                out.push((pos.clone(), neg));96            }97        }98    }99    out100}101102fn class_size(pos: &[usize], neg: &[usize], dimension: usize) -> u128 {103    let mut centralizer: u128 = 1;104    for part in [pos, neg] {105        let mut mults = BTreeMap::new();106        for &length in part {107            *mults.entry(length).or_insert(0usize) += 1;108        }109        for (length, mult) in mults {110            centralizer *= factorial(mult) * (2 * length as u128).pow(mult as u32);111        }112    }113    ((1u128 << dimension) * factorial(dimension)) / centralizer114}115116/// Returns the raw design count of a dimension, two to the number of corners.117pub fn total_designs(dimension: usize) -> u128 {118    assert!(dimension <= 6, "total exceeds u128 above dimension 6");119    1 << (1 << dimension)120}121122/// Counts designs distinct under the full symmetry group class by class, or an error when the sums break.123pub fn distinct_designs(dimension: usize) -> Result<u128> {124    let order = (1u128 << dimension) * factorial(dimension);125    let mut total: u128 = 0;126    let mut checked: u128 = 0;127    for (pos, neg) in bipartitions(dimension) {128        let size = class_size(&pos, &neg, dimension);129        checked += size;130        total += size * (1u128 << class_cycles(&pos, &neg));131    }132    if checked != order {133        return value_error("class sizes do not sum to the group order.");134    }135    if !total.is_multiple_of(order) {136        return value_error("Burnside average is not an integer.");137    }138    Ok(total / order)139}140141/// Returns the distinct-design counts for dimensions 1 through max_dimension.142///143/// ```144/// assert_eq!(mrlymath::bang::counting::sequence(4).unwrap(), vec![3, 6, 22, 402]);145/// ```146pub fn sequence(max_dimension: usize) -> Result<Vec<u128>> {147    (1..=max_dimension).map(distinct_designs).collect()148}149150/// Counts the fill classes of a dimension, the popcount profiles a base-2 design can have: one more than the corners of each weight, multiplied over the weights, A129824 at the dimension.151pub fn classes(dimension: usize) -> u128 {152    let mut row = vec![1u128];153    for _ in 0..dimension {154        let mut next = vec![1u128; row.len() + 1];155        for k in 1..row.len() {156            next[k] = row[k - 1] + row[k];157        }158        row = next;159    }160    row.iter().map(|corners| corners + 1).product()161}162163/// Returns the fill-class counts for dimensions 1 through max_dimension.164///165/// ```166/// assert_eq!(mrlymath::bang::counting::class_sequence(4), vec![4, 12, 64, 700]);167/// ```168pub fn class_sequence(max_dimension: usize) -> Vec<u128> {169    (1..=max_dimension).map(classes).collect()170}171172#[cfg(test)]173mod tests {174    use super::*;175    #[test]176    fn burnside_matches_enumeration() {177        use super::super::universe::bang;178        for d in 1..=3 {179            assert_eq!(distinct_designs(d).unwrap(), bang(d).distinct() as u128);180        }181    }182    #[test]183    fn sequence_is_a000616() {184        assert_eq!(185            sequence(6).unwrap(),186            vec![3, 6, 22, 402, 1228158, 400507806843728]187        );188    }189    #[test]190    fn classes_are_a129824() {191        assert_eq!(classes(0), 2);192        assert_eq!(193            class_sequence(7),194            vec![4, 12, 64, 700, 17424, 1053696, 160579584]195        );196    }197    #[test]198    fn totals_doubly_exponential() {199        let totals: Vec<u128> = (1..=4).map(total_designs).collect();200        assert_eq!(totals, vec![4, 16, 256, 65536]);201    }202}