ladder.rs
23.4 kB · rust · 708 lines
1use crate::core::error::{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!(mrlyrs::math::counts::ladder::digit_polynomial(3, 3).unwrap(), vec![1, 3, 3, 6, 3, 3, 1]);67/// ```68///69/// # Errors70///71/// Errors at an even base under three, or a dimension outside two to the ceiling.72pub fn digit_polynomial(base: usize, dimension: usize) -> Result<Vec<i128>> {73 let centre = middle(base)?;74 sized(dimension)?;75 let alpha: Vec<i128> = (0..base).map(|digit| i128::from(digit != centre)).collect();76 let mut poly = vec![1i128];77 for _ in 1..dimension {78 poly = convolve(&poly, &alpha)?;79 }80 let mut second = alpha;81 second[centre] += dimension as i128;82 convolve(&poly, &second)83}8485/// The count of level-one cells the design keeps, `f_D = (q - 1)^(D-1) (q - 1 + D)`.86///87/// ```88/// assert_eq!(mrlyrs::math::counts::ladder::fill(3, 3).unwrap(), 20);89/// assert_eq!(mrlyrs::math::counts::ladder::fill(5, 3).unwrap(), 112);90/// ```91///92/// # Errors93///94/// Errors at an even base under three, a dimension outside two to the ceiling, or a product past an i128.95pub fn fill(base: usize, dimension: usize) -> Result<i128> {96 middle(base)?;97 sized(dimension)?;98 let mut out = 1i128;99 for _ in 1..dimension {100 out = product(out, base as i128 - 1)?;101 }102 product(out, (base - 1 + dimension) as i128)103}104105// MATRIX106107fn coefficient(poly: &[i128], index: i128) -> i128 {108 if index < 0 || index as usize >= poly.len() {109 0110 } else {111 poly[index as usize]112 }113}114115/// The carry matrix over the reachable carries `|c| <= (D-1)/2`, rows indexed by the carry out.116///117/// A level of the design adds one base-`q` digit per coordinate, so the height of the central118/// diagonal hyperplane moves by a digit sum `s` and a carry `c -> (c + mD - s)/q` with `m` the119/// middle digit; the map contracts onto this window from every start.120///121/// # Errors122///123/// Errors at an even base under three, or a dimension outside two to the ceiling.124pub fn carry_matrix(base: usize, dimension: usize) -> Result<Vec<Vec<i128>>> {125 let poly = digit_polynomial(base, dimension)?;126 let shift = (middle(base)? * dimension) as i128;127 let q = base as i128;128 let half = ((dimension - 1) / 2) as i128;129 Ok((-half..=half)130 .map(|out| {131 (-half..=half)132 .map(|inside| coefficient(&poly, inside + shift - q * out))133 .collect()134 })135 .collect())136}137138/// The reflection-even block of the carry matrix, of size `ceil(D/2)`.139///140/// The digit polynomial is palindromic, so the reflection `c -> -c` commutes with the carry map141/// and splits it; the even half carries the central count and its characteristic polynomial is142/// the recurrence the counts obey.143///144/// ```145/// assert_eq!(mrlyrs::math::counts::ladder::even_block(3, 3).unwrap(), vec![vec![6, 6], vec![1, 3]]);146/// ```147///148/// # Errors149///150/// Errors at an even base under three, or a dimension outside two to the ceiling.151pub fn even_block(base: usize, dimension: usize) -> Result<Vec<Vec<i128>>> {152 let poly = digit_polynomial(base, dimension)?;153 let shift = (middle(base)? * dimension) as i128;154 let q = base as i128;155 let width = (dimension - 1) / 2 + 1;156 Ok((0..width)157 .map(|out| {158 (0..width)159 .map(|inside| {160 let step = shift - q * out as i128;161 let folded = if inside == 0 {162 0163 } else {164 coefficient(&poly, step - inside as i128)165 };166 coefficient(&poly, step + inside as i128) + folded167 })168 .collect()169 })170 .collect())171}172173/// The trace of a square integer matrix.174///175/// ```176/// let block = mrlyrs::math::counts::ladder::even_block(3, 3).unwrap();177/// assert_eq!(mrlyrs::math::counts::ladder::trace(&block), 9);178/// ```179pub fn trace(rows: &[Vec<i128>]) -> i128 {180 (0..rows.len()).map(|index| rows[index][index]).sum()181}182183fn multiply(left: &[Vec<i128>], right: &[Vec<i128>]) -> Result<Vec<Vec<i128>>> {184 let n = left.len();185 let mut out = vec![vec![0i128; n]; n];186 for row in 0..n {187 for step in 0..n {188 let weight = left[row][step];189 if weight == 0 {190 continue;191 }192 for col in 0..n {193 out[row][col] = total(out[row][col], product(weight, right[step][col])?)?;194 }195 }196 }197 Ok(out)198}199200/// The monic characteristic polynomial of a square integer matrix, highest power first.201///202/// Faddeev-LeVerrier in exact integers: every division lands whole, so no fraction is ever needed203/// and the answer is the recurrence's coefficients up to sign.204///205/// ```206/// let block = mrlyrs::math::counts::ladder::even_block(3, 3).unwrap();207/// assert_eq!(mrlyrs::math::counts::ladder::characteristic(&block).unwrap(), vec![1, -9, 12]);208/// ```209///210/// # Errors211///212/// Errors when the walk leaves a remainder, or a coefficient passes an i128.213pub fn characteristic(rows: &[Vec<i128>]) -> Result<Vec<i128>> {214 let n = rows.len();215 let mut held: Vec<Vec<i128>> = (0..n)216 .map(|row| (0..n).map(|col| i128::from(row == col)).collect())217 .collect();218 let mut out = vec![1i128];219 for step in 1..=n {220 let walked = multiply(rows, &held)?;221 let sum = trace(&walked);222 if sum % step as i128 != 0 {223 return value_error("the characteristic polynomial left a remainder.");224 }225 let mark = -sum / step as i128;226 out.push(mark);227 held = walked;228 for (index, row) in held.iter_mut().enumerate().take(n) {229 row[index] = total(row[index], mark)?;230 }231 }232 Ok(out)233}234235/// The determinant of a square integer matrix, read off its characteristic polynomial.236///237/// ```238/// let block = mrlyrs::math::counts::ladder::even_block(3, 3).unwrap();239/// assert_eq!(mrlyrs::math::counts::ladder::determinant(&block).unwrap(), 12);240/// ```241///242/// # Errors243///244/// Errors when the walk leaves a remainder, or a coefficient passes an i128.245pub fn determinant(rows: &[Vec<i128>]) -> Result<i128> {246 let poly = characteristic(rows)?;247 let last = poly[rows.len()];248 Ok(if rows.len().is_multiple_of(2) {249 last250 } else {251 -last252 })253}254255// LADDER256257/// The counts `a_D(L)` of level-`L` cells meeting the central diagonal hyperplane, from `L = 0`.258///259/// The count is the top-left entry of the `L`-th power of the even block. The walk stops early260/// when the next count would pass a hundred and twenty-eight bits, so the answer runs as far as261/// exact integers reach and no further.262///263/// ```264/// let terms = mrlyrs::math::counts::ladder::ladder(3, 3, 6).unwrap();265/// assert_eq!(terms, vec![1, 6, 42, 306, 2250, 16578, 122202]);266/// ```267///268/// # Errors269///270/// Errors at an even base under three, or a dimension outside two to the ceiling.271pub fn ladder(base: usize, dimension: usize, levels: usize) -> Result<Vec<i128>> {272 let block = even_block(base, dimension)?;273 let n = block.len();274 let mut carried = vec![0i128; n];275 carried[0] = 1;276 let mut out = vec![1i128];277 for _ in 0..levels {278 let mut next = vec![0i128; n];279 for row in 0..n {280 let mut sum = 0i128;281 for col in 0..n {282 match block[row][col]283 .checked_mul(carried[col])284 .and_then(|part| sum.checked_add(part))285 {286 Some(value) => sum = value,287 None => return Ok(out),288 }289 }290 next[row] = sum;291 }292 out.push(next[0]);293 carried = next;294 }295 Ok(out)296}297298// ROOT299300/// The Perron root of a nonnegative square integer matrix.301///302/// The characteristic polynomial is rescaled by the largest row sum, which bounds the root above,303/// then the largest sign change on the unit interval is hunted on a grid and closed by bisection,304/// the same walk the recurrence growth reader takes.305///306/// ```307/// let block = mrlyrs::math::counts::ladder::even_block(3, 3).unwrap();308/// let root = mrlyrs::math::counts::ladder::perron(&block).unwrap();309/// assert!((root - (9.0 + 33f64.sqrt()) / 2.0).abs() < 1e-9);310/// ```311///312/// # Errors313///314/// Errors when the block has no positive row sum, or no root inside its bound.315pub fn perron(rows: &[Vec<i128>]) -> Result<f64> {316 let poly = characteristic(rows)?;317 let bound = rows318 .iter()319 .map(|row| row.iter().sum::<i128>())320 .max()321 .unwrap_or(0) as f64;322 if bound <= 0.0 {323 return value_error("the block has no positive row sum.");324 }325 let mut scaled = Vec::with_capacity(poly.len());326 let mut power = 1.0f64;327 for &term in &poly {328 scaled.push(term as f64 / power);329 power *= bound;330 }331 let value = |y: f64| scaled.iter().fold(0.0f64, |acc, &term| acc * y + term);332 for step in (0..SCAN).rev() {333 let mut lo = step as f64 / SCAN as f64;334 let mut hi = (step + 1) as f64 / SCAN as f64;335 if value(lo) > 0.0 || value(hi) < 0.0 {336 continue;337 }338 for _ in 0..HALVINGS {339 let mid = (lo + hi) / 2.0;340 if value(mid) <= 0.0 {341 lo = mid;342 } else {343 hi = mid;344 }345 }346 return Ok(bound * (lo + hi) / 2.0);347 }348 value_error("the block shows no root inside its row-sum bound.")349}350351// SIGN352353/// The sign of `log_q rho_D - (log_q f_D - 1)`, the slice sign law's reading, in exact integers.354///355/// The characteristic polynomial is evaluated at `f_D/q` with the denominators cleared, so the356/// answer is a comparison of whole numbers and never a rounding. The law reads `(-1)^(D+1)`: the357/// slice exponent stands above the solid's dimension less one at odd `D` and below it at even `D`.358///359/// ```360/// assert_eq!(mrlyrs::math::counts::ladder::sign(3, 3).unwrap(), 1);361/// assert_eq!(mrlyrs::math::counts::ladder::sign(3, 4).unwrap(), -1);362/// ```363///364/// # Errors365///366/// Errors at an even base under three, a dimension outside two to the ceiling, or a term past an i128.367pub fn sign(base: usize, dimension: usize) -> Result<i32> {368 let block = even_block(base, dimension)?;369 let poly = characteristic(&block)?;370 let full = fill(base, dimension)?;371 let q = base as i128;372 let mut acc = 0i128;373 let mut weight = 1i128;374 for &term in &poly {375 acc = total(product(acc, full)?, product(term, weight)?)?;376 weight = product(weight, q)?;377 }378 Ok(match acc {379 acc if acc > 0 => -1,380 acc if acc < 0 => 1,381 _ => 0,382 })383}384385/// The widest dimension the exact carry arithmetic reaches at the base.386///387/// The sign is the heaviest reading: it walks the characteristic polynomial at `f_D/q` with the388/// denominators cleared, so its running value climbs like `f_D^(D/2)` and a hundred and twenty-eight389/// bits run out at dimension fifteen in base three and dimension eleven in base five.390///391/// ```392/// assert_eq!(mrlyrs::math::counts::ladder::cap(3).unwrap(), 15);393/// assert_eq!(mrlyrs::math::counts::ladder::cap(5).unwrap(), 11);394/// ```395///396/// # Errors397///398/// Errors at an even base under three, or past the dimension the ladder reaches.399pub fn cap(base: usize) -> Result<usize> {400 middle(base)?;401 let mut top = 0;402 for dimension in 2..=CEILING {403 if sign(base, dimension).is_err() {404 break;405 }406 top = dimension;407 }408 if top < 2 {409 return value_error(format!(410 "base {base} reaches no dimension in exact integers."411 ));412 }413 Ok(top)414}415416// SPECTRUM417418fn ahead(walk: &[Vec<f64>], vector: &[f64]) -> Vec<f64> {419 walk.iter()420 .map(|row| row.iter().zip(vector).map(|(&a, &x)| a * x).sum())421 .collect()422}423424fn behind(walk: &[Vec<f64>], vector: &[f64]) -> Vec<f64> {425 (0..vector.len())426 .map(|col| {427 (0..vector.len())428 .map(|row| walk[row][col] * vector[row])429 .sum()430 })431 .collect()432}433434fn top(vector: &[f64]) -> f64 {435 vector.iter().fold(0.0f64, |peak, &x| peak.max(x.abs()))436}437438fn settle(walk: &[Vec<f64>], left: bool) -> Option<Vec<f64>> {439 let mut vector = vec![1.0f64; walk.len()];440 for _ in 0..SWEEPS {441 let next = if left {442 behind(walk, &vector)443 } else {444 ahead(walk, &vector)445 };446 let peak = top(&next);447 if peak == 0.0 {448 return None;449 }450 vector = next.iter().map(|&x| x / peak).collect();451 }452 Some(vector)453}454455/// The Perron root over the modulus of the second eigenvalue, or none where the block is one wide.456///457/// Power iteration finds the leading pair, Hotelling deflation takes it out and a second walk,458/// reprojected every sweep so rounding cannot bring the leader back, reads the runner up. At base459/// three the ratio falls to `(D + 2)/(D - 2)`, so no fixed spectral gap survives the dimensions.460///461/// # Errors462///463/// Errors at an even base under three, or a dimension outside two to the ceiling.464pub fn spectral_ratio(base: usize, dimension: usize) -> Result<Option<f64>> {465 let block = even_block(base, dimension)?;466 let n = block.len();467 if n < 2 {468 return Ok(None);469 }470 let bound = block471 .iter()472 .map(|row| row.iter().sum::<i128>())473 .max()474 .unwrap_or(1) as f64;475 let walk: Vec<Vec<f64>> = block476 .iter()477 .map(|row| row.iter().map(|&x| x as f64 / bound).collect())478 .collect();479 let (right, left) = match (settle(&walk, false), settle(&walk, true)) {480 (Some(right), Some(left)) => (right, left),481 _ => return Ok(None),482 };483 let lead = top(&ahead(&walk, &right));484 let pair: f64 = left.iter().zip(&right).map(|(&u, &v)| u * v).sum();485 if pair == 0.0 || lead == 0.0 {486 return Ok(None);487 }488 let mut trail: Vec<f64> = (0..n).map(|i| 1.0 + 0.37 * ((i * 7) % 5) as f64).collect();489 let mut second = 0.0f64;490 for _ in 0..SWEEPS {491 let share: f64 = left.iter().zip(&trail).map(|(&u, &x)| u * x).sum::<f64>() / pair;492 let cleared: Vec<f64> = trail493 .iter()494 .zip(&right)495 .map(|(&x, &v)| x - share * v)496 .collect();497 let next = ahead(&walk, &cleared);498 second = top(&next);499 if second == 0.0 {500 return Ok(None);501 }502 trail = next.iter().map(|&x| x / second).collect();503 }504 Ok(Some(lead / second))505}506507#[cfg(test)]508mod tests {509 use super::*;510 use crate::math::bang::Code;511512 fn brute(base: usize, dimension: usize) -> Vec<i128> {513 let centre = (base - 1) / 2;514 let mut out = vec![0i128; (base - 1) * dimension + 1];515 for mut code in 0..base.pow(dimension as u32) {516 let mut sum = 0;517 let mut middles = 0;518 for _ in 0..dimension {519 let digit = code % base;520 code /= base;521 sum += digit;522 middles += usize::from(digit == centre);523 }524 if middles <= 1 {525 out[sum] += 1;526 }527 }528 out529 }530531 fn central(base: usize, dimension: usize, levels: usize) -> Vec<i128> {532 let poly = digit_polynomial(base, dimension).unwrap();533 let mut out = vec![1i128];534 let mut span = vec![1i128];535 let mut step = 1usize;536 for _ in 0..levels {537 let mut next = vec![0i128; span.len() + (base - 1) * dimension * step];538 for (at, &count) in span.iter().enumerate() {539 if count == 0 {540 continue;541 }542 for (weight, &many) in poly.iter().enumerate() {543 next[at + weight * step] += count * many;544 }545 }546 span = next;547 step *= base;548 out.push(span[dimension * (step - 1) / 2]);549 }550 out551 }552553 #[test]554 fn the_digit_polynomial_is_the_enumeration_of_the_kept_cells() {555 for base in [3usize, 5] {556 for dimension in 2..=6 {557 let got = digit_polynomial(base, dimension).unwrap();558 assert_eq!(559 got,560 brute(base, dimension),561 "base {base} dimension {dimension}"562 );563 assert_eq!(564 got.iter().sum::<i128>(),565 fill(base, dimension).unwrap(),566 "base {base} dimension {dimension}"567 );568 }569 }570 assert!(digit_polynomial(4, 3).is_err());571 assert!(digit_polynomial(3, 1).is_err());572 }573574 #[test]575 fn the_three_dimensional_block_is_the_published_anchor() {576 let block = even_block(3, 3).unwrap();577 assert_eq!(block, vec![vec![6, 6], vec![1, 3]]);578 assert_eq!(trace(&block), 9);579 assert_eq!(determinant(&block).unwrap(), 12);580 assert_eq!(characteristic(&block).unwrap(), vec![1, -9, 12]);581 assert_eq!(582 ladder(3, 3, 6).unwrap(),583 vec![1, 6, 42, 306, 2250, 16578, 122202]584 );585 let root = perron(&block).unwrap();586 assert!((root - (9.0 + 33f64.sqrt()) / 2.0).abs() < 1e-9);587 assert!((root.log(3.0) - 1.818_410).abs() < 1e-6);588 assert!((20f64.log(3.0) - 1.0 - 1.726_833).abs() < 1e-6);589 }590591 #[test]592 fn the_traces_read_the_two_closed_forms() {593 let want = [2i128, 9, 11, 60, 47, 336];594 for (index, &term) in want.iter().enumerate() {595 let dimension = index + 2;596 assert_eq!(trace(&even_block(3, dimension).unwrap()), term);597 }598 for dimension in 2..=cap(3).unwrap() {599 let got = trace(&even_block(3, dimension).unwrap());600 let want = if dimension % 2 == 0 {601 3 * 2i128.pow(dimension as u32 - 2) - 1602 } else {603 3 * dimension as i128 * 2i128.pow(dimension as u32 - 3)604 };605 assert_eq!(got, want, "dimension {dimension}");606 }607 }608609 #[test]610 fn the_three_generators_of_the_ladder_agree() {611 for base in [3usize, 5] {612 for dimension in 2..=6 {613 let block = ladder(base, dimension, 4).unwrap();614 assert_eq!(615 block,616 central(base, dimension, 4),617 "base {base} dimension {dimension}"618 );619 let full = carry_matrix(base, dimension).unwrap();620 let half = full.len() / 2;621 let mut walked: Vec<Vec<i128>> = (0..full.len())622 .map(|row| (0..full.len()).map(|col| i128::from(row == col)).collect())623 .collect();624 let mut got = vec![1i128];625 for _ in 0..4 {626 walked = multiply(&full, &walked).unwrap();627 got.push(walked[half][half]);628 }629 assert_eq!(got, block, "base {base} dimension {dimension}");630 }631 }632 assert_eq!(633 ladder(3, 4, 6).unwrap(),634 vec![1, 6, 132, 1848, 29040, 441408, 6772128]635 );636 assert_eq!(ladder(3, 5, 4).unwrap(), vec![1, 30, 1000, 35700, 1321600]);637 assert_eq!(638 ladder(3, 6, 4).unwrap(),639 vec![1, 20, 4030, 242300, 24642700]640 );641 assert_eq!(ladder(5, 3, 4).unwrap(), vec![1, 18, 414, 9702, 227646]);642 }643644 #[test]645 fn the_ladder_stops_where_the_exact_integers_do() {646 let terms = ladder(3, 14, 40).unwrap();647 assert_eq!(&terms[..3], &[1i128, 3432, 922_926_862]);648 assert_eq!(terms.len(), 9);649 assert_eq!(ladder(3, 10, 40).unwrap().len(), 12);650 assert_eq!(651 ladder(3, 2, 8).unwrap(),652 vec![1, 2, 4, 8, 16, 32, 64, 128, 256]653 );654 }655656 #[test]657 fn the_sign_alternates_to_the_cap_at_both_bases() {658 for base in [3usize, 5] {659 let top = cap(base).unwrap();660 for dimension in 2..=top {661 let got = sign(base, dimension).unwrap();662 let want = if dimension % 2 == 0 { -1 } else { 1 };663 assert_eq!(got, want, "base {base} dimension {dimension}");664 let block = even_block(base, dimension).unwrap();665 let root = perron(&block).unwrap();666 let edge = fill(base, dimension).unwrap() as f64 / base as f64;667 assert_eq!(668 got,669 if root > edge { 1 } else { -1 },670 "base {base} dimension {dimension}"671 );672 }673 assert!(sign(base, top + 1).is_err());674 }675 assert_eq!(cap(3).unwrap(), 15);676 assert_eq!(cap(5).unwrap(), 11);677 }678679 #[test]680 fn the_spectral_ratio_falls_to_the_free_bound() {681 assert_eq!(spectral_ratio(3, 2).unwrap(), None);682 let mut last = f64::INFINITY;683 for dimension in [6usize, 10, 20, 30, 50] {684 let got = spectral_ratio(3, dimension).unwrap().unwrap();685 let free = (dimension as f64 + 2.0) / (dimension as f64 - 2.0);686 assert!(got > 1.0 && got < last, "dimension {dimension} ratio {got}");687 assert!(688 (got - free).abs() < 0.05,689 "dimension {dimension} ratio {got}"690 );691 last = got;692 }693 let fifty = spectral_ratio(3, 50).unwrap().unwrap();694 assert!((fifty - 13.0 / 12.0).abs() < 1e-9, "ratio {fifty}");695 }696 #[test]697 fn the_ladder_is_the_sponge_diagonal_count() {698 let terms = ladder(3, 3, 5).unwrap();699 assert_eq!(terms, vec![1, 6, 42, 306, 2250, 16578]);700 let counted: Vec<u128> = (1..=5)701 .map(|level| {702 let height = 3 * (3usize.pow(level as u32) - 1) / 2;703 crate::math::three::profile(Code::from(23u64), 3, level, 2).unwrap()[height]704 })705 .collect();706 assert_eq!(counted, vec![6, 42, 306, 2250, 16578]);707 }708}