design.rs

1use num_bigint::BigUint;2use std::collections::HashSet;34const PERMS: [[usize; 3]; 6] = [5    [0, 1, 2],6    [0, 2, 1],7    [1, 0, 2],8    [1, 2, 0],9    [2, 0, 1],10    [2, 1, 0],11];1213type Matrix = [[i64; 3]; 3];1415const IDENTITY: Matrix = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];1617#[derive(Clone, Copy)]18pub struct Element {19    perm: [usize; 3],20    sign: [i64; 3],21    shift: [i64; 3],22}2324pub fn group(n: usize) -> Vec<Element> {25    let mut out = Vec::with_capacity(48 * n * n * n);26    for perm in PERMS {27        for bits in 0..8i64 {28            let sign = [29                1 - 2 * (bits >> 2 & 1),30                1 - 2 * (bits >> 1 & 1),31                1 - 2 * (bits & 1),32            ];33            for s0 in 0..n as i64 {34                for s1 in 0..n as i64 {35                    for s2 in 0..n as i64 {36                        out.push(Element {37                            perm,38                            sign,39                            shift: [s0, s1, s2],40                        });41                    }42                }43            }44        }45    }46    out47}4849fn image(n: usize, g: &Element, x: [i64; 3]) -> usize {50    let m = n as i64;51    let mut index = 0usize;52    for t in 0..3 {53        let y = (g.sign[t] * x[g.perm[t]] + g.shift[t]).rem_euclid(m);54        index = index * n + y as usize;55    }56    index57}5859fn images(n: usize, g: &Element, img: &mut [u32]) {60    let mut index = 0usize;61    for x0 in 0..n as i64 {62        for x1 in 0..n as i64 {63            for x2 in 0..n as i64 {64                img[index] = image(n, g, [x0, x1, x2]) as u32;65                index += 1;66            }67        }68    }69}7071fn walk(img: &[u32], seen: &mut [bool]) -> usize {72    seen.fill(false);73    let mut count = 0;74    for start in 0..img.len() {75        if seen[start] {76            continue;77        }78        count += 1;79        let mut j = start;80        while !seen[j] {81            seen[j] = true;82            j = img[j] as usize;83        }84    }85    count86}8788fn burnside(n: usize, histogram: &[u64]) -> BigUint {89    let mut total = BigUint::from(0u64);90    for (cycles, tally) in histogram.iter().enumerate() {91        if *tally != 0 {92            total += BigUint::from(*tally) << cycles;93        }94    }95    let order = BigUint::from(48 * (n as u64).pow(3));96    assert!((&total % &order) == BigUint::from(0u64));97    total / order98}99100pub fn by_cycles(n: usize) -> BigUint {101    let cells = n * n * n;102    let mut histogram = vec![0u64; cells + 1];103    let mut img = vec![0u32; cells];104    let mut seen = vec![false; cells];105    for g in group(n) {106        images(n, &g, &mut img);107        histogram[walk(&img, &mut seen)] += 1;108    }109    burnside(n, &histogram)110}111112fn matrix(g: &Element) -> Matrix {113    let mut rows: Matrix = [[0; 3]; 3];114    for t in 0..3 {115        rows[t][g.perm[t]] = g.sign[t];116    }117    rows118}119120fn matmul(a: &Matrix, b: &Matrix, m: i64) -> Matrix {121    let mut out: Matrix = [[0; 3]; 3];122    for i in 0..3 {123        for j in 0..3 {124            out[i][j] = (0..3).map(|k| a[i][k] * b[k][j]).sum::<i64>().rem_euclid(m);125        }126    }127    out128}129130fn affine(a: &Matrix, v: [i64; 3], t: [i64; 3], m: i64) -> [i64; 3] {131    let mut out = [0i64; 3];132    for i in 0..3 {133        out[i] = ((0..3).map(|k| a[i][k] * v[k]).sum::<i64>() + t[i]).rem_euclid(m);134    }135    out136}137138fn fixed_points(a: &Matrix, t: [i64; 3], m: i64) -> u64 {139    let mut count = 0u64;140    for x0 in 0..m {141        for x1 in 0..m {142            for x2 in 0..m {143                let x = [x0, x1, x2];144                if affine(a, x, t, m) == x {145                    count += 1;146                }147            }148        }149    }150    count151}152153fn cycles_by_powers(g: &Element, m: i64) -> usize {154    let base = matrix(g);155    let identity = matmul(&IDENTITY, &IDENTITY, m);156    let mut power = matmul(&base, &IDENTITY, m);157    let mut offset = g.shift;158    let mut total = 0u64;159    let mut k = 0u64;160    loop {161        k += 1;162        total += fixed_points(&power, offset, m);163        if power == identity && offset == [0, 0, 0] {164            break;165        }166        power = matmul(&base, &power, m);167        offset = affine(&base, offset, g.shift, m);168    }169    assert!(total % k == 0);170    (total / k) as usize171}172173pub fn by_affine(n: usize) -> BigUint {174    let cells = n * n * n;175    let mut histogram = vec![0u64; cells + 1];176    for g in group(n) {177        histogram[cycles_by_powers(&g, n as i64)] += 1;178    }179    burnside(n, &histogram)180}181182pub fn by_orbits(n: usize) -> u64 {183    let cells = n * n * n;184    let mut img = vec![0u32; cells];185    let mut maps: HashSet<Vec<u32>> = HashSet::new();186    for g in group(n) {187        images(n, &g, &mut img);188        maps.insert(img.clone());189    }190    let mut canon: HashSet<u64> = HashSet::new();191    for mask in 0u64..1u64 << cells {192        let mut best = mask;193        for map in &maps {194            let mut moved = 0u64;195            for (i, target) in map.iter().enumerate() {196                if mask >> i & 1 == 1 {197                    moved |= 1u64 << target;198                }199            }200            best = best.min(moved);201        }202        canon.insert(best);203    }204    canon.len() as u64205}