universe.rs

1use std::collections::BTreeSet;23/// The bitmask of filled corners that names a design.4pub type Code = u128;56/// Returns every permutation of 0..n in sorted order.7pub fn permutations(n: usize) -> Vec<Vec<usize>> {8    if n == 0 {9        return vec![vec![]];10    }11    let mut out = Vec::new();12    let mut items: Vec<usize> = (0..n).collect();13    heap(&mut items, n, &mut out);14    out.sort();15    out16}1718fn heap(items: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {19    if k == 1 {20        out.push(items.clone());21        return;22    }23    for i in 0..k {24        heap(items, k - 1, out);25        if k.is_multiple_of(2) {26            items.swap(i, k - 1);27        } else {28            items.swap(0, k - 1);29        }30    }31}3233/// Returns the binary corners of a dimension in code order.34pub fn corners(dimension: usize) -> Vec<Vec<u8>> {35    (0..1usize << dimension)36        .map(|i| {37            (0..dimension)38                .map(|j| ((i >> (dimension - 1 - j)) & 1) as u8)39                .collect()40        })41        .collect()42}4344/// Returns the bit position a binary corner occupies in a code.45pub fn corner_index(corner: &[u8]) -> usize {46    corner.iter().fold(0, |acc, &b| (acc << 1) | b as usize)47}4849/// Returns the full symmetry group as axis permutations paired with flip patterns.50pub fn symmetries(dimension: usize) -> Vec<(Vec<usize>, Vec<u8>)> {51    let mut out = Vec::new();52    for perm in permutations(dimension) {53        for f in 0..1usize << dimension {54            let flips: Vec<u8> = (0..dimension)55                .map(|j| ((f >> (dimension - 1 - j)) & 1) as u8)56                .collect();57            out.push((perm.clone(), flips));58        }59    }60    out61}6263/// Applies a symmetry element to a corner.64pub fn apply(element: &(Vec<usize>, Vec<u8>), corner: &[u8]) -> Vec<u8> {65    let (perm, flips) = element;66    (0..corner.len())67        .map(|i| corner[perm[i]] ^ flips[i])68        .collect()69}7071/// Returns every code a design reaches under the full symmetry group.72pub fn orbit(code: Code, dimension: usize) -> BTreeSet<Code> {73    let cells = corners(dimension);74    let group = symmetries(dimension);75    let mut out = BTreeSet::new();76    for g in &group {77        let mut image: Code = 0;78        for (i, cell) in cells.iter().enumerate() {79            if (code >> i) & 1 == 1 {80                image |= 1 << corner_index(&apply(g, cell));81            }82        }83        out.insert(image);84    }85    out86}8788/// Returns the algebraic normal form coefficients of a code, one per corner.89pub fn anf(code: Code, dimension: usize) -> Vec<u8> {90    let cells = corners(dimension);91    let mut coeff: Vec<u8> = (0..cells.len()).map(|i| ((code >> i) & 1) as u8).collect();92    for axis in 0..dimension {93        for (i, cell) in cells.iter().enumerate() {94            if cell[axis] == 1 {95                let mut lower = cell.clone();96                lower[axis] = 0;97                coeff[i] ^= coeff[corner_index(&lower)];98            }99        }100    }101    coeff102}103104/// Returns the algebraic degree of a code, or -1 for the zero design.105pub fn degree(code: Code, dimension: usize) -> i32 {106    let cells = corners(dimension);107    let coeff = anf(code, dimension);108    cells109        .iter()110        .enumerate()111        .filter(|(i, _)| coeff[*i] == 1)112        .map(|(_, c)| c.iter().map(|&b| b as i32).sum())113        .max()114        .unwrap_or(-1)115}116117/// Returns whether no two filled corners of a code sit at Hamming distance one.118///119/// Such a design buries no face at any side and any level, so its surface is six per cell.120///121/// ```122/// assert!(mrlymath::bang::universe::total_exposure(129, 3));123/// assert!(!mrlymath::bang::universe::total_exposure(23, 3));124/// ```125pub fn total_exposure(code: Code, dimension: usize) -> bool {126    let cells = corners(dimension);127    for (i, cell) in cells.iter().enumerate() {128        if (code >> i) & 1 == 0 {129            continue;130        }131        for axis in 0..dimension {132            let mut neighbor = cell.clone();133            neighbor[axis] ^= 1;134            if (code >> corner_index(&neighbor)) & 1 == 1 {135                return false;136            }137        }138    }139    true140}141142/// Returns whether a code fills the all-even corner, the rule that touches every grid corner at odd side.143///144/// ```145/// assert!(mrlymath::bang::universe::touches_every_corner(23, 3));146/// assert!(!mrlymath::bang::universe::touches_every_corner(232, 3));147/// ```148pub fn touches_every_corner(code: Code, dimension: usize) -> bool {149    let all_even: Vec<u8> = vec![0; dimension];150    (code >> corner_index(&all_even)) & 1 == 1151}152153/// Formats the algebraic normal form of a code as a sum of monomials.154pub fn anf_string(code: Code, dimension: usize) -> String {155    const NAMES: [char; 6] = ['x', 'y', 'z', 'w', 'v', 'u'];156    let cells = corners(dimension);157    let coeff = anf(code, dimension);158    let mut order: Vec<usize> = (0..cells.len()).collect();159    order.sort_by_key(|&i| {160        (161            cells[i].iter().map(|&b| b as usize).sum::<usize>(),162            cells[i].clone(),163        )164    });165    let mut terms = Vec::new();166    for i in order {167        if coeff[i] == 1 {168            let popcount: usize = cells[i].iter().map(|&b| b as usize).sum();169            if popcount == 0 {170                terms.push("1".to_string());171            } else {172                terms.push(173                    (0..dimension)174                        .filter(|&j| cells[i][j] == 1)175                        .map(|j| NAMES[j])176                        .collect(),177                );178            }179        }180    }181    if terms.is_empty() {182        "0".to_string()183    } else {184        terms.join("+")185    }186}187188/// A single design with its place in the orbit structure.189#[derive(Clone, Debug)]190pub struct Design {191    /// The design's code.192    pub i: Code,193    /// The design's dimension.194    pub dimension: usize,195    /// Whether this code is the smallest in its orbit.196    pub canonical: bool,197    /// The smallest code in the orbit.198    pub class_rep: Code,199    /// The number of codes in the orbit.200    pub orbit_size: usize,201}202203impl Design {204    /// Returns the design's name as a line of prose, `bang dim 2, code 7`.205    pub fn name(&self) -> String {206        crate::name::Named::to_mrly(&crate::name::Bang::new(self.i, self.dimension, 2))207    }208    /// Returns the design's filled corners in sorted order.209    pub fn rule(&self) -> Vec<Vec<u8>> {210        let cells = corners(self.dimension);211        let mut out: Vec<Vec<u8>> = cells212            .into_iter()213            .enumerate()214            .filter(|(i, _)| (self.i >> i) & 1 == 1)215            .map(|(_, c)| c)216            .collect();217        out.sort();218        out219    }220    /// Returns the design's algebraic degree, or -1 for the zero design.221    pub fn degree(&self) -> i32 {222        degree(self.i, self.dimension)223    }224    /// Returns the design's algebraic normal form as a string.225    pub fn anf(&self) -> String {226        anf_string(self.i, self.dimension)227    }228}229230/// The complete enumeration of one dimension's designs and orbits.231pub struct Universe {232    /// The universe's dimension.233    pub dimension: usize,234    /// The number of codes in the universe.235    pub total: usize,236    class_rep: Vec<Code>,237    orbit_size: Vec<usize>,238}239240impl Universe {241    /// Enumerates every orbit of a dimension from 1 to 4.242    pub fn new(dimension: usize) -> Self {243        assert!(244            (1..=4).contains(&dimension),245            "bang is enumerable only for dimensions 1-4"246        );247        let total = 1usize << (1usize << dimension);248        let mut class_rep = Vec::with_capacity(total);249        let mut orbit_size = Vec::with_capacity(total);250        for code in 0..total {251            let orb = orbit(code as Code, dimension);252            class_rep.push(*orb.iter().next().unwrap());253            orbit_size.push(orb.len());254        }255        Universe {256            dimension,257            total,258            class_rep,259            orbit_size,260        }261    }262    /// Returns the design at a code with its precomputed orbit facts.263    pub fn design(&self, code: Code) -> Design {264        let rep = self.class_rep[code as usize];265        Design {266            i: code,267            dimension: self.dimension,268            canonical: rep == code,269            class_rep: rep,270            orbit_size: self.orbit_size[code as usize],271        }272    }273    /// Returns every design in code order.274    pub fn all(&self) -> Vec<Design> {275        (0..self.total)276            .map(|code| self.design(code as Code))277            .collect()278    }279    /// Returns the designs whose codes lead their orbits.280    pub fn canonical(&self) -> Vec<Design> {281        self.all().into_iter().filter(|d| d.canonical).collect()282    }283    /// Returns the number of distinct orbits.284    pub fn distinct(&self) -> usize {285        let mut reps: Vec<Code> = self.class_rep.clone();286        reps.sort();287        reps.dedup();288        reps.len()289    }290}291292/// Builds the universe of a dimension.293///294/// ```295/// let u = mrlymath::bang::bang(2);296/// assert_eq!(u.total, 16);297/// assert_eq!(u.distinct(), 6);298/// ```299pub fn bang(dimension: usize) -> Universe {300    Universe::new(dimension)301}302303#[cfg(test)]304mod tests {305    use super::*;306    #[test]307    fn total_and_distinct_counts() {308        assert_eq!(bang(1).distinct(), 3);309        assert_eq!(bang(2).distinct(), 6);310        assert_eq!(bang(3).distinct(), 22);311        assert_eq!(bang(1).total, 4);312        assert_eq!(bang(2).total, 16);313        assert_eq!(bang(3).total, 256);314    }315    #[test]316    fn prefix_codes_canonical() {317        for d in 1..=3 {318            let u = bang(d);319            for k in 0..=(1usize << d) {320                let code = (1u128 << k) - 1;321                if (code as usize) < u.total {322                    assert!(u.design(code).canonical);323                }324            }325        }326    }327    #[test]328    fn anti_closure_3d() {329        let u = bang(3);330        let full: Code = (1 << (1 << 3)) - 1;331        let reps: Vec<Code> = u.canonical().iter().map(|d| d.class_rep).collect();332        for d in u.canonical() {333            let anti = full ^ d.i;334            assert!(reps.contains(&u.design(anti).class_rep));335        }336    }337    #[test]338    fn orbit_sizes_partition() {339        for d in 2..=3usize {340            let u = bang(d);341            let order = (1usize << d) * (1..=d).product::<usize>();342            let total: usize = u.canonical().iter().map(|x| x.orbit_size).sum();343            assert_eq!(total, u.total);344            for x in u.canonical() {345                assert_eq!(order % x.orbit_size, 0);346            }347        }348    }349    #[test]350    fn degree_histogram_3d() {351        let u = bang(3);352        let mut hist = std::collections::HashMap::new();353        for d in u.canonical() {354            *hist.entry(d.degree()).or_insert(0) += 1;355        }356        let expected: std::collections::HashMap<i32, i32> =357            [(-1, 1), (0, 1), (1, 3), (2, 9), (3, 8)]358                .into_iter()359                .collect();360        assert_eq!(hist, expected);361    }362    #[test]363    fn total_exposure_names_the_independent_corner_sets() {364        let exposed: Vec<Code> = (0..256).filter(|&c| total_exposure(c, 3)).collect();365        assert_eq!(exposed.len(), 35);366        let classes: BTreeSet<Code> = exposed367            .iter()368            .map(|&c| *orbit(c, 3).iter().next().unwrap())369            .collect();370        assert_eq!(371            classes.into_iter().collect::<Vec<Code>>(),372            [0, 1, 6, 22, 24, 105]373        );374        assert!(total_exposure(129, 3));375        assert!(!total_exposure(23, 3));376    }377378    #[test]379    fn half_the_rules_hold_the_all_even_corner() {380        assert_eq!(381            (0..256).filter(|&c| touches_every_corner(c, 3)).count(),382            128383        );384        for code in [23u128, 3, 129] {385            assert!(touches_every_corner(code, 3), "code={code}");386        }387        assert!(!touches_every_corner(232, 3));388    }389390    #[test]391    fn names_and_anf() {392        let u = bang(2);393        assert_eq!(u.design(0).name(), "bang dim 2, code 0");394        assert_eq!(u.design(7).name(), "bang dim 2, code 7");395        assert_eq!(u.design(0).anf(), "0");396        assert_eq!(u.design(1).anf(), "1+y+x+xy");397    }398}