fills.rs

1use crate::quasi::{degree, fit, fraction, leading, text};2use crate::tables::write_csv;3use mrlymath::bang::baseq::{4    canonical, distinct_designs, group, group_order, orbit, representatives, WALK_LIMIT,5};6use mrlymath::bang::counting;7use mrlymath::bang::factory::{corners_to_code, levels_code, residue_corners};8use mrlymath::bang::universe;9use mrlymath::bang::Code;10use mrlymath::formulas::{fill, void};11use mrlymath::rules::render;12use std::collections::{BTreeMap, BTreeSet};13use std::path::Path;1415const SIDES: usize = 12;16const CASES: [(usize, usize); 4] = [(2, 2), (2, 3), (3, 2), (3, 3)];17const SHOWN: usize = 6;1819pub fn cell_index(cell: &[u8], base: usize) -> usize {20    cell.iter().fold(0, |acc, &r| acc * base + r as usize)21}2223fn named(base: usize, dimension: usize) -> Vec<(&'static str, Code)> {24    let cells = residue_corners(dimension, base);25    let select = |keep: &dyn Fn(&[u8]) -> bool| {26        let filled: Vec<Vec<u8>> = cells.iter().filter(|cell| keep(cell)).cloned().collect();27        corners_to_code(&filled, dimension, base)28    };29    let free = if dimension == 2 { 0 } else { dimension - 1 };30    vec![31        ("carpet", levels_code(dimension, base, &[0, 1])),32        (33            "net",34            select(&|cell| cell.iter().map(|&r| r as usize).sum::<usize>() + 1 >= dimension),35        ),36        ("void", select(&|cell| cell.iter().all(|&r| r == cell[0]))),37        (38            "tree",39            select(&|cell| (0..dimension).filter(|&a| a != free).all(|a| cell[a] == 0)),40        ),41    ]42}4344fn burnside(base: usize, dimension: usize) -> u128 {45    distinct_designs(base, dimension).expect("the Burnside average is an integer")46}4748pub fn orbits_report() {49    println!("base 2 cube group: order 2^D D!, designs up to symmetry three ways");50    for dimension in 1..=4 {51        let walk = representatives(2, dimension)52            .expect("the walk stays under the code limit")53            .len();54        let order = group_order(2, dimension);55        let burnside = burnside(2, dimension);56        let classes = counting::distinct_designs(dimension).expect("the class sums close");57        println!(58            "D {dimension}: order {order} orbit walk {walk} Burnside {burnside} class sum {classes}"59        );60    }61}6263pub struct Row {64    pub base: usize,65    pub dimension: usize,66    pub code: Code,67    pub label: String,68    pub orbit: usize,69    pub popcount: u32,70    pub gf2: Option<i32>,71    pub gfq: i64,72    pub mobius: i64,73    pub poly: String,74    pub degree: String,75    pub lead: String,76    pub fills: Vec<u128>,77    pub voids: Vec<u128>,78}7980fn rendered(code: Code, number: usize, dimension: usize, base: usize) -> u128 {81    let tile = render(82        |r| code >> cell_index(r, base) & 1 == 1,83        number,84        dimension,85        base,86    )87    .expect("the tile renders");88    u128::from(tile.sum())89}9091fn live_degree(cells: &[Vec<u8>], live: impl Fn(usize) -> bool) -> i64 {92    cells93        .iter()94        .enumerate()95        .filter(|(index, _)| live(*index))96        .map(|(_, cell)| cell.iter().map(|&r| i64::from(r)).sum::<i64>())97        .max()98        .unwrap_or(-1)99}100101fn mobius_degree(code: Code, cells: &[Vec<u8>], dimension: usize, base: usize) -> i64 {102    let mut coeff: Vec<i64> = (0..cells.len()).map(|i| (code >> i & 1) as i64).collect();103    for axis in 0..dimension {104        for value in (1..base as u8).rev() {105            for (index, cell) in cells.iter().enumerate() {106                if cell[axis] == value {107                    let mut lower = cell.clone();108                    lower[axis] -= 1;109                    coeff[index] -= coeff[cell_index(&lower, base)];110                }111            }112        }113    }114    live_degree(cells, |index| coeff[index] != 0)115}116117fn inverse_modular(value: i64, modulus: i64) -> i64 {118    (1..modulus)119        .find(|candidate| candidate * value.rem_euclid(modulus) % modulus == 1)120        .expect("the pivot is invertible in a prime field")121}122123fn inverse_vandermonde(q: usize) -> Vec<Vec<i64>> {124    let modulus = q as i64;125    let mut rows: Vec<Vec<i64>> = (0..q)126        .map(|i| {127            let mut row: Vec<i64> = (0..q)128                .map(|j| (i as i64).pow(j as u32).rem_euclid(modulus))129                .collect();130            row.extend((0..q).map(|c| i64::from(c == i)));131            row132        })133        .collect();134    for column in 0..q {135        let pivot = (column..q)136            .find(|&r| rows[r][column] != 0)137            .expect("the Vandermonde matrix is nonsingular over a prime field");138        rows.swap(column, pivot);139        let scale = inverse_modular(rows[column][column], modulus);140        for value in rows[column].iter_mut() {141            *value = *value * scale % modulus;142        }143        for r in 0..q {144            if r != column && rows[r][column] != 0 {145                let factor = rows[r][column];146                let pivot_row = rows[column].clone();147                for (value, above) in rows[r].iter_mut().zip(&pivot_row) {148                    *value = (*value - factor * above).rem_euclid(modulus);149                }150            }151        }152    }153    rows.into_iter().map(|row| row[q..].to_vec()).collect()154}155156fn gfq_degree(code: Code, cells: &[Vec<u8>], dimension: usize, q: usize) -> i64 {157    let inverse = inverse_vandermonde(q);158    let mut coeff: Vec<i64> = (0..cells.len()).map(|i| (code >> i & 1) as i64).collect();159    for axis in 0..dimension {160        let mut next = vec![0i64; cells.len()];161        for cell in cells.iter().filter(|cell| cell[axis] == 0) {162            let line: Vec<i64> = (0..q)163                .map(|t| {164                    let mut key = cell.clone();165                    key[axis] = t as u8;166                    coeff[cell_index(&key, q)]167                })168                .collect();169            for (exponent, row) in inverse.iter().enumerate() {170                let acc: i64 = row.iter().zip(&line).map(|(a, b)| a * b).sum();171                let mut key = cell.clone();172                key[axis] = exponent as u8;173                next[cell_index(&key, q)] = acc.rem_euclid(q as i64);174            }175        }176        coeff = next;177    }178    live_degree(cells, |index| coeff[index] != 0)179}180181fn row(base: usize, dimension: usize, code: Code, orbit: usize, label: String) -> Row {182    let cells = residue_corners(dimension, base);183    let count =184        |n: usize, level: u32| fill(code, n, dimension, level, base).expect("the fill counts");185    let fills: Vec<u128> = (1..=SIDES).map(|n| count(n, 1)).collect();186    let voids: Vec<u128> = (1..=SIDES)187        .map(|n| void(code, n, dimension, 1, base).expect("the void counts"))188        .collect();189    for n in 1..=SIDES {190        assert!(191            rendered(code, n, dimension, base) == fills[n - 1],192            "the closed form matches the rendered tile"193        );194    }195    for n in [base, base + 1, 2 * base] {196        for level in [2u32, 3] {197            assert!(198                count(n, level) == fills[n - 1].pow(level),199                "the level law holds"200            );201        }202    }203    let (polys, collapses) = fit(code, dimension, base);204    let top = polys.iter().filter_map(|poly| degree(poly)).max();205    let poly = if collapses {206        text(&polys[0])207    } else {208        (0..base)209            .map(|class| format!("n\u{2261}{class}(mod {base}): {}", text(&polys[class])))210            .collect::<Vec<String>>()211            .join(" ; ")212    };213    Row {214        base,215        dimension,216        code,217        label,218        orbit,219        popcount: code.count_ones(),220        gf2: (base == 2).then(|| universe::degree(code, dimension)),221        gfq: gfq_degree(code, &cells, dimension, base),222        mobius: mobius_degree(code, &cells, dimension, base),223        poly,224        degree: top.map_or(String::from("-oo"), |d| d.to_string()),225        lead: fraction(&leading(&polys[0])),226        fills,227        voids,228    }229}230231fn case(base: usize, dimension: usize) -> (Vec<Row>, bool) {232    let group = group(base, dimension);233    let full = base.pow(dimension as u32) <= WALK_LIMIT;234    if full {235        let labels: BTreeMap<Code, &str> = named(base, dimension)236            .into_iter()237            .map(|(name, code)| (canonical(&group, code), name))238            .collect();239        let rows = representatives(base, dimension)240            .expect("the walk stays under the code limit")241            .into_iter()242            .map(|(code, size)| {243                let label = labels.get(&code).copied().unwrap_or_default().to_string();244                row(base, dimension, code, size, label)245            })246            .collect();247        (rows, true)248    } else {249        let rows = named(base, dimension)250            .into_iter()251            .map(|(name, code)| {252                let representative = canonical(&group, code);253                let size = orbit(&group, representative).len();254                row(base, dimension, representative, size, name.to_string())255            })256            .collect();257        (rows, false)258    }259}260261fn header() -> Vec<String> {262    let mut out: Vec<String> = [263        "base",264        "dim",265        "rep_code",266        "label",267        "orbit_size",268        "popcount",269        "gf2_degree",270        "gfq_degree",271        "intmobius_degree",272        "fill_poly",273        "fill_deg",274        "fill_lead",275    ]276    .iter()277    .map(|name| (*name).to_string())278    .collect();279    out.extend((1..=SIDES).map(|n| format!("fill_n{n}")));280    out.extend((1..=SIDES).map(|n| format!("void_n{n}")));281    out282}283284fn record(row: &Row) -> Vec<String> {285    let mut out = vec![286        row.base.to_string(),287        row.dimension.to_string(),288        row.code.to_string(),289        row.label.clone(),290        row.orbit.to_string(),291        row.popcount.to_string(),292        row.gf2.map(|d| d.to_string()).unwrap_or_default(),293        row.gfq.to_string(),294        row.mobius.to_string(),295        row.poly.clone(),296        row.degree.clone(),297        row.lead.clone(),298    ];299    out.extend(row.fills.iter().map(u128::to_string));300    out.extend(row.voids.iter().map(u128::to_string));301    out302}303304pub fn report(path: &Path) {305    println!("fill census: one row per cube-group orbit, sides 1..{SIDES}, every row rendered cell by cell");306    let mut rows: Vec<Row> = Vec::new();307    for (base, dimension) in CASES {308        let (case_rows, full) = case(base, dimension);309        let expected = burnside(base, dimension);310        let scope = if full { "full walk" } else { "named only" };311        let polys: BTreeSet<&str> = case_rows.iter().map(|row| row.poly.as_str()).collect();312        println!(313            "base {base} D {dimension}: {} designs, Burnside {expected}, {} distinct fill polynomials, {scope}",314            case_rows.len(),315            polys.len()316        );317        for row in case_rows.iter().filter(|row| !row.label.is_empty()) {318            println!(319                "  {:<6} code {:>9} orbit {:>3} popcount {:>2} degree {} lead {:<5} fill n1..n{SHOWN} {:?}",320                row.label,321                row.code,322                row.orbit,323                row.popcount,324                row.degree,325                row.lead,326                &row.fills[..SHOWN]327            );328        }329        rows.extend(case_rows);330    }331    let labels = rows.iter().filter(|row| !row.label.is_empty()).count();332    println!(333        "{} designs, {} columns, {labels} labelled rows, written sequences.csv",334        rows.len(),335        header().len()336    );337    let records: Vec<Vec<String>> = rows.iter().map(record).collect();338    write_csv(path, &header(), &records);339}