carry.rs

1use mrlycore::errors::{value_error, Result};23const CEILING: usize = 64;4const SCAN: usize = 8192;5const HALVINGS: usize = 200;6const SWEEPS: usize = 600;78fn overflow<T>() -> Result<T> {9    value_error("the carry arithmetic passes a hundred and twenty-eight bits.")10}1112fn product(a: i128, b: i128) -> Result<i128> {13    match a.checked_mul(b) {14        Some(value) => Ok(value),15        None => overflow(),16    }17}1819fn total(a: i128, b: i128) -> Result<i128> {20    match a.checked_add(b) {21        Some(value) => Ok(value),22        None => overflow(),23    }24}2526fn middle(base: usize) -> Result<usize> {27    if base < 3 || base.is_multiple_of(2) {28        return value_error("the base must be odd and at least three.");29    }30    Ok((base - 1) / 2)31}3233fn sized(dimension: usize) -> Result<()> {34    if !(2..=CEILING).contains(&dimension) {35        return value_error(format!("the dimension must be between 2 and {CEILING}."));36    }37    Ok(())38}3940// POLYNOMIAL4142fn convolve(left: &[i128], right: &[i128]) -> Result<Vec<i128>> {43    let mut out = vec![0i128; left.len() + right.len() - 1];44    for (i, &x) in left.iter().enumerate() {45        if x == 0 {46            continue;47        }48        for (j, &y) in right.iter().enumerate() {49            if y == 0 {50                continue;51            }52            out[i + j] = total(out[i + j], product(x, y)?)?;53        }54    }55    Ok(out)56}5758/// The digit polynomial of the base-`q` middle-digit design in dimension `D`, lowest power first.59///60/// The design keeps the cells whose digit vector has at most one coordinate equal to the middle61/// digit, so its weight generating function is `A(t)^(D-1) (A(t) + D t^m)`, with `A` the sum of62/// every digit power but the middle one. At base three that factors as63/// `(1 + t^2)^(D-1) (1 + D t + t^2)`, and the sum of the coefficients is the fill.64///65/// ```66/// assert_eq!(mrlymath::dim::carry::digit_polynomial(3, 3).unwrap(), vec![1, 3, 3, 6, 3, 3, 1]);67/// ```68pub fn digit_polynomial(base: usize, dimension: usize) -> Result<Vec<i128>> {69    let centre = middle(base)?;70    sized(dimension)?;71    let alpha: Vec<i128> = (0..base).map(|digit| i128::from(digit != centre)).collect();72    let mut poly = vec![1i128];73    for _ in 1..dimension {74        poly = convolve(&poly, &alpha)?;75    }76    let mut second = alpha;77    second[centre] += dimension as i128;78    convolve(&poly, &second)79}8081/// The count of level-one cells the design keeps, `f_D = (q - 1)^(D-1) (q - 1 + D)`.82///83/// ```84/// assert_eq!(mrlymath::dim::carry::fill(3, 3).unwrap(), 20);85/// assert_eq!(mrlymath::dim::carry::fill(5, 3).unwrap(), 112);86/// ```87pub fn fill(base: usize, dimension: usize) -> Result<i128> {88    middle(base)?;89    sized(dimension)?;90    let mut out = 1i128;91    for _ in 1..dimension {92        out = product(out, base as i128 - 1)?;93    }94    product(out, (base - 1 + dimension) as i128)95}9697// MATRIX9899fn coefficient(poly: &[i128], index: i128) -> i128 {100    if index < 0 || index as usize >= poly.len() {101        0102    } else {103        poly[index as usize]104    }105}106107/// The carry matrix over the reachable carries `|c| <= (D-1)/2`, rows indexed by the carry out.108///109/// A level of the design adds one base-`q` digit per coordinate, so the height of the central110/// diagonal hyperplane moves by a digit sum `s` and a carry `c -> (c + mD - s)/q` with `m` the111/// middle digit; the map contracts onto this window from every start.112pub fn carry_matrix(base: usize, dimension: usize) -> Result<Vec<Vec<i128>>> {113    let poly = digit_polynomial(base, dimension)?;114    let shift = (middle(base)? * dimension) as i128;115    let q = base as i128;116    let half = ((dimension - 1) / 2) as i128;117    Ok((-half..=half)118        .map(|out| {119            (-half..=half)120                .map(|inside| coefficient(&poly, inside + shift - q * out))121                .collect()122        })123        .collect())124}125126/// The reflection-even block of the carry matrix, of size `ceil(D/2)`.127///128/// The digit polynomial is palindromic, so the reflection `c -> -c` commutes with the carry map129/// and splits it; the even half carries the central count and its characteristic polynomial is130/// the recurrence the counts obey.131///132/// ```133/// assert_eq!(mrlymath::dim::carry::even_block(3, 3).unwrap(), vec![vec![6, 6], vec![1, 3]]);134/// ```135pub fn even_block(base: usize, dimension: usize) -> Result<Vec<Vec<i128>>> {136    let poly = digit_polynomial(base, dimension)?;137    let shift = (middle(base)? * dimension) as i128;138    let q = base as i128;139    let width = (dimension - 1) / 2 + 1;140    Ok((0..width)141        .map(|out| {142            (0..width)143                .map(|inside| {144                    let step = shift - q * out as i128;145                    let folded = if inside == 0 {146                        0147                    } else {148                        coefficient(&poly, step - inside as i128)149                    };150                    coefficient(&poly, step + inside as i128) + folded151                })152                .collect()153        })154        .collect())155}156157/// The trace of a square integer matrix.158///159/// ```160/// let block = mrlymath::dim::carry::even_block(3, 3).unwrap();161/// assert_eq!(mrlymath::dim::carry::trace(&block), 9);162/// ```163pub fn trace(rows: &[Vec<i128>]) -> i128 {164    (0..rows.len()).map(|index| rows[index][index]).sum()165}166167fn multiply(left: &[Vec<i128>], right: &[Vec<i128>]) -> Result<Vec<Vec<i128>>> {168    let n = left.len();169    let mut out = vec![vec![0i128; n]; n];170    for row in 0..n {171        for step in 0..n {172            let weight = left[row][step];173            if weight == 0 {174                continue;175            }176            for col in 0..n {177                out[row][col] = total(out[row][col], product(weight, right[step][col])?)?;178            }179        }180    }181    Ok(out)182}183184/// The monic characteristic polynomial of a square integer matrix, highest power first.185///186/// Faddeev-LeVerrier in exact integers: every division lands whole, so no fraction is ever needed187/// and the answer is the recurrence's coefficients up to sign.188///189/// ```190/// let block = mrlymath::dim::carry::even_block(3, 3).unwrap();191/// assert_eq!(mrlymath::dim::carry::characteristic(&block).unwrap(), vec![1, -9, 12]);192/// ```193pub fn characteristic(rows: &[Vec<i128>]) -> Result<Vec<i128>> {194    let n = rows.len();195    let mut held: Vec<Vec<i128>> = (0..n)196        .map(|row| (0..n).map(|col| i128::from(row == col)).collect())197        .collect();198    let mut out = vec![1i128];199    for step in 1..=n {200        let walked = multiply(rows, &held)?;201        let sum = trace(&walked);202        if sum % step as i128 != 0 {203            return value_error("the characteristic polynomial left a remainder.");204        }205        let mark = -sum / step as i128;206        out.push(mark);207        held = walked;208        for (index, row) in held.iter_mut().enumerate().take(n) {209            row[index] = total(row[index], mark)?;210        }211    }212    Ok(out)213}214215/// The determinant of a square integer matrix, read off its characteristic polynomial.216///217/// ```218/// let block = mrlymath::dim::carry::even_block(3, 3).unwrap();219/// assert_eq!(mrlymath::dim::carry::determinant(&block).unwrap(), 12);220/// ```221pub fn determinant(rows: &[Vec<i128>]) -> Result<i128> {222    let poly = characteristic(rows)?;223    let last = poly[rows.len()];224    Ok(if rows.len().is_multiple_of(2) {225        last226    } else {227        -last228    })229}230231// LADDER232233/// The counts `a_D(L)` of level-`L` cells meeting the central diagonal hyperplane, from `L = 0`.234///235/// The count is the top-left entry of the `L`-th power of the even block. The walk stops early236/// when the next count would pass a hundred and twenty-eight bits, so the answer runs as far as237/// exact integers reach and no further.238///239/// ```240/// let terms = mrlymath::dim::carry::ladder(3, 3, 6).unwrap();241/// assert_eq!(terms, vec![1, 6, 42, 306, 2250, 16578, 122202]);242/// ```243pub fn ladder(base: usize, dimension: usize, levels: usize) -> Result<Vec<i128>> {244    let block = even_block(base, dimension)?;245    let n = block.len();246    let mut carried = vec![0i128; n];247    carried[0] = 1;248    let mut out = vec![1i128];249    for _ in 0..levels {250        let mut next = vec![0i128; n];251        for row in 0..n {252            let mut sum = 0i128;253            for col in 0..n {254                match block[row][col]255                    .checked_mul(carried[col])256                    .and_then(|part| sum.checked_add(part))257                {258                    Some(value) => sum = value,259                    None => return Ok(out),260                }261            }262            next[row] = sum;263        }264        out.push(next[0]);265        carried = next;266    }267    Ok(out)268}269270// ROOT271272/// The Perron root of a nonnegative square integer matrix.273///274/// The characteristic polynomial is rescaled by the largest row sum, which bounds the root above,275/// then the largest sign change on the unit interval is hunted on a grid and closed by bisection,276/// the same walk the recurrence growth reader takes.277///278/// ```279/// let block = mrlymath::dim::carry::even_block(3, 3).unwrap();280/// let root = mrlymath::dim::carry::perron(&block).unwrap();281/// assert!((root - (9.0 + 33f64.sqrt()) / 2.0).abs() < 1e-9);282/// ```283pub fn perron(rows: &[Vec<i128>]) -> Result<f64> {284    let poly = characteristic(rows)?;285    let bound = rows286        .iter()287        .map(|row| row.iter().sum::<i128>())288        .max()289        .unwrap_or(0) as f64;290    if bound <= 0.0 {291        return value_error("the block has no positive row sum.");292    }293    let mut scaled = Vec::with_capacity(poly.len());294    let mut power = 1.0f64;295    for &term in &poly {296        scaled.push(term as f64 / power);297        power *= bound;298    }299    let value = |y: f64| scaled.iter().fold(0.0f64, |acc, &term| acc * y + term);300    for step in (0..SCAN).rev() {301        let mut lo = step as f64 / SCAN as f64;302        let mut hi = (step + 1) as f64 / SCAN as f64;303        if value(lo) > 0.0 || value(hi) < 0.0 {304            continue;305        }306        for _ in 0..HALVINGS {307            let mid = (lo + hi) / 2.0;308            if value(mid) <= 0.0 {309                lo = mid;310            } else {311                hi = mid;312            }313        }314        return Ok(bound * (lo + hi) / 2.0);315    }316    value_error("the block shows no root inside its row-sum bound.")317}318319// SIGN320321/// The sign of `log_q rho_D - (log_q f_D - 1)`, the slice sign law's reading, in exact integers.322///323/// The characteristic polynomial is evaluated at `f_D/q` with the denominators cleared, so the324/// answer is a comparison of whole numbers and never a rounding. The law reads `(-1)^(D+1)`: the325/// slice exponent stands above the solid's dimension less one at odd `D` and below it at even `D`.326///327/// ```328/// assert_eq!(mrlymath::dim::carry::sign(3, 3).unwrap(), 1);329/// assert_eq!(mrlymath::dim::carry::sign(3, 4).unwrap(), -1);330/// ```331pub fn sign(base: usize, dimension: usize) -> Result<i32> {332    let block = even_block(base, dimension)?;333    let poly = characteristic(&block)?;334    let full = fill(base, dimension)?;335    let q = base as i128;336    let mut acc = 0i128;337    let mut weight = 1i128;338    for &term in &poly {339        acc = total(product(acc, full)?, product(term, weight)?)?;340        weight = product(weight, q)?;341    }342    Ok(match acc {343        acc if acc > 0 => -1,344        acc if acc < 0 => 1,345        _ => 0,346    })347}348349/// The widest dimension the exact carry arithmetic reaches at the base.350///351/// The sign is the heaviest reading: it walks the characteristic polynomial at `f_D/q` with the352/// denominators cleared, so its running value climbs like `f_D^(D/2)` and a hundred and twenty-eight353/// bits run out at dimension fifteen in base three and dimension eleven in base five.354///355/// ```356/// assert_eq!(mrlymath::dim::carry::cap(3).unwrap(), 15);357/// assert_eq!(mrlymath::dim::carry::cap(5).unwrap(), 11);358/// ```359pub fn cap(base: usize) -> Result<usize> {360    middle(base)?;361    let mut top = 0;362    for dimension in 2..=CEILING {363        if sign(base, dimension).is_err() {364            break;365        }366        top = dimension;367    }368    if top < 2 {369        return value_error(format!(370            "base {base} reaches no dimension in exact integers."371        ));372    }373    Ok(top)374}375376// SPECTRUM377378fn ahead(walk: &[Vec<f64>], vector: &[f64]) -> Vec<f64> {379    walk.iter()380        .map(|row| row.iter().zip(vector).map(|(&a, &x)| a * x).sum())381        .collect()382}383384fn behind(walk: &[Vec<f64>], vector: &[f64]) -> Vec<f64> {385    (0..vector.len())386        .map(|col| {387            (0..vector.len())388                .map(|row| walk[row][col] * vector[row])389                .sum()390        })391        .collect()392}393394fn top(vector: &[f64]) -> f64 {395    vector.iter().fold(0.0f64, |peak, &x| peak.max(x.abs()))396}397398fn settle(walk: &[Vec<f64>], left: bool) -> Option<Vec<f64>> {399    let mut vector = vec![1.0f64; walk.len()];400    for _ in 0..SWEEPS {401        let next = if left {402            behind(walk, &vector)403        } else {404            ahead(walk, &vector)405        };406        let peak = top(&next);407        if peak == 0.0 {408            return None;409        }410        vector = next.iter().map(|&x| x / peak).collect();411    }412    Some(vector)413}414415/// The Perron root over the modulus of the second eigenvalue, or none where the block is one wide.416///417/// Power iteration finds the leading pair, Hotelling deflation takes it out and a second walk,418/// reprojected every sweep so rounding cannot bring the leader back, reads the runner up. At base419/// three the ratio falls to `(D + 2)/(D - 2)`, so no fixed spectral gap survives the dimensions.420pub fn spectral_ratio(base: usize, dimension: usize) -> Result<Option<f64>> {421    let block = even_block(base, dimension)?;422    let n = block.len();423    if n < 2 {424        return Ok(None);425    }426    let bound = block427        .iter()428        .map(|row| row.iter().sum::<i128>())429        .max()430        .unwrap_or(1) as f64;431    let walk: Vec<Vec<f64>> = block432        .iter()433        .map(|row| row.iter().map(|&x| x as f64 / bound).collect())434        .collect();435    let (right, left) = match (settle(&walk, false), settle(&walk, true)) {436        (Some(right), Some(left)) => (right, left),437        _ => return Ok(None),438    };439    let lead = top(&ahead(&walk, &right));440    let pair: f64 = left.iter().zip(&right).map(|(&u, &v)| u * v).sum();441    if pair == 0.0 || lead == 0.0 {442        return Ok(None);443    }444    let mut trail: Vec<f64> = (0..n).map(|i| 1.0 + 0.37 * ((i * 7) % 5) as f64).collect();445    let mut second = 0.0f64;446    for _ in 0..SWEEPS {447        let share: f64 = left.iter().zip(&trail).map(|(&u, &x)| u * x).sum::<f64>() / pair;448        let cleared: Vec<f64> = trail449            .iter()450            .zip(&right)451            .map(|(&x, &v)| x - share * v)452            .collect();453        let next = ahead(&walk, &cleared);454        second = top(&next);455        if second == 0.0 {456            return Ok(None);457        }458        trail = next.iter().map(|&x| x / second).collect();459    }460    Ok(Some(lead / second))461}462463#[cfg(test)]464mod tests {465    use super::*;466467    fn brute(base: usize, dimension: usize) -> Vec<i128> {468        let centre = (base - 1) / 2;469        let mut out = vec![0i128; (base - 1) * dimension + 1];470        for mut code in 0..base.pow(dimension as u32) {471            let mut sum = 0;472            let mut middles = 0;473            for _ in 0..dimension {474                let digit = code % base;475                code /= base;476                sum += digit;477                middles += usize::from(digit == centre);478            }479            if middles <= 1 {480                out[sum] += 1;481            }482        }483        out484    }485486    fn central(base: usize, dimension: usize, levels: usize) -> Vec<i128> {487        let poly = digit_polynomial(base, dimension).unwrap();488        let mut out = vec![1i128];489        let mut span = vec![1i128];490        let mut step = 1usize;491        for _ in 0..levels {492            let mut next = vec![0i128; span.len() + (base - 1) * dimension * step];493            for (at, &count) in span.iter().enumerate() {494                if count == 0 {495                    continue;496                }497                for (weight, &many) in poly.iter().enumerate() {498                    next[at + weight * step] += count * many;499                }500            }501            span = next;502            step *= base;503            out.push(span[dimension * (step - 1) / 2]);504        }505        out506    }507508    #[test]509    fn the_digit_polynomial_is_the_enumeration_of_the_kept_cells() {510        for base in [3usize, 5] {511            for dimension in 2..=6 {512                let got = digit_polynomial(base, dimension).unwrap();513                assert_eq!(514                    got,515                    brute(base, dimension),516                    "base {base} dimension {dimension}"517                );518                assert_eq!(519                    got.iter().sum::<i128>(),520                    fill(base, dimension).unwrap(),521                    "base {base} dimension {dimension}"522                );523            }524        }525        assert!(digit_polynomial(4, 3).is_err());526        assert!(digit_polynomial(3, 1).is_err());527    }528529    #[test]530    fn the_three_dimensional_block_is_the_published_anchor() {531        let block = even_block(3, 3).unwrap();532        assert_eq!(block, vec![vec![6, 6], vec![1, 3]]);533        assert_eq!(trace(&block), 9);534        assert_eq!(determinant(&block).unwrap(), 12);535        assert_eq!(characteristic(&block).unwrap(), vec![1, -9, 12]);536        assert_eq!(537            ladder(3, 3, 6).unwrap(),538            vec![1, 6, 42, 306, 2250, 16578, 122202]539        );540        let root = perron(&block).unwrap();541        assert!((root - (9.0 + 33f64.sqrt()) / 2.0).abs() < 1e-9);542        assert!((root.log(3.0) - 1.818_410).abs() < 1e-6);543        assert!((20f64.log(3.0) - 1.0 - 1.726_833).abs() < 1e-6);544    }545546    #[test]547    fn the_traces_read_the_two_closed_forms() {548        let want = [2i128, 9, 11, 60, 47, 336];549        for (index, &term) in want.iter().enumerate() {550            let dimension = index + 2;551            assert_eq!(trace(&even_block(3, dimension).unwrap()), term);552        }553        for dimension in 2..=cap(3).unwrap() {554            let got = trace(&even_block(3, dimension).unwrap());555            let want = if dimension % 2 == 0 {556                3 * 2i128.pow(dimension as u32 - 2) - 1557            } else {558                3 * dimension as i128 * 2i128.pow(dimension as u32 - 3)559            };560            assert_eq!(got, want, "dimension {dimension}");561        }562    }563564    #[test]565    fn the_three_generators_of_the_ladder_agree() {566        for base in [3usize, 5] {567            for dimension in 2..=6 {568                let block = ladder(base, dimension, 4).unwrap();569                assert_eq!(570                    block,571                    central(base, dimension, 4),572                    "base {base} dimension {dimension}"573                );574                let full = carry_matrix(base, dimension).unwrap();575                let half = full.len() / 2;576                let mut walked: Vec<Vec<i128>> = (0..full.len())577                    .map(|row| (0..full.len()).map(|col| i128::from(row == col)).collect())578                    .collect();579                let mut got = vec![1i128];580                for _ in 0..4 {581                    walked = multiply(&full, &walked).unwrap();582                    got.push(walked[half][half]);583                }584                assert_eq!(got, block, "base {base} dimension {dimension}");585            }586        }587        assert_eq!(588            ladder(3, 4, 6).unwrap(),589            vec![1, 6, 132, 1848, 29040, 441408, 6772128]590        );591        assert_eq!(ladder(3, 5, 4).unwrap(), vec![1, 30, 1000, 35700, 1321600]);592        assert_eq!(593            ladder(3, 6, 4).unwrap(),594            vec![1, 20, 4030, 242300, 24642700]595        );596        assert_eq!(ladder(5, 3, 4).unwrap(), vec![1, 18, 414, 9702, 227646]);597    }598599    #[test]600    fn the_ladder_stops_where_the_exact_integers_do() {601        let terms = ladder(3, 14, 40).unwrap();602        assert_eq!(&terms[..3], &[1i128, 3432, 922_926_862]);603        assert_eq!(terms.len(), 9);604        assert_eq!(ladder(3, 10, 40).unwrap().len(), 12);605        assert_eq!(606            ladder(3, 2, 8).unwrap(),607            vec![1, 2, 4, 8, 16, 32, 64, 128, 256]608        );609    }610611    #[test]612    fn the_sign_alternates_to_the_cap_at_both_bases() {613        for base in [3usize, 5] {614            let top = cap(base).unwrap();615            for dimension in 2..=top {616                let got = sign(base, dimension).unwrap();617                let want = if dimension % 2 == 0 { -1 } else { 1 };618                assert_eq!(got, want, "base {base} dimension {dimension}");619                let block = even_block(base, dimension).unwrap();620                let root = perron(&block).unwrap();621                let edge = fill(base, dimension).unwrap() as f64 / base as f64;622                assert_eq!(623                    got,624                    if root > edge { 1 } else { -1 },625                    "base {base} dimension {dimension}"626                );627            }628            assert!(sign(base, top + 1).is_err());629        }630        assert_eq!(cap(3).unwrap(), 15);631        assert_eq!(cap(5).unwrap(), 11);632    }633634    #[test]635    fn the_spectral_ratio_falls_to_the_free_bound() {636        assert_eq!(spectral_ratio(3, 2).unwrap(), None);637        let mut last = f64::INFINITY;638        for dimension in [6usize, 10, 20, 30, 50] {639            let got = spectral_ratio(3, dimension).unwrap().unwrap();640            let free = (dimension as f64 + 2.0) / (dimension as f64 - 2.0);641            assert!(got > 1.0 && got < last, "dimension {dimension} ratio {got}");642            assert!(643                (got - free).abs() < 0.05,644                "dimension {dimension} ratio {got}"645            );646            last = got;647        }648        let fifty = spectral_ratio(3, 50).unwrap().unwrap();649        assert!((fifty - 13.0 / 12.0).abs() < 1e-9, "ratio {fifty}");650    }651}