diagonal.rs

1use mrlycore::errors::{value_error, Result};2use mrlycore::tensor::Tensor;3use std::collections::BTreeMap;45/// The widest side a profile spans.6pub const WIDEST: usize = 1 << 18;78fn histogram(tile: &Tensor) -> BTreeMap<usize, u128> {9    let shape = &tile.shape;10    let mut out = BTreeMap::new();11    for (flat, &cell) in tile.bytes().iter().enumerate() {12        if cell == 0 {13            continue;14        }15        let mut rest = flat;16        let mut weight = 0;17        for &side in shape.iter().rev() {18            weight += rest % side;19            rest /= side;20        }21        *out.entry(weight).or_insert(0) += 1;22    }23    out24}2526/// Counts the filled cells of the tile's level-fold power on every diagonal plane `x_1 + ... + x_D = s`.27///28/// The count is the coefficient of the digit polynomial, the level-fold product of the tile's29/// weight sums, so no cell of the power is ever built; the tile must be a hypercube.30///31/// ```32/// let gasket = mrlymath::bang::factory::create(126, 2, 3, 2, 1).unwrap();33/// let counts = mrlymath::formulas::profile_of_tile(&gasket, 4).unwrap();34/// assert_eq!(counts[15..=30].iter().copied().collect::<Vec<u128>>(), vec![81u128; 16]);35/// ```36pub fn profile_of_tile(tile: &Tensor, level: u32) -> Result<Vec<u128>> {37    let dimension = tile.shape.len();38    let number = tile.shape.first().copied().unwrap_or(0);39    if tile.shape.iter().any(|&side| side != number) {40        return value_error("the profile wants a hypercube tile.");41    }42    if level < 1 {43        return value_error("level must be at least 1.");44    }45    let side = match number.checked_pow(level) {46        Some(side) if side <= WIDEST => side,47        _ => return value_error(format!("the side must stay at or below {WIDEST}.")),48    };49    let weights = histogram(tile);50    let span = dimension * (side - 1) + 1;51    let mut poly = vec![0u128; span];52    poly[0] = 1;53    let mut step = 1usize;54    for _ in 0..level {55        let mut next = vec![0u128; span];56        for (exponent, &count) in poly.iter().enumerate() {57            if count == 0 {58                continue;59            }60            for (&weight, &multiplicity) in &weights {61                let slot = exponent + step * weight;62                match count63                    .checked_mul(multiplicity)64                    .and_then(|added| next[slot].checked_add(added))65                {66                    Some(total) => next[slot] = total,67                    None => return value_error("the slice counts overflow a u128."),68                }69            }70        }71        poly = next;72        step *= number;73    }74    Ok(poly)75}7677#[cfg(test)]78mod tests {79    use super::*;80    use crate::bang::factory;81    #[test]82    fn the_profile_sums_to_the_fill_and_matches_a_rendered_count() {83        for code in [7u128, 9, 11] {84            let tile = factory::create(code, 3, 2, 2, 1).unwrap();85            let counts = profile_of_tile(&tile, 3).unwrap();86            let rendered = factory::create(code, 3, 2, 2, 3).unwrap();87            assert_eq!(counts.iter().sum::<u128>(), u128::from(rendered.sum()));88            let mut direct = vec![0u128; counts.len()];89            for (flat, &cell) in rendered.bytes().iter().enumerate() {90                if cell != 0 {91                    direct[flat / 27 + flat % 27] += 1;92                }93            }94            assert_eq!(counts, direct, "code={code}");95        }96        let tile = factory::create(23, 3, 3, 2, 1).unwrap();97        assert_eq!(profile_of_tile(&tile, 1).unwrap(), [1, 3, 3, 6, 3, 3, 1]);98        assert!(profile_of_tile(&tile, 0).is_err());99        assert!(profile_of_tile(&mrlycore::atoms::ones_3d(2), 20).is_err());100    }101}