blend.rs

1const PRIMES: [u64; 3] = [2_147_483_647, 2_147_483_629, 2_147_483_587];23/// Adds two sequences term by term over their shared length.4pub fn add(a: &[i128], b: &[i128]) -> Vec<i128> {5    a.iter().zip(b).map(|(&x, &y)| x + y).collect()6}78/// Subtracts the second sequence from the first over their shared length.9pub fn sub(a: &[i128], b: &[i128]) -> Vec<i128> {10    a.iter().zip(b).map(|(&x, &y)| x - y).collect()11}1213/// Multiplies two sequences term by term over their shared length.14pub fn hadamard(a: &[i128], b: &[i128]) -> Vec<i128> {15    a.iter().zip(b).map(|(&x, &y)| x * y).collect()16}1718/// Convolves two sequences, keeping the exact prefix their shared length affords.19///20/// Panics when a convolution sum passes a signed hundred and twenty-eight bits.21pub fn cauchy(a: &[i128], b: &[i128]) -> Vec<i128> {22    let length = a.len().min(b.len());23    (0..length)24        .map(|n| {25            (0..=n).fold(0i128, |sum, i| {26                a[i].checked_mul(b[n - i])27                    .and_then(|v| sum.checked_add(v))28                    .expect("the convolution passes a hundred and twenty-eight bits")29            })30        })31        .collect()32}3334/// Drops the first terms of a sequence.35pub fn shift(a: &[i128], count: usize) -> Vec<i128> {36    a.iter().skip(count).copied().collect()37}3839/// Keeps every step-th term from the offset onward.40///41/// Panics at a step of zero.42pub fn decimate(a: &[i128], step: usize, offset: usize) -> Vec<i128> {43    assert!(step > 0, "decimate needs a step above zero");44    a.iter().skip(offset).step_by(step).copied().collect()45}4647/// Returns the first differences of a sequence, one term shorter.48pub fn delta(a: &[i128]) -> Vec<i128> {49    a.windows(2).map(|w| w[1] - w[0]).collect()50}5152/// Returns the partial sums of a sequence.53///54/// Panics when a partial sum passes a signed hundred and twenty-eight bits.55pub fn sigma(a: &[i128]) -> Vec<i128> {56    let mut total = 0i128;57    a.iter()58        .map(|&x| {59            total = total60                .checked_add(x)61                .expect("the partial sums pass a hundred and twenty-eight bits");62            total63        })64        .collect()65}6667/// Multiplies every term of a sequence by the factor.68pub fn scale(a: &[i128], factor: i128) -> Vec<i128> {69    a.iter().map(|&x| x * factor).collect()70}7172fn residues(terms: &[i128], prime: u64) -> Vec<u64> {73    let p = prime as i128;74    terms.iter().map(|&t| (t.rem_euclid(p)) as u64).collect()75}7677fn inverse(value: u64, prime: u64) -> u64 {78    let mut power = prime - 2;79    let mut out = 1u64;80    let mut factor = value % prime;81    while power > 0 {82        if power & 1 == 1 {83            out = out * factor % prime;84        }85        factor = factor * factor % prime;86        power >>= 1;87    }88    out89}9091fn massey(seq: &[u64], prime: u64) -> usize {92    let mut connection = vec![1u64];93    let mut previous = vec![1u64];94    let mut order = 0usize;95    let mut gap = 1usize;96    let mut last_delta = 1u64;97    for n in 0..seq.len() {98        let mut delta = 0u64;99        for (i, &c) in connection.iter().enumerate() {100            if i <= n {101                delta = (delta + c * seq[n - i]) % prime;102            }103        }104        if delta == 0 {105            gap += 1;106            continue;107        }108        if 2 * order > n {109            let scale = delta * inverse(last_delta, prime) % prime;110            for (i, &b) in previous.iter().enumerate() {111                let slot = i + gap;112                if slot >= connection.len() {113                    connection.resize(slot + 1, 0);114                }115                connection[slot] = (connection[slot] + prime * prime - scale * b) % prime;116            }117            gap += 1;118            continue;119        }120        let held = connection.clone();121        let scale = delta * inverse(last_delta, prime) % prime;122        for (i, &b) in previous.iter().enumerate() {123            let slot = i + gap;124            if slot >= connection.len() {125                connection.resize(slot + 1, 0);126            }127            connection[slot] = (connection[slot] + prime * prime - scale * b) % prime;128        }129        previous = held;130        last_delta = delta;131        gap = 1;132        order = n + 1 - order;133    }134    order135}136137fn fit(seq: &[u64], order: usize, prime: u64) -> Option<Vec<u64>> {138    let rows = seq.len() - order;139    let mut matrix: Vec<Vec<u64>> = (0..rows)140        .map(|r| {141            let mut row: Vec<u64> = (0..order).map(|c| seq[r + order - 1 - c]).collect();142            row.push(seq[r + order]);143            row144        })145        .collect();146    let mut pivots = Vec::new();147    let mut lead = 0usize;148    for col in 0..order {149        let Some(pivot) = (lead..rows).find(|&r| matrix[r][col] != 0) else {150            continue;151        };152        matrix.swap(lead, pivot);153        let inv = inverse(matrix[lead][col], prime);154        for value in &mut matrix[lead][col..=order] {155            *value = *value * inv % prime;156        }157        let cleared = matrix[lead].clone();158        for (row, entries) in matrix.iter_mut().enumerate() {159            if row != lead && entries[col] != 0 {160                let factor = entries[col];161                for (value, &top) in entries[col..=order].iter_mut().zip(&cleared[col..=order]) {162                    *value = (*value + prime - factor * top % prime) % prime;163                }164            }165        }166        pivots.push((lead, col));167        lead += 1;168    }169    if matrix[lead..].iter().any(|row| row[order] != 0) {170        return None;171    }172    let mut out = vec![0u64; order];173    for &(row, col) in &pivots {174        out[col] = matrix[row][order];175    }176    Some(out)177}178179fn combine(parts: [u64; 3]) -> u128 {180    let (p0, p1, p2) = (PRIMES[0] as u128, PRIMES[1] as u128, PRIMES[2] as u128);181    let m01 = p0 * p1;182    let step = (parts[1] as u128 + p1 - parts[0] as u128 % p1) % p1;183    let lift = inverse((p0 % p1) as u64, PRIMES[1]) as u128;184    let x01 = parts[0] as u128 + p0 * (step * lift % p1);185    let step2 = (parts[2] as u128 + p2 - x01 % p2) % p2;186    let lift2 = inverse((m01 % p2) as u64, PRIMES[2]) as u128;187    x01 + m01 * (step2 * lift2 % p2)188}189190fn reconstruct(residue: u128) -> Option<(i128, i128)> {191    let modulus = (PRIMES[0] as u128) * (PRIMES[1] as u128) * (PRIMES[2] as u128);192    let bound: u128 = 1 << 45;193    let (mut r0, mut r1) = (modulus as i128, residue as i128);194    let (mut t0, mut t1) = (0i128, 1i128);195    while r1.unsigned_abs() > bound {196        let q = r0 / r1;197        (r0, r1) = (r1, r0 - q * r1);198        (t0, t1) = (t1, t0 - q * t1);199    }200    if t1 == 0 || t1.unsigned_abs() > bound {201        return None;202    }203    let (num, den) = if t1 < 0 { (-r1, -t1) } else { (r1, t1) };204    let (mut a, mut b) = (num.unsigned_abs(), den.unsigned_abs());205    while b != 0 {206        (a, b) = (b, a % b);207    }208    if a > 1 {209        return None;210    }211    Some((num, den))212}213214fn verify(terms: &[i128], coefficients: &[(i128, i128)]) -> bool {215    let order = coefficients.len();216    let mut clear = 1i128;217    for &(_, den) in coefficients {218        clear = clear / gcd(clear, den) * den;219    }220    let weights: Vec<i128> = coefficients221        .iter()222        .map(|&(num, den)| num * (clear / den))223        .collect();224    let mut exact = true;225    for n in order..terms.len() {226        let mut sum = Some(0i128);227        for (i, &w) in weights.iter().enumerate() {228            sum = sum.and_then(|s| {229                w.checked_mul(terms[n - 1 - i])230                    .and_then(|v| s.checked_add(v))231            });232        }233        match (sum, terms[n].checked_mul(clear)) {234            (Some(left), Some(right)) => {235                if left != right {236                    return false;237                }238            }239            _ => {240                exact = false;241                break;242            }243        }244    }245    if exact {246        return true;247    }248    PRIMES.iter().all(|&prime| {249        let seq = residues(terms, prime);250        let p = prime as u128;251        let c = (clear.rem_euclid(prime as i128)) as u128;252        let ws: Vec<u128> = weights253            .iter()254            .map(|&w| (w.rem_euclid(prime as i128)) as u128)255            .collect();256        (order..terms.len()).all(|n| {257            let sum = ws258                .iter()259                .enumerate()260                .fold(0u128, |s, (i, &w)| (s + w * seq[n - 1 - i] as u128) % p);261            sum == c * seq[n] as u128 % p262        })263    })264}265266fn gcd(a: i128, b: i128) -> i128 {267    let (mut a, mut b) = (a.abs(), b.abs());268    while b != 0 {269        (a, b) = (b, a % b);270    }271    a.max(1)272}273274/// Finds the smallest linear constant-coefficient recurrence that fits every supplied term.275///276/// The coefficients come back as reduced fractions with the newest term first, so a277/// result of one and one twice is the Fibonacci rule. The hunt runs modulo three278/// primes, rebuilds the rationals by remainder reconstruction, and verifies the rule279/// on every term before answering; too few terms to trust an order returns nothing,280/// and the zero sequence returns the empty rule.281///282/// ```283/// let fib = vec![0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89];284/// assert_eq!(mrlynum::blend::recurrence(&fib), Some(vec![(1, 1), (1, 1)]));285/// ```286pub fn recurrence(terms: &[i128]) -> Option<Vec<(i128, i128)>> {287    if terms.len() < 4 {288        return None;289    }290    if terms.iter().all(|&t| t == 0) {291        return Some(Vec::new());292    }293    let order = PRIMES294        .iter()295        .map(|&p| massey(&residues(terms, p), p))296        .max()?;297    if order == 0 || 2 * order + 2 > terms.len() {298        return None;299    }300    let mut parts = Vec::new();301    for &prime in &PRIMES {302        parts.push(fit(&residues(terms, prime), order, prime)?);303    }304    let coefficients: Option<Vec<(i128, i128)>> = (0..order)305        .map(|i| reconstruct(combine([parts[0][i], parts[1][i], parts[2][i]])))306        .collect();307    let coefficients = coefficients?;308    verify(terms, &coefficients).then_some(coefficients)309}310311/// Returns the monic characteristic polynomial of a recurrence, highest power first.312pub fn characteristic(coefficients: &[(i128, i128)]) -> Vec<(i128, i128)> {313    let mut out = vec![(1, 1)];314    out.extend(coefficients.iter().map(|&(num, den)| (-num, den)));315    out316}317318/// Returns the largest positive real root of a recurrence's characteristic polynomial, the growth rate, or a not-a-number where no real root lands.319///320/// A simple root lands at machine precision; a repeated root lands to a few decimals only.321///322/// ```323/// let rate = mrlynum::blend::growth(&[(1, 1), (1, 1)]);324/// assert!((rate - 1.618033988749895).abs() < 1e-12);325/// ```326pub fn growth(coefficients: &[(i128, i128)]) -> f64 {327    let poly: Vec<f64> = characteristic(coefficients)328        .iter()329        .map(|&(num, den)| num as f64 / den as f64)330        .collect();331    let value = |x: f64| poly.iter().fold(0.0, |acc, &c| acc * x + c);332    let bound = 1.0 + poly.iter().skip(1).map(|c| c.abs()).fold(0.0f64, f64::max);333    let steps = 8192;334    let mut best = f64::NAN;335    let mut gap = f64::INFINITY;336    for i in (0..steps).rev() {337        let (lo, hi) = (338            bound * i as f64 / steps as f64,339            bound * (i + 1) as f64 / steps as f64,340        );341        if value(lo) <= 0.0 && value(hi) >= 0.0 {342            let (mut lo, mut hi) = (lo, hi);343            for _ in 0..200 {344                let mid = (lo + hi) / 2.0;345                if value(mid) <= 0.0 {346                    lo = mid;347                } else {348                    hi = mid;349                }350            }351            return (lo + hi) / 2.0;352        }353        let mid = (lo + hi) / 2.0;354        if value(mid).abs() < gap {355            gap = value(mid).abs();356            best = mid;357        }358    }359    let mut x = best;360    for _ in 0..100 {361        let h = 1e-7 * x.abs().max(1.0);362        let slope = (value(x + h) - value(x - h)) / (2.0 * h);363        if slope == 0.0 {364            break;365        }366        x -= value(x) / slope;367    }368    if value(x).abs() < 1e-6 * (1.0 + x.abs().powi(poly.len() as i32 - 1)) {369        x370    } else {371        f64::NAN372    }373}374375#[cfg(test)]376mod tests {377    use super::*;378    use crate::classics::{catalan, primes};379380    fn fibonacci(count: usize) -> Vec<i128> {381        let mut out = vec![0i128, 1];382        while out.len() < count {383            out.push(out[out.len() - 1] + out[out.len() - 2]);384        }385        out.truncate(count);386        out387    }388389    #[test]390    fn the_fibonacci_rule_and_its_golden_growth() {391        let fib = fibonacci(30);392        let rule = recurrence(&fib).unwrap();393        assert_eq!(rule, vec![(1, 1), (1, 1)]);394        assert!((growth(&rule) - 1.618_033_988_749_895).abs() < 1e-12);395    }396397    #[test]398    fn a_cubic_polynomial_sequence_needs_order_four() {399        let terms: Vec<i128> = (0..20).map(|n| n * n * (4 * n + 3)).collect();400        let rule = recurrence(&terms).unwrap();401        assert_eq!(rule, vec![(4, 1), (-6, 1), (4, 1), (-1, 1)]);402        assert!((growth(&rule) - 1.0).abs() < 1e-3);403    }404405    #[test]406    fn the_hexagram_rule_grows_by_its_perron_root() {407        let mut terms = vec![1i128, 6];408        while terms.len() < 14 {409            let n = terms.len();410            terms.push(9 * terms[n - 1] - 12 * terms[n - 2]);411        }412        let rule = recurrence(&terms).unwrap();413        assert_eq!(rule, vec![(9, 1), (-12, 1)]);414        assert!((growth(&rule) - (9.0 + 33f64.sqrt()) / 2.0).abs() < 1e-12);415    }416417    #[test]418    fn a_halving_sequence_carries_a_rational_rule() {419        let terms = vec![1024i128, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1];420        assert_eq!(recurrence(&terms), Some(vec![(1, 2)]));421    }422423    #[test]424    fn the_primes_and_the_catalan_numbers_refuse_every_rule() {425        let ps: Vec<i128> = primes(200).iter().map(|&p| p as i128).collect();426        assert_eq!(recurrence(&ps), None);427        let cs: Vec<i128> = catalan(40_000_000).iter().map(|&c| c as i128).collect();428        assert_eq!(recurrence(&cs), None);429    }430431    #[test]432    fn blends_of_c_finite_sequences_stay_c_finite() {433        let fib = fibonacci(40);434        let squared = hadamard(&fib, &fib);435        assert_eq!(recurrence(&squared).unwrap().len(), 3);436        let summed = sigma(&fib);437        assert_eq!(recurrence(&summed).unwrap().len(), 3);438        let paired = decimate(&fib, 2, 0);439        let rule = recurrence(&paired).unwrap();440        assert_eq!(rule, vec![(3, 1), (-1, 1)]);441        let convolved = cauchy(&fib, &fib);442        assert_eq!(recurrence(&convolved).unwrap().len(), 4);443    }444445    #[test]446    fn the_term_ops_keep_their_hand_checked_values() {447        let a = vec![1i128, 2, 3, 4, 5];448        let b = vec![10i128, 20, 30];449        assert_eq!(add(&a, &b), vec![11, 22, 33]);450        assert_eq!(sub(&b, &a), vec![9, 18, 27]);451        assert_eq!(hadamard(&a, &b), vec![10, 40, 90]);452        assert_eq!(cauchy(&a, &b), vec![10, 40, 100]);453        assert_eq!(shift(&a, 2), vec![3, 4, 5]);454        assert_eq!(decimate(&a, 2, 1), vec![2, 4]);455        assert_eq!(delta(&a), vec![1, 1, 1, 1]);456        assert_eq!(sigma(&a), vec![1, 3, 6, 10, 15]);457        assert_eq!(scale(&a, -2), vec![-2, -4, -6, -8, -10]);458    }459460    #[test]461    fn the_zero_sequence_returns_the_empty_rule() {462        assert_eq!(recurrence(&[0, 0, 0, 0, 0, 0]), Some(vec![]));463        assert_eq!(recurrence(&[1, 2]), None);464    }465466    #[test]467    fn a_short_prefix_is_not_enough_evidence() {468        assert_eq!(recurrence(&[1, 2, 4]), None);469        assert_eq!(recurrence(&[1, 2, 4, 8]), Some(vec![(2, 1)]));470    }471472    #[test]473    #[should_panic(expected = "decimate needs a step above zero")]474    fn decimate_refuses_a_zero_step() {475        let _ = decimate(&[1, 2, 3], 0, 0);476    }477}