press.rs

1use mrlycore::errors::{value_error, Result};2use mrlymath::bang::factory::{code_to_corners, MagicLayer};3use mrlymath::bang::universe::Code;45/// The largest corner count the tally press accepts, keeping its table a million rows.6pub const CORNERS: usize = 20;78fn corner_count(dimension: usize, base: usize) -> usize {9    let count = base.pow(dimension as u32);10    assert!(count < 128, "too many corners for a u128 code");11    count12}1314/// Returns the corner-usage mask of a number, one bit per digit vector its expansion uses.15///16/// A number is read in base `base` to the power of the dimension, so each digit is one17/// residue corner of the design cube, and zero uses exactly the zero corner.18///19/// ```20/// assert_eq!(mrlylab::press::usage(0, 2, 2), 1);21/// assert_eq!(mrlylab::press::usage(6, 2, 2), 0b0110);22/// ```23pub fn usage(number: u128, dimension: usize, base: usize) -> Code {24    let radix = corner_count(dimension, base) as u128;25    if number == 0 {26        return 1;27    }28    let mut out: Code = 0;29    let mut rest = number;30    while rest > 0 {31        out |= 1 << (rest % radix);32        rest /= radix;33    }34    out35}3637/// Returns whether every digit vector of the number lies in the design.38///39/// This is the scalar membership rule of the sequence press: the number's base40/// `base` to the dimension digits, read as residue corners, must all be filled.41/// At dimension one it is the classic restricted-digit set.42///43/// ```44/// let members: Vec<u128> = (0..30).filter(|&n| mrlylab::press::member(0b0111, n, 2, 2)).collect();45/// assert_eq!(members, vec![0, 1, 2, 4, 5, 6, 8, 9, 10, 16, 17, 18, 20, 21, 22, 24, 25, 26]);46/// ```47pub fn member(code: Code, number: u128, dimension: usize, base: usize) -> bool {48    usage(number, dimension, base) & !code == 049}5051/// Returns the count of distinct digit vectors the number uses.52pub fn distinct(number: u128, dimension: usize, base: usize) -> u32 {53    usage(number, dimension, base).count_ones()54}5556/// Returns the number of designs of the dimension and base that contain the number.57///58/// A design contains the number exactly when it fills every used corner, so the count59/// is two to the free corners, and the average membership over all designs is one over60/// two to the distinct-vector count.61///62/// ```63/// assert_eq!(mrlylab::press::containing(6, 2, 2), 4);64/// ```65pub fn containing(number: u128, dimension: usize, base: usize) -> u128 {66    1 << (corner_count(dimension, base) as u32 - distinct(number, dimension, base))67}6869/// Splits a number into its dimension coordinates, one base digit peeled per axis in parallel.70///71/// ```72/// assert_eq!(mrlylab::press::coordinates(6, 2, 2), vec![1, 2]);73/// ```74pub fn coordinates(number: u128, dimension: usize, base: usize) -> Vec<u128> {75    let radix = corner_count(dimension, base) as u128;76    let mut out = vec![0u128; dimension];77    let mut rest = number;78    let mut place: u128 = 1;79    while rest > 0 {80        let mut corner = rest % radix;81        for axis in (0..dimension).rev() {82            out[axis] += (corner % base as u128) * place;83            corner /= base as u128;84        }85        rest /= radix;86        place *= base as u128;87    }88    out89}9091/// Weaves dimension coordinates back into their single interleaved number.92///93/// Panics when the woven number passes a hundred and twenty-eight bits.94pub fn interleave(coords: &[u128], base: usize) -> u128 {95    assert!(96        !coords.is_empty(),97        "interleave needs at least one coordinate"98    );99    let dimension = coords.len();100    let radix = corner_count(dimension, base) as u128;101    let mut digits = Vec::new();102    let mut rest: Vec<u128> = coords.to_vec();103    while rest.iter().any(|&c| c > 0) {104        let mut corner: u128 = 0;105        for value in rest.iter_mut() {106            corner = corner * base as u128 + *value % base as u128;107            *value /= base as u128;108        }109        digits.push(corner);110    }111    let mut out: u128 = 0;112    for &corner in digits.iter().rev() {113        out = out114            .checked_mul(radix)115            .and_then(|v| v.checked_add(corner))116            .expect("the woven number passes a hundred and twenty-eight bits");117    }118    out119}120121/// Returns the first members of a design in ascending order.122///123/// Stops early where the next member would pass a hundred and twenty-eight bits.124///125/// ```126/// assert_eq!(mrlylab::press::members(0b10, 1, 2, 5), vec![1, 3, 7, 15, 31]);127/// ```128pub fn members(code: Code, dimension: usize, base: usize, count: usize) -> Vec<u128> {129    let radix = corner_count(dimension, base) as u128;130    let allowed: Vec<u128> = (0..radix).filter(|&i| (code >> i) & 1 == 1).collect();131    let mut out = Vec::with_capacity(count);132    if count == 0 || allowed.is_empty() {133        return out;134    }135    if allowed[0] == 0 {136        out.push(0);137    }138    let mut length = 1usize;139    while out.len() < count {140        let mut slots = vec![0usize; length];141        if allowed[0] == 0 {142            slots[0] = 1;143            if allowed.len() == 1 {144                return out;145            }146        }147        'level: loop {148            let mut value: u128 = 0;149            let mut fits = true;150            for &slot in &slots {151                match value152                    .checked_mul(radix)153                    .and_then(|v| v.checked_add(allowed[slot]))154                {155                    Some(next) => value = next,156                    None => {157                        fits = false;158                        break;159                    }160                }161            }162            if !fits {163                return out;164            }165            out.push(value);166            if out.len() == count {167                return out;168            }169            for place in (0..length).rev() {170                slots[place] += 1;171                if slots[place] < allowed.len() {172                    continue 'level;173                }174                slots[place] = usize::from(place == 0 && allowed[0] == 0);175            }176            break;177        }178        length += 1;179        if radix.checked_pow(length as u32 - 1).is_none() {180            return out;181        }182    }183    out184}185186/// Counts the members of a design below the limit.187///188/// ```189/// assert_eq!(mrlylab::press::count_below(0b0111, 2, 2, 27), 18);190/// ```191pub fn count_below(code: Code, dimension: usize, base: usize, limit: u128) -> u128 {192    let radix = corner_count(dimension, base) as u128;193    let allowed: Vec<u128> = (0..radix).filter(|&i| (code >> i) & 1 == 1).collect();194    if limit == 0 || allowed.is_empty() {195        return 0;196    }197    let mut digits = Vec::new();198    let mut rest = limit;199    while rest > 0 {200        digits.push(rest % radix);201        rest /= radix;202    }203    digits.reverse();204    let k = allowed.len() as u128;205    let lead = allowed.iter().filter(|&&v| v != 0).count() as u128;206    let mut total: u128 = if allowed[0] == 0 { 1 } else { 0 };207    let mut power: u128 = 1;208    for _ in 1..digits.len() {209        total += lead * power;210        power *= k;211    }212    for (place, &digit) in digits.iter().enumerate() {213        let below = allowed214            .iter()215            .filter(|&&v| v < digit && (place > 0 || v != 0))216            .count() as u128;217        let tail = digits.len() - place - 1;218        total += below * k.pow(tail as u32);219        if !allowed.contains(&digit) {220            break;221        }222    }223    total224}225226/// The tally press: one pass over the integers weighs every design of a universe at once.227///228/// Each added number lands its weight in the bucket of its corner-usage mask, and a229/// design's total is the sum over the submasks of its code, so a single sweep prices230/// a Mertens sum, a member count or a prime count for all two to the corners designs.231pub struct Press {232    /// The design dimension of the universe.233    pub dimension: usize,234    /// The numeral base of the universe.235    pub base: usize,236    corners: usize,237    tallies: Vec<i128>,238}239240impl Press {241    /// Builds an empty press over every design of the dimension and base.242    ///243    /// Panics past twenty corners, where the bucket table leaves a million rows.244    pub fn new(dimension: usize, base: usize) -> Press {245        let corners = corner_count(dimension, base);246        assert!(247            corners <= CORNERS,248            "the tally press holds at most twenty corners"249        );250        Press {251            dimension,252            base,253            corners,254            tallies: vec![0; 1 << corners],255        }256    }257    /// Adds a weighted number to its usage bucket.258    pub fn add(&mut self, number: u128, weight: i128) {259        self.tallies[usage(number, self.dimension, self.base) as usize] += weight;260    }261    /// Returns the total weight the design at a code has collected.262    pub fn total(&self, code: Code) -> i128 {263        let mut sum = self.tallies[0];264        let mut sub = code;265        while sub != 0 {266            sum += self.tallies[sub as usize];267            sub = (sub - 1) & code;268        }269        sum270    }271    /// Returns every design's total in code order by one subset-sum transform.272    pub fn totals(&self) -> Vec<i128> {273        let mut out = self.tallies.clone();274        for bit in 0..self.corners {275            for mask in 0..out.len() {276                if mask >> bit & 1 == 1 {277                    out[mask] += out[mask ^ (1 << bit)];278                }279            }280        }281        out282    }283}284285fn layer_radix(layer: &MagicLayer) -> u128 {286    (layer.number as u128).pow(layer.design.dim as u32)287}288289/// Returns the allowed digit table of one magic layer, one flag per cell of its tile.290///291/// A cell is allowed when its coordinate residues form a filled corner, which is the292/// tile the layer renders read as a digit alphabet.293pub fn layer_table(layer: &MagicLayer) -> Result<Vec<bool>> {294    let corners = code_to_corners(layer.design.code, layer.design.dim, layer.design.base)?;295    let dimension = layer.design.dim;296    let base = layer.design.base;297    let side = layer.number;298    let mut out = Vec::with_capacity(side.pow(dimension as u32));299    for cell in 0..side.pow(dimension as u32) {300        let mut rest = cell;301        let mut residue = vec![0u8; dimension];302        for axis in (0..dimension).rev() {303            residue[axis] = ((rest % side) % base) as u8;304            rest /= side;305        }306        out.push(corners.contains(&residue));307    }308    Ok(out)309}310311fn word_tables(layers: &[MagicLayer]) -> Result<Vec<Vec<bool>>> {312    if layers.is_empty() {313        return value_error("a word needs at least one layer.");314    }315    let dimension = layers[0].design.dim;316    if layers.iter().any(|l| l.design.dim != dimension) {317        return value_error("all word layers must have the same dimension.");318    }319    layers.iter().map(layer_table).collect()320}321322/// Counts the members of a magic word from its layer fills, without enumeration.323pub fn word_count(layers: &[MagicLayer]) -> Result<u128> {324    let tables = word_tables(layers)?;325    Ok(tables326        .iter()327        .map(|t| t.iter().filter(|&&b| b).count() as u128)328        .product())329}330331/// Returns whether the number lies in the magic word's composed design.332///333/// The number is read in the word's mixed radix, one digit per layer with the first334/// layer most significant, and every digit must land on an allowed cell of its tile.335/// A number past the word's domain is an error.336pub fn word_member(layers: &[MagicLayer], number: u128) -> Result<bool> {337    let tables = word_tables(layers)?;338    let mut rest = number;339    let mut ok = true;340    for (layer, table) in layers.iter().zip(&tables).rev() {341        let radix = layer_radix(layer);342        ok &= table[(rest % radix) as usize];343        rest /= radix;344    }345    if rest > 0 {346        return value_error(format!("number {number} lies past the word's domain."));347    }348    Ok(ok)349}350351/// Enumerates every member of the magic word in ascending order.352///353/// The member count is the product of the layer fills, so measure with `word_count`354/// before pressing a word too rich to hold.355pub fn word_members(layers: &[MagicLayer]) -> Result<Vec<u128>> {356    let tables = word_tables(layers)?;357    let alphabets: Vec<Vec<u128>> = tables358        .iter()359        .map(|t| {360            t.iter()361                .enumerate()362                .filter(|(_, &b)| b)363                .map(|(i, _)| i as u128)364                .collect()365        })366        .collect();367    if alphabets.iter().any(|a| a.is_empty()) {368        return Ok(Vec::new());369    }370    let radixes: Vec<u128> = layers.iter().map(layer_radix).collect();371    let mut out = Vec::new();372    let mut slots = vec![0usize; layers.len()];373    loop {374        let mut value: u128 = 0;375        for (place, &slot) in slots.iter().enumerate() {376            value = value * radixes[place] + alphabets[place][slot];377        }378        out.push(value);379        let mut place = layers.len();380        loop {381            if place == 0 {382                return Ok(out);383            }384            place -= 1;385            slots[place] += 1;386            if slots[place] < alphabets[place].len() {387                break;388            }389            slots[place] = 0;390        }391    }392}393394/// Returns the diagonal slice profile of a magic word by the substitution product.395///396/// The profile of a tile lists, per coordinate sum, its filled cells, and the profile397/// of a Kronecker word is the product of its layer profiles with strides, so no cell398/// of the composed design is ever enumerated.399pub fn word_profile(layers: &[MagicLayer]) -> Result<Vec<u128>> {400    let tables = word_tables(layers)?;401    let dimension = layers[0].design.dim;402    let mut out = vec![1u128];403    let mut stride: usize = 1;404    for (layer, table) in layers.iter().zip(&tables).rev() {405        let side = layer.number;406        let mut profile = vec![0u128; dimension * (side - 1) + 1];407        for (cell, &filled) in table.iter().enumerate() {408            if filled {409                let mut rest = cell;410                let mut total = 0usize;411                for _ in 0..dimension {412                    total += rest % side;413                    rest /= side;414                }415                profile[total] += 1;416            }417        }418        let mut next = vec![0u128; (profile.len() - 1) * stride + out.len()];419        for (i, &a) in profile.iter().enumerate() {420            if a == 0 {421                continue;422            }423            for (j, &b) in out.iter().enumerate() {424                next[i * stride + j] += a * b;425            }426        }427        out = next;428        stride *= side;429    }430    Ok(out)431}432433/// Returns the diagonal slice profile of one design pressed to a fractal level.434pub fn profile(code: Code, dimension: usize, base: usize, level: usize) -> Result<Vec<u128>> {435    if level < 1 {436        return value_error("level must be at least 1.");437    }438    let layer = MagicLayer::new(mrlymath::name::Bang::new(code, dimension, base), base);439    word_profile(&vec![layer; level])440}441442#[cfg(test)]443mod tests {444    use super::*;445    use mrlymath::bang::factory::create;446    use mrlymath::name::Bang;447448    #[test]449    fn usage_of_zero_is_the_zero_corner() {450        assert_eq!(usage(0, 2, 2), 1);451        assert_eq!(usage(0, 1, 10), 1);452        assert_eq!(distinct(0, 2, 2), 1);453    }454455    #[test]456    fn membership_matches_a_digit_check_at_base_ten() {457        let no_seven: Code = !(1 << 7) & ((1 << 10) - 1);458        for n in 0..10_000u128 {459            let digits_clean = !n.to_string().contains('7');460            assert_eq!(member(no_seven, n, 1, 10), digits_clean, "{n}");461        }462    }463464    #[test]465    fn members_of_the_repunit_design_are_the_mersenne_numbers() {466        assert_eq!(members(0b10, 1, 2, 6), vec![1, 3, 7, 15, 31, 63]);467    }468469    #[test]470    fn members_walk_ascending_and_agree_with_membership() {471        for code in [0b0111u128, 0b0110, 0b1001, 0b1111] {472            let list = members(code, 2, 2, 40);473            for pair in list.windows(2) {474                assert!(pair[0] < pair[1]);475            }476            let scanned: Vec<u128> = (0..200).filter(|&n| member(code, n, 2, 2)).collect();477            let shared = list.len().min(scanned.len());478            assert_eq!(list[..shared], scanned[..shared], "{code}");479        }480    }481482    #[test]483    fn count_below_agrees_with_the_member_walk() {484        for code in [0b0111u128, 0b0110, 0b1011, 0b0001, 0b0000] {485            let list = members(code, 2, 2, 60);486            for limit in 0..300u128 {487                let walked = list.iter().filter(|&&m| m < limit).count() as u128;488                if list.len() < 60 || walked < 60 {489                    assert_eq!(count_below(code, 2, 2, limit), walked, "{code} {limit}");490                }491            }492        }493    }494495    #[test]496    fn the_member_count_at_a_level_boundary_is_the_geometric_sum() {497        let code: Code = 0b0111;498        let k: u128 = 3;499        for level in 1..6u32 {500            let boundary = 4u128.pow(level);501            let full: u128 = 1 + (k - 1) * (k.pow(level) - 1) / (k - 1);502            assert_eq!(count_below(code, 2, 2, boundary), full);503        }504    }505506    #[test]507    fn membership_matches_the_rendered_fractal() {508        for code in [7u128, 6, 9, 11] {509            let level = 3;510            let tile = create(code, 2, 2, 2, level).unwrap();511            let side = 1u128 << level;512            for n in 0..4u128.pow(level as u32) {513                let coords = coordinates(n, 2, 2);514                let flat = (coords[0] * side + coords[1]) as usize;515                let filled = tile.bytes()[flat] == 1;516                let padded =517                    member(code, n, 2, 2) && (code & 1 == 1 || n >= 4u128.pow(level as u32 - 1));518                assert_eq!(padded, filled, "{code} {n}");519            }520        }521    }522523    #[test]524    fn coordinates_and_interleave_round_trip() {525        for n in 0..5_000u128 {526            assert_eq!(interleave(&coordinates(n, 2, 2), 2), n);527            assert_eq!(interleave(&coordinates(n, 3, 2), 2), n);528            assert_eq!(interleave(&coordinates(n, 2, 3), 3), n);529        }530    }531532    #[test]533    fn containing_counts_the_designs_that_hold_the_number() {534        for n in 0..500u128 {535            let direct = (0..16u128).filter(|&code| member(code, n, 2, 2)).count();536            assert_eq!(containing(n, 2, 2), direct as u128, "{n}");537        }538    }539540    #[test]541    fn the_membership_average_over_all_designs_is_two_to_minus_distinct() {542        for n in 0..2_000u128 {543            assert_eq!(containing(n, 2, 2), 1 << (4 - distinct(n, 2, 2)));544            assert_eq!(containing(n, 3, 2), 1 << (8 - distinct(n, 3, 2)));545        }546    }547548    #[test]549    fn the_press_totals_agree_with_direct_member_sums() {550        let mut press = Press::new(2, 2);551        let weights: Vec<i128> = (0..600).map(|n| (n as i128 % 7) - 3).collect();552        for (n, &w) in weights.iter().enumerate() {553            press.add(n as u128, w);554        }555        let totals = press.totals();556        for code in 0..16u128 {557            let direct: i128 = weights558                .iter()559                .enumerate()560                .filter(|(n, _)| member(code, *n as u128, 2, 2))561                .map(|(_, &w)| w)562                .sum();563            assert_eq!(press.total(code), direct, "{code}");564            assert_eq!(totals[code as usize], direct, "{code}");565        }566    }567568    #[test]569    fn a_native_word_is_the_stationary_press() {570        let layer = MagicLayer::new(Bang::new(7, 2, 2), 2);571        let word = vec![layer; 3];572        for n in 0..64u128 {573            let padded = member(7, n, 2, 2);574            assert_eq!(word_member(&word, n).unwrap(), padded, "{n}");575        }576        assert_eq!(word_count(&word).unwrap(), 27);577    }578579    #[test]580    fn word_members_match_the_magic_tensor() {581        let word = [582            MagicLayer::new(Bang::new(7, 2, 2), 3),583            MagicLayer::new(Bang::new(14, 2, 2), 5),584        ];585        let tensor = mrlymath::bang::factory::magic(&word).unwrap();586        let side = 15u128;587        let list = word_members(&word).unwrap();588        assert_eq!(list.len() as u128, word_count(&word).unwrap());589        for n in 0..word.iter().map(layer_radix).product::<u128>() {590            let mut rest = n;591            let mut x = 0u128;592            let mut y = 0u128;593            let mut place = 1u128;594            for layer in word.iter().rev() {595                let cell = rest % layer_radix(layer);596                rest /= layer_radix(layer);597                let s = layer.number as u128;598                x += (cell / s) * place;599                y += (cell % s) * place;600                place *= s;601            }602            let filled = tensor.bytes()[(x * side + y) as usize] == 1;603            assert_eq!(word_member(&word, n).unwrap(), filled, "{n}");604            assert_eq!(list.contains(&n), filled, "{n}");605        }606    }607608    #[test]609    fn word_profiles_match_the_rendered_diagonal_sums() {610        let word = [611            MagicLayer::new(Bang::new(7, 2, 2), 3),612            MagicLayer::new(Bang::new(14, 2, 2), 5),613        ];614        let tensor = mrlymath::bang::factory::magic(&word).unwrap();615        let side = 15usize;616        let mut direct = vec![0u128; 2 * side - 1];617        for (flat, &b) in tensor.bytes().iter().enumerate() {618            if b == 1 {619                direct[flat / side + flat % side] += 1;620            }621        }622        assert_eq!(word_profile(&word).unwrap(), direct);623    }624625    #[test]626    fn profile_totals_are_the_fill_powers() {627        let native = profile(7, 2, 2, 4).unwrap();628        assert_eq!(native.iter().sum::<u128>(), 3u128.pow(4));629        let sponge = vec![MagicLayer::new(Bang::new(23, 3, 2), 3); 3];630        let classic = word_profile(&sponge).unwrap();631        assert_eq!(classic.iter().sum::<u128>(), 20u128.pow(3));632    }633634    #[test]635    #[should_panic(expected = "the tally press holds at most twenty corners")]636    fn the_press_refuses_a_universe_past_twenty_corners() {637        let _ = Press::new(5, 2);638    }639640    #[test]641    fn a_word_refuses_mismatched_dimensions_and_numbers_past_its_domain() {642        let plane = MagicLayer::new(Bang::new(7, 2, 2), 3);643        let cube = MagicLayer::new(Bang::new(23, 3, 2), 3);644        assert!(word_member(&[plane.clone(), cube], 0).is_err());645        assert!(word_member(&[plane], 9).is_err());646        assert!(word_member(&[], 0).is_err());647    }648}