baseq.rs

1use super::factory::residue_corners;2use super::universe::{permutations, Code};3use mrlycore::errors::{value_error, Result};4use mrlynum::classics::factorial;5use std::collections::{BTreeSet, HashMap};67/// The most cells a code walk visits, so that the walk stays within `2^20` codes.8pub const WALK_LIMIT: usize = 20;910/// Returns the distinct rotation and reflection maps of a base-q axis.11pub fn axis_maps(base: usize) -> Vec<Vec<usize>> {12    let mut out: Vec<Vec<usize>> = Vec::new();13    for b in 0..base {14        let rot: Vec<usize> = (0..base).map(|r| (r + b) % base).collect();15        let ref_: Vec<usize> = (0..base).map(|r| (base + b - r) % base).collect();16        if !out.contains(&rot) {17            out.push(rot);18        }19        if !out.contains(&ref_) {20            out.push(ref_);21        }22    }23    out24}2526/// Returns the symmetry group order counted from the enumerated axis maps.27pub fn group_order(base: usize, dimension: usize) -> u128 {28    (axis_maps(base).len() as u128).pow(dimension as u32) * factorial(dimension)29}3031/// Returns the closed-form group order the axis-map count must match.32pub fn predicted_group_order(base: usize, dimension: usize) -> u128 {33    let per_axis = if base == 2 { 2u128 } else { 2 * base as u128 };34    per_axis.pow(dimension as u32) * factorial(dimension)35}3637fn choices(axis: &[Vec<usize>], dimension: usize) -> Vec<Vec<usize>> {38    let mut out = vec![vec![]];39    for _ in 0..dimension {40        let mut next = Vec::new();41        for prefix in &out {42            for (i, _) in axis.iter().enumerate() {43                let mut item = prefix.clone();44                item.push(i);45                next.push(item);46            }47        }48        out = next;49    }50    out51}5253fn cycles(54    perm: &[usize],55    choice: &[usize],56    axis: &[Vec<usize>],57    cells: &[Vec<u8>],58    index: &HashMap<Vec<u8>, usize>,59) -> u32 {60    let apply = |corner: &[u8]| -> Vec<u8> {61        (0..corner.len())62            .map(|i| axis[choice[i]][corner[perm[i]] as usize] as u8)63            .collect()64    };65    let mut seen = vec![false; cells.len()];66    let mut count = 0;67    for start in 0..cells.len() {68        if seen[start] {69            continue;70        }71        count += 1;72        let mut j = start;73        while !seen[j] {74            seen[j] = true;75            j = index[&apply(&cells[j])];76        }77    }78    count79}8081/// Counts base-q designs distinct under symmetry, or an error when the Burnside average breaks.82pub fn distinct_designs(base: usize, dimension: usize) -> Result<u128> {83    if base < 1 {84        return value_error("base must be at least 1.");85    }86    if dimension < 1 {87        return value_error("dimension must be at least 1.");88    }89    let cells = residue_corners(dimension, base);90    let index: HashMap<Vec<u8>, usize> = cells91        .iter()92        .enumerate()93        .map(|(i, c)| (c.clone(), i))94        .collect();95    let axis = axis_maps(base);96    let mut order: u128 = 0;97    let mut total: u128 = 0;98    for perm in permutations(dimension) {99        for choice in choices(&axis, dimension) {100            order += 1;101            total += 1u128 << cycles(&perm, &choice, &axis, &cells, &index);102        }103    }104    if !total.is_multiple_of(order) {105        return value_error("Burnside average is not an integer.");106    }107    Ok(total / order)108}109110/// Returns the symmetry group as cell maps, each sending the cell at index `i` to `element[i]`.111pub fn group(base: usize, dimension: usize) -> Vec<Vec<usize>> {112    let cells = residue_corners(dimension, base);113    let axis = axis_maps(base);114    let mut out = Vec::new();115    for perm in permutations(dimension) {116        for choice in choices(&axis, dimension) {117            out.push(118                cells119                    .iter()120                    .map(|cell| {121                        (0..dimension).fold(0, |acc, i| {122                            acc * base + axis[choice[i]][cell[perm[i]] as usize]123                        })124                    })125                    .collect(),126            );127        }128    }129    out130}131132/// Carries a code through one group element.133pub fn carry(element: &[usize], code: Code) -> Code {134    element135        .iter()136        .enumerate()137        .filter(|(index, _)| code >> index & 1 == 1)138        .map(|(_, &image)| 1u128 << image)139        .sum()140}141142/// Returns every code a design reaches under the group.143pub fn orbit(group: &[Vec<usize>], code: Code) -> BTreeSet<Code> {144    group.iter().map(|element| carry(element, code)).collect()145}146147/// Returns the least code of the design's orbit.148pub fn canonical(group: &[Vec<usize>], code: Code) -> Code {149    orbit(group, code)150        .into_iter()151        .next()152        .expect("the group is not empty")153}154155/// Walks every code of a base and dimension and returns each orbit's least code with the orbit's size, or an error past the walk limit.156///157/// ```158/// assert_eq!(mrlymath::bang::baseq::representatives(3, 1).unwrap().len(), 4);159/// ```160pub fn representatives(base: usize, dimension: usize) -> Result<Vec<(Code, usize)>> {161    let cells = base.pow(dimension as u32);162    if cells > WALK_LIMIT {163        return value_error(format!(164            "base {base} dimension {dimension} has {cells} cells, past the walk limit of {WALK_LIMIT}."165        ));166    }167    let group = group(base, dimension);168    let mut seen = vec![false; 1 << cells];169    let mut out = Vec::new();170    for code in 0..1u128 << cells {171        if seen[code as usize] {172            continue;173        }174        let orbit = orbit(&group, code);175        out.push((code, orbit.len()));176        for member in orbit {177            seen[member as usize] = true;178        }179    }180    Ok(out)181}182183/// Returns the raw design count before symmetry, two to the number of cells.184pub fn total_designs(base: usize, dimension: usize) -> u128 {185    let cells = base.pow(dimension as u32);186    assert!(cells < 128, "too many cells for a u128 count");187    1 << cells188}189190/// Returns the distinct-design counts for dimensions 1 through max_dimension.191pub fn sequence(base: usize, max_dimension: usize) -> Result<Vec<u128>> {192    (1..=max_dimension)193        .map(|d| distinct_designs(base, d))194        .collect()195}196197/// Returns the distinct one-dimensional design counts for bases 1 through max_base.198pub fn bracelets(max_base: usize) -> Result<Vec<u128>> {199    (1..=max_base).map(|q| distinct_designs(q, 1)).collect()200}201202/// Returns the filled-cell count of a binary design at a side number, folded from its filled corners.203pub fn fill_from_corners(filled: &[Vec<u8>], number: usize, dimension: usize) -> u128 {204    let even = number.div_ceil(2) as u128;205    let odd = (number / 2) as u128;206    filled207        .iter()208        .map(|corner| {209            let popcount = corner.iter().filter(|&&b| b != 0).count();210            even.pow((dimension - popcount) as u32) * odd.pow(popcount as u32)211        })212        .sum()213}214215/// Returns the collapsed fill count at an even side number, or an error at odd.216pub fn even_fill_is_balanced(number: usize, dimension: usize, popcount: u128) -> Result<u128> {217    if !number.is_multiple_of(2) {218        return value_error("the duality collapse holds only at even number.");219    }220    Ok(((number / 2) as u128).pow(dimension as u32) * popcount)221}222223#[cfg(test)]224mod tests {225    use super::*;226    #[test]227    fn group_orders_match_prediction() {228        for base in 2..=5 {229            for dimension in 1..=3 {230                assert_eq!(231                    group_order(base, dimension),232                    predicted_group_order(base, dimension)233                );234            }235        }236    }237    #[test]238    fn base2_matches_bang() {239        use super::super::universe::bang;240        for d in 1..=3 {241            assert_eq!(distinct_designs(2, d).unwrap(), bang(d).distinct() as u128);242        }243    }244    #[test]245    fn the_walk_agrees_with_burnside_and_the_base_two_universe() {246        use super::super::catalog::universe_codes;247        for (base, dimension) in [(3usize, 1usize), (3, 2), (4, 1), (4, 2), (5, 1)] {248            let walk = representatives(base, dimension).unwrap();249            assert_eq!(250                walk.len() as u128,251                distinct_designs(base, dimension).unwrap(),252                "q={base} d={dimension}"253            );254            let total: usize = walk.iter().map(|&(_, size)| size).sum();255            assert_eq!(total, 1 << base.pow(dimension as u32));256        }257        for dimension in 1..=4 {258            let codes: Vec<Code> = representatives(2, dimension)259                .unwrap()260                .into_iter()261                .map(|(code, _)| code)262                .collect();263            assert_eq!(codes, universe_codes(dimension));264        }265        assert_eq!(representatives(3, 2).unwrap().len(), 26);266        assert_eq!(representatives(4, 2).unwrap().len(), 805);267        assert!(representatives(3, 3).is_err());268    }269    #[test]270    fn bracelet_sequence_is_a000029() {271        assert_eq!(bracelets(8).unwrap(), vec![2, 3, 4, 6, 8, 13, 18, 30]);272    }273}