memory.rs

1use mrlycore::errors::{value_error, Result};23/// The largest digit span a rule may read, so its code fits a `u64`.4pub const SPAN: usize = 6;56/// A rule on `k` consecutive digits of a design word.7///8/// The alphabet is the `2^D` digit vectors `d` of `{0, 1}^D`, each read as the corner integer `c = sum_i d[i] 2^i`, with `x` bit `0`, `y` bit `1` and `z` bit `2`.9/// A window is `k` digits `(c_1, ..., c_k)` read as `w = sum_j c_j 2^(D (k - j))`, the first digit most significant, and the code has bit `w` set exactly when that window is allowed.10/// The code therefore lives in `[0, 2^(2^(k D)))`, which forces `k D <= SPAN`.11/// A word `d_1 d_2 ... d_L` is read coarsest digit first and is accepted when every window of `k` consecutive digits is allowed; for `L < k` there is no window, so every word is accepted.12/// Width 1 is the memoryless design of the same code: the cells of `bang dim <dim>, code <code>` at that level.13#[derive(Clone, Copy, Debug, PartialEq, Eq)]14pub struct Rule {15    /// The dimension `D`, one to three.16    pub dimension: usize,17    /// The window width `k`, at least one.18    pub width: usize,19    /// The window code, bit `w` set when window `w` is allowed.20    pub code: u64,21}2223impl Rule {24    /// Builds a rule, or an error when the dimension, the width or the code is out of range.25    ///26    /// ```27    /// assert!(mrlynum::memory::Rule::new(1, 2, 7).is_ok());28    /// assert!(mrlynum::memory::Rule::new(2, 4, 0).is_err());29    /// ```30    pub fn new(dimension: usize, width: usize, code: u64) -> Result<Rule> {31        if !(1..=3).contains(&dimension) {32            return value_error(format!("dimension {dimension} is not one, two or three."));33        }34        if width < 1 {35            return value_error("width must be at least one.");36        }37        if width * dimension > SPAN {38            return value_error(format!(39                "width {width} at dimension {dimension} reads {} digit bits, over the span of {SPAN}.",40                width * dimension41            ));42        }43        let rule = Rule {44            dimension,45            width,46            code: 0,47        };48        let bound = rule.codes();49        if u128::from(code) >= bound {50            return value_error(format!(51                "code {code} is out of range for width {width} at dimension {dimension} (0..{}).",52                bound - 153            ));54        }55        Ok(Rule {56            dimension,57            width,58            code,59        })60    }6162    /// Returns the rule that allows every window.63    ///64    /// ```65    /// assert_eq!(mrlynum::memory::Rule::full(1, 2).code, 15);66    /// ```67    pub fn full(dimension: usize, width: usize) -> Rule {68        let rule = Rule::new(dimension, width, 0).expect("a zero code is always in range");69        Rule {70            code: (rule.codes() - 1) as u64,71            ..rule72        }73    }7475    /// Returns the letter count `2^D`, the digit vectors of the cube's corners.76    pub fn letters(&self) -> usize {77        1 << self.dimension78    }7980    /// Returns the window count `2^(k D)`.81    pub fn windows(&self) -> usize {82        1 << (self.width * self.dimension)83    }8485    /// Returns the count of rules of this shape, `2^(2^(k D))`.86    pub fn codes(&self) -> u128 {87        1u128 << self.windows()88    }8990    /// Returns the state count `2^((k - 1) D)`, the windows of one digit less that the transfer matrix runs on.91    pub fn states(&self) -> usize {92        1 << ((self.width - 1) * self.dimension)93    }9495    /// Returns whether the window is allowed, and false for any window out of range.96    ///97    /// ```98    /// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();99    /// assert!(golden.allowed(2));100    /// assert!(!golden.allowed(3));101    /// ```102    pub fn allowed(&self, window: usize) -> bool {103        window < self.windows() && (self.code >> window) & 1 == 1104    }105106    /// Returns whether a word, coarsest digit first, is accepted.107    ///108    /// ```109    /// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();110    /// assert!(golden.accepts(&[1, 0, 1]));111    /// assert!(!golden.accepts(&[1, 1, 0]));112    /// ```113    pub fn accepts(&self, word: &[usize]) -> bool {114        if word.len() < self.width {115            return word.iter().all(|&c| c < self.letters());116        }117        word.windows(self.width).all(|slice| {118            let mut w = 0;119            for &c in slice {120                w = (w << self.dimension) | c;121            }122            self.allowed(w)123        })124    }125126    /// Returns the letters that stand in at least one allowed window.127    ///128    /// ```129    /// assert_eq!(mrlynum::memory::Rule::new(1, 2, 7).unwrap().alphabet(), vec![0, 1]);130    /// assert!(mrlynum::memory::Rule::new(1, 2, 0).unwrap().alphabet().is_empty());131    /// ```132    pub fn alphabet(&self) -> Vec<usize> {133        let mut seen = vec![false; self.letters()];134        for w in 0..self.windows() {135            if !self.allowed(w) {136                continue;137            }138            for j in 0..self.width {139                seen[(w >> (self.dimension * j)) & (self.letters() - 1)] = true;140            }141        }142        (0..self.letters()).filter(|&c| seen[c]).collect()143    }144}145146// THE TRANSFER MATRIX147148/// Returns the transfer matrix on the `(k - 1)`-windows: entry `(s, t)` is one when the window that overlaps state `s` onto state `t` is allowed.149///150/// At `k = 1` the one state is the empty window and the single entry counts the allowed letters.151///152/// ```153/// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();154/// assert_eq!(mrlynum::memory::transfer(&golden), vec![vec![1, 1], vec![1, 0]]);155/// ```156pub fn transfer(rule: &Rule) -> Vec<Vec<u64>> {157    let states = rule.states();158    let letters = rule.letters();159    let mut matrix = vec![vec![0u64; states]; states];160    if rule.width == 1 {161        matrix[0][0] = (0..letters).filter(|&c| rule.allowed(c)).count() as u64;162        return matrix;163    }164    for (s, row) in matrix.iter_mut().enumerate() {165        for c in 0..letters {166            let window = (s << rule.dimension) | c;167            if !rule.allowed(window) {168                continue;169            }170            row[window & (states - 1)] = 1;171        }172    }173    matrix174}175176fn step(matrix: &[Vec<u64>], vector: &[u64]) -> Option<Vec<u64>> {177    let mut out = vec![0u64; vector.len()];178    for (s, row) in matrix.iter().enumerate() {179        for (t, &entry) in row.iter().enumerate() {180            if entry == 0 || vector[t] == 0 {181                continue;182            }183            out[s] = out[s].checked_add(entry.checked_mul(vector[t])?)?;184        }185    }186    Some(out)187}188189// THE COUNTS190191/// Returns `N_W(L)`, the count of accepted words, for `L = 1 ..= levels`, and stops early on the level whose count overruns a `u64`.192///193/// For `L < k` every word of that length is accepted, so the count is `2^(D L)`.194///195/// ```196/// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();197/// assert_eq!(mrlynum::memory::counts(&golden, 6), vec![2, 3, 5, 8, 13, 21]);198/// ```199pub fn counts(rule: &Rule, levels: usize) -> Vec<u64> {200    let matrix = transfer(rule);201    let mut vector = vec![1u64; rule.states()];202    let mut out = Vec::with_capacity(levels);203    for level in 1..=levels {204        if level + 1 < rule.width {205            match 1u64.checked_shl((rule.dimension * level) as u32) {206                Some(count) => out.push(count),207                None => break,208            }209            continue;210        }211        if level + 1 > rule.width {212            match step(&matrix, &vector) {213                Some(next) => vector = next,214                None => break,215            }216        }217        out.push(vector.iter().sum());218    }219    out220}221222/// Returns the accepted words of the level as cell indices of the `2^L` grid, `x` from bit `0` of every digit, `y` from bit `1`, `z` from bit `2`, coarsest digit first.223///224/// The index is row major, `z side^2 + y side + x` with `side = 2^L`, so at `k = 1` the set is the level-`L` fill of the design of the same code.225///226/// ```227/// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();228/// assert_eq!(mrlynum::memory::cells(&golden, 2), vec![0, 1, 2]);229/// ```230pub fn cells(rule: &Rule, level: usize) -> Vec<u64> {231    let side = 1u64 << level;232    let mut out = Vec::new();233    let mut word = vec![0usize; level];234    walk(rule, level, 0, &mut word, &mut |word| {235        let mut place = vec![0u64; rule.dimension];236        for (j, &c) in word.iter().enumerate() {237            let weight = 1u64 << (level - 1 - j);238            for (axis, seat) in place.iter_mut().enumerate() {239                *seat += ((c >> axis) & 1) as u64 * weight;240            }241        }242        out.push(place.iter().rev().fold(0, |run, &seat| run * side + seat));243    });244    out245}246247fn walk(248    rule: &Rule,249    level: usize,250    at: usize,251    word: &mut Vec<usize>,252    emit: &mut impl FnMut(&[usize]),253) {254    if at == level {255        emit(word);256        return;257    }258    for c in 0..rule.letters() {259        word[at] = c;260        if at + 1 >= rule.width {261            let start = at + 1 - rule.width;262            let mut w = 0;263            for &digit in &word[start..=at] {264                w = (w << rule.dimension) | digit;265            }266            if !rule.allowed(w) {267                continue;268            }269        }270        walk(rule, level, at + 1, word, emit);271    }272}273274// THE GROWTH275276/// The absolute `l^1` move of the normalised iterate that stops the power iteration, counted only when it holds over three consecutive sweeps.277pub const TOLERANCE: f64 = 1e-14;278279/// The sweep cap of the power iteration.280pub const SWEEPS: usize = 100_000;281282/// Returns the Perron root of the transfer matrix, the count's growth per level.283///284/// The digraph is split into its strongly connected components first; a component carrying no cycle contributes nothing, and a rule whose components all die past the window returns exactly `0`.285/// Each cyclic component is irreducible, so `I + B` is primitive and the power iteration on it converges geometrically; the iteration stops on a sustained Cauchy move of the normalised vector, never on a plateau of one scalar, and the Collatz-Wielandt ratios bracket the root.286/// The estimate is then made exact against the component's characteristic polynomial, taken by integer Faddeev-LeVerrier: an integer root is returned exactly, and otherwise the simple Perron root is bisected to the resolution of an `f64`.287/// The answer is the largest component root.288///289/// ```290/// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();291/// assert!((mrlynum::memory::perron(&golden) - 1.618_033_988_749_895).abs() < 1e-12);292/// assert_eq!(mrlynum::memory::perron(&mrlynum::memory::Rule::new(1, 2, 0).unwrap()), 0.0);293/// ```294pub fn perron(rule: &Rule) -> f64 {295    let matrix = transfer(rule);296    let reach = reachability(&matrix);297    let states = matrix.len();298    let mut seen = vec![false; states];299    let mut best = 0.0f64;300    for s in 0..states {301        if seen[s] || (reach[s] >> s) & 1 == 0 {302            continue;303        }304        let part: Vec<usize> = (0..states)305            .filter(|&t| t == s || ((reach[s] >> t) & 1 == 1 && (reach[t] >> s) & 1 == 1))306            .collect();307        for &t in &part {308            seen[t] = true;309        }310        let block: Vec<Vec<i128>> = part311            .iter()312            .map(|&i| part.iter().map(|&j| matrix[i][j] as i128).collect())313            .collect();314        let root = component_root(&block);315        if root > best {316            best = root;317        }318    }319    best320}321322fn reachability(matrix: &[Vec<u64>]) -> Vec<u64> {323    let states = matrix.len();324    let mut reach: Vec<u64> = matrix325        .iter()326        .map(|row| {327            row.iter()328                .enumerate()329                .filter(|(_, &entry)| entry > 0)330                .map(|(t, _)| 1u64 << t)331                .sum()332        })333        .collect();334    for _ in 0..states {335        let old = reach.clone();336        for s in 0..states {337            for t in 0..states {338                if (old[s] >> t) & 1 == 1 {339                    reach[s] |= old[t];340                }341            }342        }343        if reach == old {344            break;345        }346    }347    reach348}349350fn multiply(left: &[Vec<i128>], right: &[Vec<i128>]) -> Vec<Vec<i128>> {351    let n = left.len();352    let mut out = vec![vec![0i128; n]; n];353    for i in 0..n {354        for t in 0..n {355            if left[i][t] == 0 {356                continue;357            }358            for j in 0..n {359                out[i][j] += left[i][t] * right[t][j];360            }361        }362    }363    out364}365366fn characteristic(block: &[Vec<i128>]) -> Vec<i128> {367    let n = block.len();368    let mut poly = vec![0i128; n + 1];369    poly[n] = 1;370    let mut carry = vec![vec![0i128; n]; n];371    for step in 1..=n {372        let mut product = multiply(block, &carry);373        for (i, row) in product.iter_mut().enumerate() {374            row[i] += poly[n - step + 1];375        }376        carry = product;377        let next = multiply(block, &carry);378        let trace: i128 = (0..n).map(|i| next[i][i]).sum();379        poly[n - step] = -trace / step as i128;380    }381    poly382}383384fn value(poly: &[i128], x: f64) -> f64 {385    poly.iter().rev().fold(0.0f64, |run, &c| run * x + c as f64)386}387388fn vanishes(poly: &[i128], x: i128) -> bool {389    let mut run = 0i128;390    for &c in poly.iter().rev() {391        run = match run.checked_mul(x).and_then(|v| v.checked_add(c)) {392            Some(v) => v,393            None => return false,394        };395    }396    run == 0397}398399/// Returns the seed itself, unpolished, when no sign change of the characteristic polynomial is found within `0.25` of it, so a component whose bracket search fails reads the Collatz-Wielandt midpoint and never a wrong root.400fn component_root(block: &[Vec<i128>]) -> f64 {401    let n = block.len();402    let mut vector = vec![1.0f64 / n as f64; n];403    let mut settled = 0usize;404    for _ in 0..SWEEPS {405        let mut next = vec![0.0f64; n];406        for (s, row) in block.iter().enumerate() {407            let mut sum = vector[s];408            for (t, &entry) in row.iter().enumerate() {409                sum += entry as f64 * vector[t];410            }411            next[s] = sum;412        }413        let mass: f64 = next.iter().sum();414        for seat in next.iter_mut() {415            *seat /= mass;416        }417        let move_size: f64 = next418            .iter()419            .zip(vector.iter())420            .map(|(a, b)| (a - b).abs())421            .sum();422        vector = next;423        if move_size <= TOLERANCE {424            settled += 1;425            if settled >= 3 {426                break;427            }428        } else {429            settled = 0;430        }431    }432    let mut low = f64::INFINITY;433    let mut high = 0.0f64;434    for (s, row) in block.iter().enumerate() {435        let image: f64 = row436            .iter()437            .enumerate()438            .map(|(t, &entry)| entry as f64 * vector[t])439            .sum();440        let ratio = image / vector[s];441        low = low.min(ratio);442        high = high.max(ratio);443    }444    let seed = 0.5 * (low + high);445    let poly = characteristic(block);446    let rounded = seed.round();447    if (0.0..1e18).contains(&rounded) && vanishes(&poly, rounded as i128) {448        return rounded;449    }450    let mut under;451    let mut over;452    let mut delta = 1e-13f64;453    loop {454        under = seed - delta;455        over = seed + delta;456        if value(&poly, under) < 0.0 && value(&poly, over) > 0.0 {457            break;458        }459        delta *= 4.0;460        if delta > 0.25 {461            return seed;462        }463    }464    for _ in 0..200 {465        let middle = 0.5 * (under + over);466        if middle <= under || middle >= over {467            break;468        }469        if value(&poly, middle) < 0.0 {470            under = middle;471        } else {472            over = middle;473        }474    }475    0.5 * (under + over)476}477478/// Returns the growth exponent `log_2 rho`, the growth per digit of the accepted word count.479///480/// ```481/// let full = mrlynum::memory::Rule::full(2, 1);482/// assert!((mrlynum::memory::exponent(&full) - 2.0).abs() < 1e-12);483/// ```484pub fn exponent(rule: &Rule) -> f64 {485    let rho = perron(rule);486    if rho <= 0.0 {487        return f64::NEG_INFINITY;488    }489    rho.log2()490}491492/// Returns the count of allowed windows `card W`, the bits the code sets inside its window range.493///494/// ```495/// assert_eq!(mrlynum::memory::allowed_windows(&mrlynum::memory::Rule::new(1, 2, 7).unwrap()), 3);496/// ```497pub fn allowed_windows(rule: &Rule) -> usize {498    let windows = rule.windows();499    let mask = if windows >= 64 {500        u64::MAX501    } else {502        (1u64 << windows) - 1503    };504    (rule.code & mask).count_ones() as usize505}506507/// Returns the memory number `kappa(W) = log_2(card W) / k - log_2 rho`, the bits a digit spends on memory.508///509/// The window budget `log_2(card W) / k` is what a word of length `mk` could carry if its `m` disjoint windows were free of each other, and `log_2 rho` is what it really carries, so `kappa >= 0` and the gap is the interference the rule pays for.510/// It is zero on every product rule `W = F^k`, where `card W = card F^k` and `rho = card F`, so every `k = 1` rule reads zero; it is `f64::INFINITY` when a window is allowed and nothing survives past the window, and zero when no window is allowed.511///512/// ```513/// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();514/// assert!((mrlynum::memory::kappa(&golden) - 0.098_239_336_730).abs() < 1e-9);515/// assert_eq!(mrlynum::memory::kappa(&mrlynum::memory::Rule::full(2, 2)), 0.0);516/// ```517pub fn kappa(rule: &Rule) -> f64 {518    let windows = allowed_windows(rule);519    if windows == 0 {520        return 0.0;521    }522    let rho = perron(rule);523    if rho <= 0.0 {524        return f64::INFINITY;525    }526    (windows as f64).log2() / rule.width as f64 - rho.log2()527}528529#[cfg(test)]530mod tests {531    use super::*;532533    #[test]534    fn width_one_is_the_memoryless_design() {535        for code in [1u64, 7, 11, 13, 14] {536            let rule = Rule::new(2, 1, code).unwrap();537            let fills = code.count_ones() as u64;538            assert_eq!(539                counts(&rule, 4),540                vec![fills, fills.pow(2), fills.pow(3), fills.pow(4)]541            );542            assert_eq!(kappa(&rule), 0.0);543            let level = cells(&rule, 1);544            let want: Vec<u64> = (0..4).filter(|c| (code >> c) & 1 == 1).collect();545            assert_eq!(level, want);546        }547    }548549    #[test]550    fn the_golden_rule_counts_the_fibonacci_numbers() {551        let golden = Rule::new(1, 2, 7).unwrap();552        assert_eq!(counts(&golden, 8), vec![2, 3, 5, 8, 13, 21, 34, 55]);553        assert!((perron(&golden) - 1.618_033_988_749_895).abs() < 1e-12);554        assert_eq!(format!("{:.6}", kappa(&golden)), "0.098239");555    }556557    #[test]558    fn the_perron_root_is_the_spectral_radius_on_every_reducible_rule() {559        let plastic = 1.324_717_957_244_746;560        let golden = 1.618_033_988_749_895;561        let quartic = 1.380_277_569_097_614;562        for (code, want) in [563            (5u64, 1.0),564            (62, plastic),565            (91, quartic),566            (95, golden),567            (125, plastic),568            (190, plastic),569        ] {570            let rule = Rule::new(1, 3, code).unwrap();571            let root = perron(&rule);572            assert!(573                (root - want).abs() < 1e-12,574                "code {code} reads {root} against {want}"575            );576        }577    }578579    #[test]580    fn the_supergolden_rule_counts_the_narayana_cows() {581        let rule = Rule::new(1, 3, 23).unwrap();582        assert_eq!(counts(&rule, 8), vec![2, 4, 4, 6, 9, 13, 19, 28]);583        let terms = counts(&rule, 12);584        for l in 5..terms.len() {585            assert_eq!(terms[l], terms[l - 1] + terms[l - 3]);586        }587    }588589    #[test]590    fn the_full_rule_counts_every_word() {591        for dimension in 1..=3 {592            for width in 1..=SPAN / dimension {593                let rule = Rule::full(dimension, width);594                let want: Vec<u64> = (1..=5).map(|l| 1u64 << (dimension * l)).collect();595                assert_eq!(counts(&rule, 5), want, "d{dimension} k{width}");596                assert!((exponent(&rule) - dimension as f64).abs() < 1e-12);597            }598        }599    }600601    #[test]602    fn the_empty_rule_dies_past_its_window() {603        for width in 1..=3 {604            let rule = Rule::new(1, width, 0).unwrap();605            let want: Vec<u64> = (1..=5)606                .map(|l| if l < width { 1 << l } else { 0 })607                .collect();608            assert_eq!(counts(&rule, 5), want, "k{width}");609            assert_eq!(perron(&rule), 0.0);610            assert!(cells(&rule, 5).is_empty());611        }612    }613614    #[test]615    fn the_cells_are_the_accepted_words() {616        let rule = Rule::new(2, 2, 0xf1e2).unwrap();617        let level = 4;618        let side = 1u64 << level;619        let drawn = cells(&rule, level);620        assert_eq!(drawn.len() as u64, counts(&rule, level)[level - 1]);621        let mut seen = std::collections::BTreeSet::new();622        for index in &drawn {623            assert!(seen.insert(*index));624            assert!(*index < side * side);625        }626    }627628    #[test]629    fn a_rule_out_of_range_is_refused() {630        assert!(Rule::new(0, 1, 0).is_err());631        assert!(Rule::new(4, 1, 0).is_err());632        assert!(Rule::new(1, 0, 0).is_err());633        assert!(Rule::new(3, 3, 0).is_err());634        assert!(Rule::new(1, 2, 16).is_err());635    }636}