diagonal.rs

1use crate::core::error::{value_error, Result};2use crate::core::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 in 0..tile.size() {12        if tile.at(flat) == 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/// use mrlyrs::math::bang::Code;33/// let gasket = mrlyrs::math::bang::factory::create(Code::from(126u64), 2, 3, 2, 1).unwrap();34/// let counts = mrlyrs::math::counts::profile_of_tile(&gasket, 4).unwrap();35/// assert_eq!(counts[15..=30].iter().copied().collect::<Vec<u128>>(), vec![81u128; 16]);36/// ```37///38/// # Errors39///40/// Errors for a tile that is not a hypercube, a level below one, a side past the widest, or counts past a u128.41pub fn profile_of_tile(tile: &Tensor, level: u32) -> Result<Vec<u128>> {42    let dimension = tile.shape.len();43    let number = tile.shape.first().copied().unwrap_or(0);44    if tile.shape.iter().any(|&side| side != number) {45        return value_error("the profile wants a hypercube tile.");46    }47    if level < 1 {48        return value_error("level must be at least 1.");49    }50    let side = match number.checked_pow(level) {51        Some(side) if side <= WIDEST => side,52        _ => return value_error(format!("the side must stay at or below {WIDEST}.")),53    };54    let weights = histogram(tile);55    let span = dimension * (side - 1) + 1;56    let mut poly = vec![0u128; span];57    poly[0] = 1;58    let mut step = 1usize;59    for _ in 0..level {60        let mut next = vec![0u128; span];61        for (exponent, &count) in poly.iter().enumerate() {62            if count == 0 {63                continue;64            }65            for (&weight, &multiplicity) in &weights {66                let slot = exponent + step * weight;67                match count68                    .checked_mul(multiplicity)69                    .and_then(|added| next[slot].checked_add(added))70                {71                    Some(total) => next[slot] = total,72                    None => return value_error("the slice counts overflow a u128."),73                }74            }75        }76        poly = next;77        step *= number;78    }79    Ok(poly)80}8182#[cfg(test)]83mod tests {84    use super::*;85    use crate::math::bang::{factory, Code};86    #[test]87    fn the_profile_sums_to_the_fill_and_matches_a_rendered_count() {88        for code in [Code(7), Code(9), Code(11)] {89            let tile = factory::create(code, 3, 2, 2, 1).unwrap();90            let counts = profile_of_tile(&tile, 3).unwrap();91            let rendered = factory::create(code, 3, 2, 2, 3).unwrap();92            assert_eq!(counts.iter().sum::<u128>(), u128::from(rendered.sum()));93            let mut direct = vec![0u128; counts.len()];94            for (flat, &cell) in rendered.bytes().unwrap().iter().enumerate() {95                if cell != 0 {96                    direct[flat / 27 + flat % 27] += 1;97                }98            }99            assert_eq!(counts, direct, "code={code}");100        }101        let tile = factory::create(Code(23), 3, 3, 2, 1).unwrap();102        assert_eq!(profile_of_tile(&tile, 1).unwrap(), [1, 3, 3, 6, 3, 3, 1]);103        assert!(profile_of_tile(&tile, 0).is_err());104        assert!(profile_of_tile(&crate::math::atoms::ones_3d(2), 20).is_err());105    }106}