count.rs

1use std::collections::BTreeMap;23pub struct Dsu {4    parent: Vec<u32>,5}67impl Dsu {8    pub fn new(size: usize) -> Dsu {9        Dsu {10            parent: (0..size as u32).collect(),11        }12    }1314    pub fn find(&mut self, mut i: u32) -> u32 {15        while self.parent[i as usize] != i {16            let up = self.parent[self.parent[i as usize] as usize];17            self.parent[i as usize] = up;18            i = up;19        }20        i21    }2223    pub fn union(&mut self, a: u32, b: u32) -> bool {24        let (ra, rb) = (self.find(a), self.find(b));25        if ra == rb {26            return false;27        }28        self.parent[ra as usize] = rb;29        true30    }31}3233pub fn arcs(side: usize, on: &[bool]) -> (u64, u64) {34    let rows = side * (side + 1);35    let h = |x: usize, y: usize| (y * side + x) as u32;36    let v = |x: usize, y: usize| (rows + y * (side + 1) + x) as u32;37    let mut dsu = Dsu::new(2 * rows);38    let mut degree = vec![0u8; 2 * rows];39    for y in 0..side {40        for x in 0..side {41            let (bottom, top, left, right) = (h(x, y), h(x, y + 1), v(x, y), v(x + 1, y));42            let pairs = if on[y * side + x] {43                [(left, bottom), (right, top)]44            } else {45                [(bottom, right), (left, top)]46            };47            for (a, b) in pairs {48                degree[a as usize] += 1;49                degree[b as usize] += 1;50                dsu.union(a, b);51            }52        }53    }54    let mut closed = vec![1u8; 2 * rows];55    let mut root = vec![false; 2 * rows];56    for i in 0..2 * rows {57        let r = dsu.find(i as u32) as usize;58        root[r] = true;59        if degree[i] != 2 {60            closed[r] = 0;61        }62    }63    let curves = root.iter().filter(|&&r| r).count() as u64;64    let loops = (0..2 * rows).filter(|&i| root[i] && closed[i] == 1).count() as u64;65    (loops, curves - loops)66}6768pub fn mirrors(side: usize, on: &[bool]) -> i64 {69    let w = side + 1;70    let p = |x: usize, y: usize| (y * w + x) as u32;71    let mut dsu = Dsu::new(w * w);72    let mut comp = (w * w) as i64;73    for y in 0..side {74        for x in 0..side {75            let (a, b) = if on[y * side + x] {76                (p(x + 1, y), p(x, y + 1))77            } else {78                (p(x, y), p(x + 1, y + 1))79            };80            if dsu.union(a, b) {81                comp -= 1;82            }83        }84    }85    comp - 2 * side as i64 - 186}8788pub struct Block {89    pub side: usize,90    pub partner: Vec<u32>,91    pub loops: u128,92    pub lengths: BTreeMap<u64, u64>,93}9495pub fn unit(on: bool) -> Block {96    let partner = if on {97        vec![3, 2, 1, 0]98    } else {99        vec![1, 0, 3, 2]100    };101    Block {102        side: 1,103        partner,104        loops: 0,105        lengths: BTreeMap::new(),106    }107}108109pub fn void_partner(s: usize, p: usize) -> usize {110    let t = p % s;111    match p / s {112        0 => 2 * s - 1 - t,113        1 => s - 1 - t,114        2 => 4 * s - 1 - t,115        _ => 3 * s - 1 - t,116    }117}118119pub fn void_block(s: usize) -> Block {120    Block {121        side: s,122        partner: (0..4 * s).map(|p| void_partner(s, p) as u32).collect(),123        loops: 0,124        lengths: BTreeMap::new(),125    }126}127128pub fn glue(n: usize, tile: &[bool], kept: &Block, keep: bool) -> Block {129    let s = kept.side;130    let big = n * s;131    let s4 = 4 * s;132    let mut seen = vec![0u64; (n * n * s4).div_ceil(64)];133    let mut partner = if keep {134        vec![u32::MAX; 4 * big]135    } else {136        Vec::new()137    };138    let inside = |b: usize, p: usize| {139        if tile[b] {140            kept.partner[p] as usize141        } else {142            void_partner(s, p)143        }144    };145    let outer = |i: usize, j: usize, p: usize| -> Option<usize> {146        let t = p % s;147        match p / s {148            0 => (i == 0).then_some(j * s + t),149            1 => (j == n - 1).then_some(big + i * s + t),150            2 => (i == n - 1).then_some(2 * big + j * s + t),151            _ => (j == 0).then_some(3 * big + i * s + t),152        }153    };154    let across = |i: usize, j: usize, p: usize| -> (usize, usize, usize) {155        let t = p % s;156        match p / s {157            0 => (i - 1, j, 2 * s + t),158            1 => (i, j + 1, 3 * s + t),159            2 => (i + 1, j, t),160            _ => (i, j - 1, s + t),161        }162    };163    let key = |i: usize, j: usize, p: usize| (i * n + j) * s4 + p;164    let mark = |seen: &mut Vec<u64>, k: usize| seen[k / 64] |= 1 << (k % 64);165    let is_seen = |seen: &Vec<u64>, k: usize| seen[k / 64] >> (k % 64) & 1 == 1;166    for i in 0..n {167        for j in 0..n {168            for p in 0..s4 {169                let Some(start) = outer(i, j, p) else {170                    continue;171                };172                if is_seen(&seen, key(i, j, p)) {173                    continue;174                }175                let (mut a, mut b, mut c) = (i, j, p);176                loop {177                    mark(&mut seen, key(a, b, c));178                    let q = inside(a * n + b, c);179                    mark(&mut seen, key(a, b, q));180                    if let Some(end) = outer(a, b, q) {181                        if keep {182                            partner[start] = end as u32;183                            partner[end] = start as u32;184                        }185                        break;186                    }187                    (a, b, c) = across(a, b, q);188                }189            }190        }191    }192    let mut cycles = 0u128;193    let mut lengths = BTreeMap::new();194    for i in 0..n {195        for j in 0..n {196            for p in 0..s4 {197                if is_seen(&seen, key(i, j, p)) {198                    continue;199                }200                cycles += 1;201                let mut length = 0u64;202                let (mut a, mut b, mut c) = (i, j, p);203                loop {204                    length += 1;205                    mark(&mut seen, key(a, b, c));206                    let q = inside(a * n + b, c);207                    mark(&mut seen, key(a, b, q));208                    (a, b, c) = across(a, b, q);209                    if (a, b, c) == (i, j, p) {210                        break;211                    }212                }213                *lengths.entry(length).or_insert(0) += 1;214            }215        }216    }217    let inner = tile.iter().filter(|&&on| on).count() as u128 * kept.loops;218    Block {219        side: big,220        partner,221        loops: inner + cycles,222        lengths,223    }224}225226pub fn series(n: usize, tile: &[bool], top: usize) -> Vec<u128> {227    let mut kept = unit(true);228    let mut out = vec![0];229    for level in 0..top {230        kept = glue(n, tile, &kept, level + 1 < top);231        out.push(kept.loops);232    }233    out234}235236pub fn blocks(n: usize, tile: &[bool], top: usize) -> Vec<Block> {237    let mut out = vec![unit(true)];238    for level in 0..top {239        let next = glue(n, tile, &out[level], true);240        out.push(next);241    }242    out243}244245#[cfg(test)]246mod tests {247    use super::*;248    use crate::design::Design;249250    #[test]251    fn three_counts_agree_at_base_two() {252        for code in 0..16u128 {253            let d = Design::full(code, 2);254            let blocks = series(2, &d.tile, 5);255            for level in 0..=5 {256                let (side, on) = d.cells(level);257                let (loops, strands) = arcs(side, &on);258                assert_eq!(loops as u128, blocks[level]);259                assert_eq!(mirrors(side, &on), loops as i64);260                assert_eq!(strands, 2 * side as u64);261            }262        }263    }264265    #[test]266    fn carpet_law() {267        let s = series(3, &Design::new(7, 3, 2).tile, 8);268        for (n, &x) in s.iter().enumerate() {269            let n = n as u32;270            assert_eq!(271                x as i128,272                (8i128.pow(n) - 1) / 7 - 3i128.pow(n) + n as i128 + 1273            );274        }275    }276277    #[test]278    fn void_block_glues_to_its_formula() {279        for n in 2..6 {280            let glued = glue(n, &vec![true; n * n], &void_block(n), true);281            assert_eq!(glued.partner, void_block(n * n).partner);282            assert_eq!(glued.loops, 0);283        }284    }285}