automaton.rs

1use crate::memory::{transfer, Rule};2use crate::zeta::{raise, Complex};3use mrlycore::errors::{value_error, Result};45/// The relative rounding allowance the double-precision matrix ladder charges against the scale it carries.6///7/// The truncation bound is the proved one, carried entrywise by the same recursion that carries the value vector; this constant is a measured allowance for the arithmetic itself and is not a proof.8pub const ROUNDING: f64 = 1e-13;910const PEEL_TARGET: f64 = 100.0;11const PEEL_CAP: usize = 26;12const SHIFT_START: f64 = 6.0;13const SHIFT_STEP: f64 = 6.0;14const SHIFT_CAP: f64 = 400.0;15const CUT_START: usize = 8;16const CUT_STEP: usize = 6;17const CUT_CAP: usize = 120;18const RATIO_CAP: f64 = 0.9;19const COLLATZ_WIELANDT: usize = 60;20const PIVOT_FLOOR: f64 = 1e-12;2122// THE STATE SPACE2324/// A memory design read as a matrix ladder: the rule, the transfer matrix on its `(k-1)`-window states, and the peel depth its Dirichlet series is continued from.25///26/// The elements are the positive integers whose minimal base-`q` string, `q = 2^D`, is accepted by the rule, and the object is `zeta_W(s) = sum n^(-s)` over them.27/// [`crate::memory::Rule::accepts`] holds a word shorter than the width `k` to be accepted, there being no window in it, so every integer below `q^(k-1)` with a nonzero leading digit sits in `S_W` for every rule of that width, and the head polynomial carries them.28/// With `E_j(w)` the vector whose entry `u` sums `n^(-w)` over accepted words of exactly `j` digits ending in state `u`, and `G_P = sum_(j >= P) E_j`, splitting a word on its last digit gives `(I - q^(-w) T) G_P(w) = E_P(w) + sum_(l >= 1) binom(-w, l) q^(-w-l) Gamma_l G_P(w+l)`, where `Gamma_l(u', u) = sum a^l` over the letters `a` with the window `(u, a)` allowed and `u' = shift(u, a)`, and `T = Gamma_0` is [`crate::memory::transfer`] transposed.29/// Then `zeta_W(w) = 1^T D_(P-1)(w) + 1^T G_P(w)` with `D_(P-1)` the Dirichlet polynomial over the accepted words of at most `P-1` digits, and the scalar ladder of [`crate::ladder`] is the width one case with `card F` where `T` stands.30#[derive(Clone, Debug, PartialEq)]31pub struct Automaton {32    rule: Rule,33    base: u64,34    states: usize,35    peel: usize,36    largest: u64,37    arcs: Vec<(usize, usize, u64)>,38    power: Vec<Vec<f64>>,39    low: Vec<f64>,40    mid: Vec<Vec<f64>>,41    seats: Vec<f64>,42    guide: Vec<f64>,43    perron_high: f64,44    perron_low: f64,45    charpoly: Vec<f64>,46    adjugate: Vec<Vec<Vec<f64>>>,47}4849fn walk_words(50    rule: &Rule,51    base: u64,52    depth: usize,53    word: &mut Vec<usize>,54    visit: &mut impl FnMut(usize, usize, u64),55) {56    if word.len() >= depth {57        return;58    }59    let start = usize::from(word.is_empty());60    for c in start..rule.letters() {61        word.push(c);62        if rule.accepts(word) {63            let value = word.iter().fold(0u64, |acc, &d| acc * base + d as u64);64            let span = word.len() - (rule.width - 1).min(word.len());65            let state = word[span..]66                .iter()67                .fold(0usize, |acc, &d| acc * base as usize + d);68            visit(word.len(), state, value);69            walk_words(rule, base, depth, word, visit);70        }71        word.pop();72    }73}7475impl Automaton {76    /// Builds the ladder of a rule, choosing the peel depth, or an error when the rule admits no element or its state space overruns the exact integers.77    ///78    /// ```79    /// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();80    /// let ladder = mrlynum::automaton::Automaton::new(&golden).unwrap();81    /// assert_eq!(ladder.states(), 2);82    /// assert!((ladder.abscissa() - 0.694_241_913_630_617_4).abs() < 1e-15);83    /// ```84    pub fn new(rule: &Rule) -> Result<Automaton> {85        let mut peel = rule.width.max(2);86        let base = rule.letters() as u64;87        while peel < PEEL_CAP && (base as f64).powi(peel as i32 + 1) < 9.0e15 {88            let mut seen = 0f64;89            let mut word = Vec::new();90            walk_words(rule, base, peel, &mut word, &mut |len, _, _| {91                if len == peel {92                    seen += 1.0;93                }94            });95            if seen > PEEL_TARGET {96                break;97            }98            peel += 1;99        }100        Automaton::with_peel(rule, peel)101    }102103    /// Builds the ladder at an explicit peel depth, at least the rule width and at least two, or an error when the rule admits no element or the depth overruns the exact integers.104    ///105    /// ```106    /// let full = mrlynum::memory::Rule::full(1, 2);107    /// let ladder = mrlynum::automaton::Automaton::with_peel(&full, 7).unwrap();108    /// assert_eq!(ladder.peel(), 7);109    /// assert_eq!(ladder.abscissa(), 1.0);110    /// ```111    pub fn with_peel(rule: &Rule, peel: usize) -> Result<Automaton> {112        let base = rule.letters() as u64;113        let states = rule.states();114        if peel < rule.width.max(2) {115            return value_error(format!(116                "peel depth {peel} is under the rule width {} or under two.",117                rule.width118            ));119        }120        if (base as f64).powi(peel as i32) > 9.007_199_254_740_992e15 {121            return value_error(format!("peel depth {peel} overruns the exact integers."));122        }123        let source = transfer(rule);124        let mut arcs: Vec<(usize, usize, u64)> = Vec::new();125        if rule.width == 1 {126            for c in 0..rule.letters() {127                if rule.allowed(c) {128                    arcs.push((0, 0, c as u64));129                }130            }131        } else {132            for s in 0..source.len() {133                for c in 0..rule.letters() {134                    let window = (s << rule.dimension) | c;135                    if rule.allowed(window) {136                        arcs.push((window & (states - 1), s, c as u64));137                    }138                }139            }140        }141        let largest = arcs.iter().map(|&(_, _, a)| a).max().unwrap_or(0);142        let power: Vec<Vec<f64>> = (0..base)143            .map(|a| (0..=CUT_CAP).map(|l| (a as f64).powi(l as i32)).collect())144            .collect();145        let mut low: Vec<f64> = Vec::new();146        let mut mid: Vec<Vec<f64>> = vec![Vec::new(); states];147        let mut seats = vec![0.0f64; states];148        let mut word = Vec::new();149        walk_words(rule, base, peel, &mut word, &mut |len, state, value| {150            if len < peel {151                low.push((value as f64).ln());152            } else {153                mid[state].push((value as f64).ln());154                seats[state] += 1.0;155            }156        });157        if low.is_empty() && seats.iter().all(|&c| c == 0.0) {158            return value_error("the rule accepts no element with a nonzero leading digit.");159        }160        let mut ladder = Automaton {161            rule: *rule,162            base,163            states,164            peel,165            largest,166            arcs,167            power,168            low,169            mid,170            seats,171            guide: vec![1.0; states],172            perron_high: 0.0,173            perron_low: 0.0,174            charpoly: Vec::new(),175            adjugate: Vec::new(),176        };177        ladder.tighten();178        ladder.faddeev();179        Ok(ladder)180    }181182    fn tighten(&mut self) {183        let mut guide = vec![1.0f64; self.states];184        for _ in 0..COLLATZ_WIELANDT {185            let next = self.step_real(0, &guide);186            let peak = next.iter().cloned().fold(0.0f64, f64::max);187            if peak > 1e120 {188                break;189            }190            for (seat, &gain) in guide.iter_mut().zip(next.iter()) {191                *seat += gain;192            }193        }194        let image = self.step_real(0, &guide);195        let mut high = 0.0f64;196        let mut low = f64::INFINITY;197        for (&gain, &seat) in image.iter().zip(guide.iter()) {198            high = high.max(gain / seat);199            low = low.min(gain / seat);200        }201        self.guide = guide;202        self.perron_high = high * (1.0 + 1e-14) + 1e-300;203        self.perron_low = low * (1.0 - 1e-14);204    }205206    fn step_real(&self, power: usize, vector: &[f64]) -> Vec<f64> {207        let mut out = vec![0.0f64; self.states];208        for &(t, s, a) in &self.arcs {209            out[t] += self.power[a as usize][power] * vector[s];210        }211        out212    }213214    fn step_complex(&self, power: usize, vector: &[Complex]) -> Vec<Complex> {215        let mut out = vec![Complex::new(0.0, 0.0); self.states];216        for &(t, s, a) in &self.arcs {217            out[t] = out[t] + vector[s] * self.power[a as usize][power];218        }219        out220    }221222    /// Returns the rule.223    pub fn rule(&self) -> Rule {224        self.rule225    }226227    /// Returns the base `q = 2^D`.228    pub fn base(&self) -> u64 {229        self.base230    }231232    /// Returns the state count `q^(k-1)`.233    pub fn states(&self) -> usize {234        self.states235    }236237    /// Returns the peel depth `P`.238    pub fn peel(&self) -> usize {239        self.peel240    }241242    /// Returns the transfer matrix `T = Gamma_0` the ladder runs on, the transpose of [`crate::memory::transfer`], entry `(u', u)` counting the letters carrying `u` to `u'`.243    ///244    /// ```245    /// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();246    /// let ladder = mrlynum::automaton::Automaton::new(&golden).unwrap();247    /// assert_eq!(ladder.matrix(), vec![vec![1.0, 1.0], vec![1.0, 0.0]]);248    /// ```249    pub fn matrix(&self) -> Vec<Vec<f64>> {250        let mut out = vec![vec![0.0f64; self.states]; self.states];251        for &(t, s, _) in &self.arcs {252            out[t][s] += 1.0;253        }254        out255    }256257    /// Returns the Collatz-Wielandt bracket `(low, high)` of the Perron root of the transfer matrix, the ratios the ladder divides with.258    ///259    /// For any nonnegative `T` and any `v > 0`, `min_u (T v)_u / v_u <= rho <= max_u (T v)_u / v_u`, so the pair brackets the root with no primitivity hypothesis; the guide is `v = (I + T)^m 1` for `m` sixty, or fewer when the entries would overrun the exponent range.260    /// The ends carry a relative `1e-14` allowance of the same measured kind as [`ROUNDING`] and are not a rounding-certified interval, and the bracket is loose exactly when `T` is defective: at code `13` and `k = 2`, a Jordan block of root `1`, the upper end is `64/62`. Only the upper end is load bearing, as the `theta < 1` test that opens the Neumann branch; the root itself is [`crate::memory::perron`].261    pub fn perron(&self) -> (f64, f64) {262        (self.perron_low, self.perron_high)263    }264265    /// Returns the abscissa `alpha = log_q rho`, with `rho` the exact Perron root of [`crate::memory::perron`].266    ///267    /// The series converges absolutely on `Re s > alpha`. The bracket of [`Automaton::perron`] is a bound and not the root, so it is not read here.268    ///269    /// ```270    /// let jordan = mrlynum::memory::Rule::new(1, 2, 13).unwrap();271    /// let ladder = mrlynum::automaton::Automaton::new(&jordan).unwrap();272    /// assert_eq!(ladder.abscissa(), 0.0);273    /// ```274    pub fn abscissa(&self) -> f64 {275        let root = crate::memory::perron(&self.rule);276        if root <= 0.0 {277            return f64::NEG_INFINITY;278        }279        root.ln() / (self.base as f64).ln()280    }281282    /// Returns the pole spacing `2 pi / log q`.283    pub fn period(&self) -> f64 {284        2.0 * std::f64::consts::PI / (self.base as f64).ln()285    }286}287288// THE RESOLVENT289290fn invert(matrix: &[Vec<Complex>]) -> Option<Vec<Vec<Complex>>> {291    let n = matrix.len();292    let mut a: Vec<Vec<Complex>> = matrix.to_vec();293    let mut b: Vec<Vec<Complex>> = (0..n)294        .map(|i| {295            (0..n)296                .map(|j| Complex::new(f64::from(u8::from(i == j)), 0.0))297                .collect()298        })299        .collect();300    let scale = matrix301        .iter()302        .flat_map(|row| row.iter())303        .map(|z| z.abs())304        .fold(0.0f64, f64::max)305        .max(1.0);306    for col in 0..n {307        let mut pick = col;308        for row in col + 1..n {309            if a[row][col].abs() > a[pick][col].abs() {310                pick = row;311            }312        }313        if a[pick][col].abs() < PIVOT_FLOOR * scale {314            return None;315        }316        a.swap(col, pick);317        b.swap(col, pick);318        let inv = Complex::new(1.0, 0.0) / a[col][col];319        for j in 0..n {320            a[col][j] = a[col][j] * inv;321            b[col][j] = b[col][j] * inv;322        }323        let prow = a[col].clone();324        let qrow = b[col].clone();325        for row in 0..n {326            if row == col {327                continue;328            }329            let factor = a[row][col];330            if factor.abs() == 0.0 {331                continue;332            }333            for j in 0..n {334                a[row][j] = a[row][j] - factor * prow[j];335                b[row][j] = b[row][j] - factor * qrow[j];336            }337        }338    }339    Some(b)340}341342fn log_binomial(top: f64, pick: usize) -> f64 {343    (1..=pick)344        .map(|i| ((top - pick as f64 + i as f64) / i as f64).ln())345        .sum()346}347348fn poly(logs: &[f64], w: Complex) -> Complex {349    let mut acc = Complex::new(0.0, 0.0);350    for &l in logs {351        acc = acc + ((-w) * l).exp();352    }353    acc354}355356fn poly_scale(logs: &[f64], sigma: f64) -> f64 {357    logs.iter().map(|&l| (-sigma * l).exp()).sum()358}359360struct Rung {361    value: Vec<Complex>,362    bound: Vec<f64>,363    scale: Vec<f64>,364}365366impl Automaton {367    fn neumann(&self, sigma: f64, seed: &[f64]) -> Option<Vec<f64>> {368        let base = self.base as f64;369        let theta = base.powf(-sigma) * self.perron_high;370        if theta >= 1.0 {371            return None;372        }373        let factor = base.powf(-sigma);374        let beta = seed375            .iter()376            .zip(self.guide.iter())377            .map(|(&y, &v)| y / v)378            .fold(0.0f64, f64::max);379        let mut out = vec![0.0f64; self.states];380        let mut term = seed.to_vec();381        let mut power = 1.0f64;382        for _ in 0..4000 {383            for (seat, &gain) in out.iter_mut().zip(term.iter()) {384                *seat += gain;385            }386            term = self387                .step_real(0, &term)388                .iter()389                .map(|&v| v * factor)390                .collect();391            power *= theta;392            let peak = out.iter().cloned().fold(0.0f64, f64::max);393            if beta * power <= 1e-18 * (1.0 - theta) * peak || term.iter().all(|&v| v == 0.0) {394                break;395            }396        }397        let rest = beta * power / (1.0 - theta);398        for (seat, &v) in out.iter_mut().zip(self.guide.iter()) {399            *seat += rest * v;400        }401        Some(out)402    }403404    fn tail_parts(&self, sigma: f64) -> Option<(f64, Vec<f64>)> {405        let mass = self.neumann(sigma, &self.seats)?;406        let log = -(self.peel as f64 - 1.0) * sigma * (self.base as f64).ln();407        Some((log, mass))408    }409410    fn tail(&self, sigma: f64) -> Option<Vec<f64>> {411        let (log, mass) = self.tail_parts(sigma)?;412        Some(mass.iter().map(|&m| log.exp() * m).collect())413    }414415    fn cut(&self, w: Complex, depth: usize) -> Vec<f64> {416        let blown = vec![f64::INFINITY; self.states];417        if self.largest == 0 {418            return vec![0.0; self.states];419        }420        let (span, sigma) = (w.abs(), w.re);421        let base = self.base as f64;422        let largest = self.largest as f64;423        let ratio = ((span + depth as f64 + 1.0) / (depth as f64 + 2.0)).max(1.0) * largest424            / base.powi(self.peel as i32);425        if ratio >= RATIO_CAP {426            return blown;427        }428        let Some((rest, mass)) = self.tail_parts(sigma + depth as f64 + 1.0) else {429            return blown;430        };431        let log = log_binomial(span + depth as f64, depth + 1)432            + (-sigma - depth as f64 - 1.0) * base.ln()433            + (depth as f64 + 1.0) * largest.ln()434            + rest;435        let front = self.step_real(0, &mass);436        front437            .iter()438            .map(|&m| log.exp() * m / (1.0 - ratio))439            .collect()440    }441442    fn divide(&self, w: Complex, acc: &[Complex], err: &[f64], wide: &[f64]) -> Option<Rung> {443        let base = self.base as f64;444        let x = raise(base, -w);445        let front: Vec<Vec<Complex>> = self446            .matrix()447            .iter()448            .enumerate()449            .map(|(i, row)| {450                row.iter()451                    .enumerate()452                    .map(|(j, &t)| Complex::new(f64::from(u8::from(i == j)), 0.0) - x * t)453                    .collect()454            })455            .collect();456        let inverse = invert(&front)?;457        let value: Vec<Complex> = inverse458            .iter()459            .map(|row| {460                row.iter()461                    .zip(acc.iter())462                    .fold(Complex::new(0.0, 0.0), |a, (&c, &y)| a + c * y)463            })464            .collect();465        let carry = |seed: &[f64]| -> Option<Vec<f64>> {466            if let Some(out) = self.neumann(w.re, seed) {467                return Some(out);468            }469            let mut residual = 0.0f64;470            for (i, source) in front.iter().enumerate() {471                let mut row = 0.0f64;472                for j in 0..self.states {473                    let mut entry = Complex::new(f64::from(u8::from(i == j)), 0.0);474                    for (&near, image) in source.iter().zip(inverse.iter()) {475                        entry = entry - near * image[j];476                    }477                    row += entry.abs();478                }479                residual = residual.max(row);480            }481            if residual >= 1.0 {482                return None;483            }484            let peak = seed.iter().cloned().fold(0.0f64, f64::max);485            let lift = peak * residual / (1.0 - residual);486            Some(487                inverse488                    .iter()489                    .map(|row| {490                        row.iter()491                            .zip(seed.iter())492                            .map(|(&c, &y)| c.abs() * (y + lift))493                            .sum()494                    })495                    .collect(),496            )497        };498        Some(Rung {499            value,500            bound: carry(err)?,501            scale: carry(wide)?,502        })503    }504505    fn rung(&self, s: Complex, shift: f64, depth: usize, whole: bool) -> Option<Rung> {506        let base = self.base as f64;507        let levels = (shift - s.re).ceil().max(0.0) as usize;508        if !whole && levels == 0 {509            let rest = self.tail(s.re)?;510            let lift = self.step_real(0, &rest);511            let factor = base.powf(-s.re);512            let out: Vec<f64> = rest513                .iter()514                .zip(lift.iter())515                .map(|(&r, &l)| r + factor * l)516                .collect();517            return Some(Rung {518                value: vec![Complex::new(0.0, 0.0); self.states],519                bound: out.clone(),520                scale: out,521            });522        }523        let rows = levels + depth + 2;524        let mut value = vec![vec![Complex::new(0.0, 0.0); self.states]; rows];525        let mut bound = vec![vec![0.0f64; self.states]; rows];526        let mut scale = vec![vec![0.0f64; self.states]; rows];527        for j in levels..rows {528            let rest = self.tail(s.re + j as f64)?;529            bound[j] = rest.clone();530            scale[j] = rest;531        }532        for j in (0..levels).rev() {533            let w = s + j as f64;534            let mut acc: Vec<Complex> = self535                .mid536                .iter()537                .map(|logs| poly(logs, w))538                .collect::<Vec<Complex>>();539            let mut err = self.cut(w, depth);540            let mut wide: Vec<f64> = self541                .mid542                .iter()543                .map(|logs| poly_scale(logs, w.re))544                .zip(err.iter())545                .map(|(p, &e)| p + e)546                .collect();547            let mut coefficient = Complex::new(1.0, 0.0);548            for l in 1..=depth {549                coefficient = coefficient * ((-w) - (l as f64 - 1.0)) * (1.0 / l as f64);550                let weight = raise(base, (-w) - l as f64) * coefficient;551                let lift = self.step_complex(l, &value[j + l]);552                let heavy = self.step_real(l, &bound[j + l]);553                let broad = self.step_real(l, &scale[j + l]);554                for u in 0..self.states {555                    acc[u] = acc[u] + weight * lift[u];556                    err[u] += weight.abs() * heavy[u];557                    wide[u] += weight.abs() * broad[u];558                }559            }560            if !whole && j == 0 {561                return Some(Rung {562                    value: acc,563                    bound: err,564                    scale: wide,565                });566            }567            let step = self.divide(w, &acc, &err, &wide)?;568            value[j] = step.value;569            bound[j] = step.bound;570            scale[j] = step.scale;571        }572        Some(Rung {573            value: value[0].clone(),574            bound: bound[0].clone(),575            scale: scale[0].clone(),576        })577    }578579    fn tune(580        &self,581        s: Complex,582        tolerance: f64,583        whole: bool,584        head: f64,585    ) -> Result<(Vec<Complex>, Vec<f64>)> {586        let mut shift = SHIFT_START;587        let mut depth = CUT_START;588        let mut best = f64::INFINITY;589        let mut walked = false;590        while shift <= SHIFT_CAP && depth <= CUT_CAP {591            if let Some(step) = self.rung(s, shift, depth, whole) {592                let carried: Vec<f64> = step593                    .bound594                    .iter()595                    .zip(step.scale.iter())596                    .map(|(&b, &c)| b + ROUNDING * c)597                    .collect();598                let total: f64 = carried.iter().sum::<f64>() + ROUNDING * head;599                if total < tolerance {600                    return Ok((step.value, carried));601                }602                best = best.min(total);603            } else {604                walked = true;605            }606            shift += SHIFT_STEP;607            depth += CUT_STEP;608        }609        if walked {610            return value_error(format!(611                "the ladder cannot invert I - q^(-w) T on the walk up from s = {} + {}i: a level sits on a pole of the resolvent or left of the abscissa with no certified inverse.",612                s.re, s.im613            ));614        }615        value_error(format!(616            "tolerance {tolerance:e} is out of reach at s = {} + {}i, best bound {best:e}.",617            s.re, s.im618        ))619    }620}621622// THE DETERMINANT AND THE ADJUGATE623624impl Automaton {625    fn faddeev(&mut self) {626        let n = self.states;627        let matrix = self.matrix();628        let mut blocks: Vec<Vec<Vec<f64>>> = Vec::with_capacity(n);629        let mut coefficients = vec![1.0f64];630        let mut current: Vec<Vec<f64>> = (0..n)631            .map(|i| (0..n).map(|j| f64::from(u8::from(i == j))).collect())632            .collect();633        for k in 1..=n {634            blocks.push(current.clone());635            let mut product = vec![vec![0.0f64; n]; n];636            for i in 0..n {637                for t in 0..n {638                    if matrix[i][t] == 0.0 {639                        continue;640                    }641                    for j in 0..n {642                        product[i][j] += matrix[i][t] * current[t][j];643                    }644                }645            }646            let trace: f64 = (0..n).map(|i| product[i][i]).sum();647            let coefficient = -trace / k as f64;648            coefficients.push(coefficient);649            for (i, row) in product.iter_mut().enumerate() {650                row[i] += coefficient;651            }652            current = product;653        }654        self.charpoly = coefficients;655        self.adjugate = blocks;656    }657658    /// Returns the coefficients `c_0 .. c_n` of `det(I - x T) = sum c_i x^i`, the ladder denominator read as a polynomial in `x = q^(-s)`.659    ///660    /// The roots in `x` are the reciprocals of the nonzero eigenvalues of `T`, so the poles of `zeta_W` sit at `s = log_q lambda_i - m + (arg lambda_i + 2 pi j) i / log q`: one comb per distinct eigenvalue, offset by the eigenvalue's argument.661    ///662    /// ```663    /// let golden = mrlynum::memory::Rule::new(1, 2, 7).unwrap();664    /// let ladder = mrlynum::automaton::Automaton::new(&golden).unwrap();665    /// let coefficients = ladder.denominator();666    /// assert!((coefficients[1] + 1.0).abs() < 1e-12 && (coefficients[2] + 1.0).abs() < 1e-12);667    /// ```668    pub fn denominator(&self) -> Vec<f64> {669        self.charpoly.clone()670    }671672    fn det(&self, x: Complex) -> Complex {673        let mut acc = Complex::new(0.0, 0.0);674        let mut power = Complex::new(1.0, 0.0);675        for &c in &self.charpoly {676            acc = acc + power * c;677            power = power * x;678        }679        acc680    }681682    fn det_slope(&self, x: Complex) -> Complex {683        let mut acc = Complex::new(0.0, 0.0);684        let mut power = Complex::new(1.0, 0.0);685        for (k, &c) in self.charpoly.iter().enumerate().skip(1) {686            acc = acc + power * (c * k as f64);687            power = power * x;688        }689        acc690    }691692    fn adjugate_apply(&self, x: Complex, vector: &[Complex]) -> Vec<Complex> {693        let mut out = vec![Complex::new(0.0, 0.0); self.states];694        let mut power = Complex::new(1.0, 0.0);695        for block in &self.adjugate {696            for (i, row) in block.iter().enumerate() {697                for (j, &entry) in row.iter().enumerate() {698                    if entry != 0.0 {699                        out[i] = out[i] + power * entry * vector[j];700                    }701                }702            }703            power = power * x;704        }705        out706    }707708    fn adjugate_weights(&self, radius: f64) -> Vec<f64> {709        let mut out = vec![0.0f64; self.states];710        let mut power = 1.0f64;711        for block in &self.adjugate {712            for row in block.iter() {713                for (j, &entry) in row.iter().enumerate() {714                    out[j] += power * entry.abs();715                }716            }717            power *= radius;718        }719        out720    }721}722723// READINGS724725impl Automaton {726    /// Returns `zeta_W(s)` and the bound it is known to, or an error when the tolerance is out of reach or a level of the walk sits on a pole of the resolvent.727    ///728    /// The bound is the propagated truncation bound, carried entrywise as a nonnegative vector through `(I - q^(-w) T)^(-1)`, plus [`ROUNDING`] times the scale the ladder carries, which is a measured allowance and not a proof.729    /// Right of the abscissa the inverse is majorised by its Neumann series, `sum_i (q^(-Re w) T)^i`, which is entrywise nonnegative and needs no norm and no primitivity; left of it the bound runs through the computed inverse certified by its own residual, `abs(C) (y + norm(y) r/(1-r) 1)` with `r = norm(I - (I - q^(-w) T) C)`, and the module raises rather than return when `r >= 1`.730    ///731    /// ```732    /// let full = mrlynum::memory::Rule::full(1, 2);733    /// let ladder = mrlynum::automaton::Automaton::new(&full).unwrap();734    /// let s = mrlynum::zeta::Complex::new(2.0, 0.0);735    /// let (value, bound) = ladder.zeta(s, 1e-10).unwrap();736    /// assert!((value.re - std::f64::consts::PI * std::f64::consts::PI / 6.0).abs() < bound);737    /// ```738    pub fn zeta(&self, s: Complex, tolerance: f64) -> Result<(Complex, f64)> {739        let head = poly_scale(&self.low, s.re);740        let (value, carried) = self.tune(s, tolerance, true, head)?;741        let sum = value.iter().fold(Complex::new(0.0, 0.0), |acc, &z| acc + z);742        let bound = carried.iter().sum::<f64>() + ROUNDING * head;743        Ok((poly(&self.low, s) + sum, bound))744    }745746    /// Returns the matrix Lyndon cofactor `Z_W(s) = det(I - q^(-s) T) zeta_W(s)` and the bound it is known to.747    ///748    /// The scalar `1 - k q^(-s)` becomes the determinant, and the cofactor is carried as `det(I - q^(-s) T) D_(P-1)(s) + 1^T adj(I - q^(-s) T) N(s)` with `N` the ladder numerator, so it is read on the whole `m = 0` pole comb where `zeta_W` itself is singular.749    pub fn cofactor(&self, s: Complex, tolerance: f64) -> Result<(Complex, f64)> {750        let x = raise(self.base as f64, -s);751        let det = self.det(x);752        let head = det.abs() * poly_scale(&self.low, s.re);753        let (value, carried) = self.tune(s, tolerance, false, head)?;754        let lifted = self.adjugate_apply(x, &value);755        let sum = lifted756            .iter()757            .fold(Complex::new(0.0, 0.0), |acc, &z| acc + z);758        let weights = self.adjugate_weights(x.abs());759        let bound = weights760            .iter()761            .zip(carried.iter())762            .map(|(&m, &b)| m * b)763            .sum::<f64>()764            + ROUNDING * head;765        Ok((det * poly(&self.low, s) + sum, bound))766    }767768    /// Returns the residue of `zeta_W` at a simple pole `w0` of the resolvent and the bound it is known to, or an error when `w0` is not a simple root of `det(I - q^(-w) T)`.769    ///770    /// At such a point the resolvent is `adj(I - x T) / det(I - x T)` with `x = q^(-w)`, so the residue is `1^T adj(I - x0 T) N(w0)` over `-x0 log q det'(x0)`: the adjugate is the spectral projector in polynomial form, and no eigenvector is solved for.771    /// The head term `det(I - q^(-s) T) D_(P-1)(s)` is dropped, since the determinant vanishes at the pole.772    pub fn residue(&self, w0: Complex, tolerance: f64) -> Result<(Complex, f64)> {773        let base = self.base as f64;774        let x = raise(base, -w0);775        let slope = self.det_slope(x);776        let den = -x * base.ln() * slope;777        if self.det(x).abs() > 1e-8 * self.charpoly.iter().map(|c| c.abs()).sum::<f64>().max(1.0) {778            return value_error(format!(779                "w0 = {} + {}i is not a root of det(I - q^(-w) T).",780                w0.re, w0.im781            ));782        }783        if den.abs() < 1e-12 {784            return value_error(format!(785                "the root at w0 = {} + {}i is not simple, so the pole order exceeds one.",786                w0.re, w0.im787            ));788        }789        let (value, carried) = self.tune(w0, tolerance, false, 0.0)?;790        let lifted = self.adjugate_apply(x, &value);791        let sum = lifted792            .iter()793            .fold(Complex::new(0.0, 0.0), |acc, &z| acc + z);794        let weights = self.adjugate_weights(x.abs());795        let bound = weights796            .iter()797            .zip(carried.iter())798            .map(|(&m, &b)| m * b)799            .sum::<f64>()800            / den.abs();801        Ok((sum / den, bound))802    }803}804805#[cfg(test)]806mod tests {807    use super::*;808    use crate::ladder::Design;809810    const PINNED: f64 = 1e-9;811812    fn near(value: Complex, re: f64, im: f64, bound: f64) {813        assert!(814            (value.re - re).abs() < bound && (value.im - im).abs() < bound,815            "{value:?} is not {re} + {im}i inside {bound:e}"816        );817    }818819    fn golden() -> Automaton {820        Automaton::new(&Rule::new(1, 2, 7).unwrap()).unwrap()821    }822823    #[test]824    fn the_full_rule_is_the_base_two_design() {825        let design = Design::new(2, &[0, 1]).unwrap();826        for width in [2usize, 3] {827            let ladder = Automaton::new(&Rule::full(1, width)).unwrap();828            assert_eq!(ladder.states(), 1 << (width - 1));829            assert_eq!(ladder.abscissa(), 1.0);830            let mut wanted = vec![0.0; 1 << (width - 1)];831            wanted[0] = 1.0;832            if width > 1 {833                wanted[1] = -2.0;834            }835            wanted.push(0.0);836            assert_eq!(ladder.denominator(), wanted);837            for s in [Complex::new(2.0, 0.0), Complex::new(0.3, 40.0)] {838                let (mine, one) = ladder.zeta(s, 1e-8).unwrap();839                let (theirs, two) = crate::ladder::zeta(&design, s, 1e-8).unwrap();840                assert!(841                    (mine - theirs).abs() < one + two,842                    "width {width} at {s:?}: {mine:?} {theirs:?} {one:e} {two:e}"843                );844            }845        }846    }847848    #[test]849    fn the_product_rule_is_the_mersenne_series() {850        let ladder = Automaton::new(&Rule::new(1, 2, 8).unwrap()).unwrap();851        assert!(ladder.abscissa().abs() < 1e-12);852        for s in [Complex::new(2.0, 0.0), Complex::new(1.5, 3.0)] {853            let (value, bound) = ladder.zeta(s, 1e-9).unwrap();854            let mut want = Complex::new(0.0, 0.0);855            for m in 1..=60u32 {856                want = want + ((-s) * ((2f64.powi(m as i32) - 1.0).ln())).exp();857            }858            assert!(859                (value - want).abs() < bound + 1e-15,860                "{s:?} {value:?} {want:?}"861            );862        }863    }864865    #[test]866    fn the_golden_rule_is_the_fibbinary_series() {867        let ladder = golden();868        assert_eq!(ladder.peel(), 12);869        let (value, bound) = ladder.zeta(Complex::new(3.0, 0.0), PINNED).unwrap();870        let (mut want, mut drift) = (0.0f64, 0.0f64);871        for n in 1u64..1 << 22 {872            if n & (n << 1) == 0 {873                let term = (n as f64).powf(-3.0) - drift;874                let total = want + term;875                drift = (total - want) - term;876                want = total;877            }878        }879        let (mut a, mut b) = (1.0f64, 1.0f64);880        let mut rest = 0.0f64;881        for _ in 1..=22u32 {882            let next = a + b;883            a = b;884            b = next;885        }886        for step in 0..400u32 {887            rest += a * 2f64.powf(-3.0 * (22 + step) as f64);888            let next = a + b;889            a = b;890            b = next;891        }892        assert!(893            value.re > want - bound && value.re < want + rest + bound,894            "{value:?} {want} {rest:e} {bound:e}"895        );896    }897898    #[test]899    fn the_golden_transfer_carries_the_two_combs() {900        let ladder = golden();901        let phi = (1.0 + 5f64.sqrt()) / 2.0;902        let (low, high) = ladder.perron();903        assert!(low <= phi && phi <= high, "{low} {phi} {high}");904        assert!(high - low < 1e-9, "{low} {high}");905        let coefficients = ladder.denominator();906        assert_eq!(coefficients.len(), 3);907        assert!((coefficients[1] + 1.0).abs() < 1e-12);908        assert!((coefficients[2] + 1.0).abs() < 1e-12);909        assert_eq!(ladder.abscissa(), phi.log2());910    }911912    #[test]913    fn the_residues_on_the_first_comb_are_the_contour_averages() {914        let ladder = golden();915        let phi = (1.0 + 5f64.sqrt()) / 2.0;916        for j in 0..3i64 {917            let w0 = Complex::new(phi.log2(), ladder.period() * j as f64);918            let (want, bound) = ladder.residue(w0, PINNED).unwrap();919            let (radius, nodes) = (0.05, 8);920            let mut got = Complex::new(0.0, 0.0);921            let mut carried = 0.0;922            for k in 0..nodes {923                let turn = Complex::turn(std::f64::consts::TAU * k as f64 / nodes as f64);924                let (value, other) = ladder.zeta(w0 + turn * radius, 1e-9).unwrap();925                got = got + turn * radius * value * (1.0 / nodes as f64);926                carried += radius * other / nodes as f64;927            }928            assert!(929                (got - want).abs() < bound + carried + 1e-6,930                "j {j} {got:?} {want:?}"931            );932        }933    }934935    #[test]936    fn the_cofactor_is_the_series_times_the_determinant() {937        let ladder = golden();938        for s in [Complex::new(2.0, 0.0), Complex::new(1.2, 9.0)] {939            let x = raise(2.0, -s);940            let mut det = Complex::new(1.0, 0.0);941            let mut power = Complex::new(1.0, 0.0);942            for &c in ladder.denominator().iter().skip(1) {943                power = power * x;944                det = det + power * c;945            }946            let (value, bound) = ladder.cofactor(s, 1e-8).unwrap();947            let (series, other) = ladder.zeta(s, 1e-8).unwrap();948            assert!(949                (value - det * series).abs() < bound + det.abs() * other,950                "{s:?} {value:?} {:?}",951                det * series952            );953        }954    }955956    #[test]957    fn the_golden_residues_on_the_first_comb_are_pinned() {958        let ladder = golden();959        let phi = (1.0 + 5f64.sqrt()) / 2.0;960        let wanted = [961            (0.946_743_395_641_970_1, 0.0),962            (0.210_170_579_042_707_6, -0.581_938_842_807_736_5),963            (0.062_192_494_764_691_71, -0.051_732_573_832_334_72),964        ];965        for (j, (re, im)) in wanted.iter().enumerate() {966            let w0 = Complex::new(phi.log2(), ladder.period() * j as f64);967            let (value, bound) = ladder.residue(w0, PINNED).unwrap();968            assert!(bound < 1.3e-13, "j {j} bound {bound:e}");969            near(value, *re, *im, bound);970        }971    }972973    #[test]974    fn the_second_comb_is_genuine() {975        let ladder = golden();976        let phi = (1.0 + 5f64.sqrt()) / 2.0;977        let half = std::f64::consts::PI / 2f64.ln();978        let wanted = [979            (-0.259_501_222_742_937_13, -0.592_535_006_433_179_4),980            (0.896_350_590_641_920_6, 1.403_072_744_223_695_2),981            (0.491_379_388_883_392_94, -3.790_264_223_035_055_4),982        ];983        for (j, (re, im)) in wanted.iter().enumerate() {984            let w0 = Complex::new(-phi.log2(), half * (2 * j + 1) as f64);985            let (value, bound) = ladder.residue(w0, 1e-7).unwrap();986            assert!(value.abs() > 0.6, "j {j} residue {value:?}");987            near(value, *re, *im, bound);988            let (radius, nodes) = (0.05, 12);989            let mut got = Complex::new(0.0, 0.0);990            let mut carried = 0.0;991            for k in 0..nodes {992                let turn = Complex::turn(std::f64::consts::TAU * k as f64 / nodes as f64);993                let (point, other) = ladder.zeta(w0 + turn * radius, 1e-5).unwrap();994                got = got + turn * radius * point * (1.0 / nodes as f64);995                carried += radius * other / nodes as f64;996            }997            assert!(998                (got - value).abs() < bound + carried + 1e-6,999                "j {j} {got:?} {value:?}"1000            );1001        }1002    }10031004    #[test]1005    fn the_residue_at_the_abscissa_is_the_limit_of_the_digit_sums() {1006        let ladder = golden();1007        let alpha = ladder.abscissa();1008        let (value, bound) = ladder.residue(Complex::new(alpha, 0.0), PINNED).unwrap();1009        let wanted = [1010            (16usize, 0.946_747_043_404_283_4, 4e-6),1011            (20, 0.946_743_630_023_051_8, 3e-7),1012            (24, 0.946_743_410_426_742, 2e-8),1013        ];1014        for (length, want, reach) in wanted {1015            let (mut sum, mut drift) = (0.0f64, 0.0f64);1016            for n in (1u64 << (length - 1))..(1u64 << length) {1017                if n & (n << 1) == 0 {1018                    let term = (n as f64).powf(-alpha) - drift;1019                    let total = sum + term;1020                    drift = (total - sum) - term;1021                    sum = total;1022                }1023            }1024            let got = sum / 2f64.ln();1025            assert!((got - want).abs() < 1e-12, "L {length} {got}");1026            assert!(1027                (got - value.re).abs() < reach + bound,1028                "L {length} {got} {value:?}"1029            );1030        }1031    }10321033    #[test]1034    fn the_golden_readings_are_pinned() {1035        let ladder = golden();1036        let pinned = [1037            (Complex::new(3.0, 0.0), 1.154_012_963_277_640_1, 0.0),1038            (Complex::new(2.0, 0.0), 1.415_825_532_884_778_4, 0.0),1039            (1040                Complex::new(1.2, 9.0),1041                1.906_409_024_243_908_5,1042                -0.453_243_424_778_265_46,1043            ),1044            (Complex::new(0.8, 0.0), 9.536_379_694_275_015, 0.0),1045        ];1046        for (s, re, im) in pinned {1047            let (value, bound) = ladder.zeta(s, PINNED).unwrap();1048            near(value, re, im, bound);1049        }1050        let cofactors = [1051            (Complex::new(3.0, 0.0), 0.991_729_890_316_722),1052            (Complex::new(2.0, 0.0), 0.973_380_053_858_285_1),1053            (Complex::new(0.8, 0.0), 0.913_335_748_872_126_1),1054        ];1055        for (s, re) in cofactors {1056            let (value, bound) = ladder.cofactor(s, PINNED).unwrap();1057            near(value, re, 0.0, bound);1058        }1059    }10601061    #[test]1062    fn the_fibbinary_set_carries_no_euler_product() {1063        let fib = |n: u64| n & (n << 1) == 0;1064        let mut found = None;1065        'outer: for product in 2u64..1 << 16 {1066            for a in 2u64..=(product as f64).sqrt() as u64 + 1 {1067                if product % a != 0 {1068                    continue;1069                }1070                let b = product / a;1071                if a >= b || !fib(a) || !fib(b) {1072                    continue;1073                }1074                let (mut x, mut y) = (a, b);1075                while y != 0 {1076                    let t = x % y;1077                    x = y;1078                    y = t;1079                }1080                if x != 1 || fib(product) {1081                    continue;1082                }1083                found = Some((a, b, product));1084                break 'outer;1085            }1086        }1087        assert_eq!(found, Some((5, 9, 45)));1088    }10891090    #[test]1091    fn the_abscissa_is_the_exact_perron_root_and_not_the_bracket() {1092        let jordan = Automaton::new(&Rule::new(1, 2, 13).unwrap()).unwrap();1093        assert_eq!(jordan.matrix(), vec![vec![1.0, 1.0], vec![0.0, 1.0]]);1094        assert_eq!(jordan.abscissa(), 0.0);1095        let (_, high) = jordan.perron();1096        assert!(high > 1.03 && high < 1.04, "{high}");1097    }10981099    #[test]1100    fn the_ladder_refuses_a_level_it_cannot_invert() {1101        let ladder = golden();1102        let phi = (1.0 + 5f64.sqrt()) / 2.0;1103        let w0 = Complex::new(phi.log2(), 0.0);1104        assert!(ladder.zeta(w0, 1e-6).is_err());1105        assert!(ladder.residue(Complex::new(2.0, 0.0), 1e-6).is_err());1106    }1107}