quasi.rs

1use mrlymath::bang::Code;2use mrlymath::formulas::fill;3use num_bigint::{BigInt, Sign};4use num_rational::BigRational;56pub type Poly = Vec<BigRational>;78const EXTRA: usize = 5;910fn rational(value: i128) -> BigRational {11    BigRational::from_integer(BigInt::from(value))12}1314fn is_zero(value: &BigRational) -> bool {15    value.numer().sign() == Sign::NoSign16}1718pub fn trim(mut poly: Poly) -> Poly {19    while poly.last().is_some_and(is_zero) {20        poly.pop();21    }22    poly23}2425fn add(a: &[BigRational], b: &[BigRational]) -> Poly {26    let zero = rational(0);27    let width = a.len().max(b.len());28    (0..width)29        .map(|i| a.get(i).unwrap_or(&zero) + b.get(i).unwrap_or(&zero))30        .collect()31}3233fn multiply(a: &[BigRational], b: &[BigRational]) -> Poly {34    let mut out = vec![rational(0); a.len() + b.len() - 1];35    for (i, x) in a.iter().enumerate() {36        for (j, y) in b.iter().enumerate() {37            out[i + j] += x * y;38        }39    }40    out41}4243pub fn evaluate(poly: &[BigRational], x: i128) -> BigRational {44    let point = rational(x);45    poly.iter()46        .rev()47        .fold(rational(0), |acc, coefficient| acc * &point + coefficient)48}4950pub fn interpolate(points: &[(i128, u128)]) -> Poly {51    let mut out = Poly::new();52    for (i, &(xi, yi)) in points.iter().enumerate() {53        let mut term = vec![rational(yi as i128)];54        for (j, &(xj, _)) in points.iter().enumerate() {55            if j != i {56                let scale = rational(xi - xj);57                term = multiply(&term, &[-rational(xj) / &scale, rational(1) / &scale]);58            }59        }60        out = add(&out, &term);61    }62    trim(out)63}6465pub fn degree(poly: &[BigRational]) -> Option<usize> {66    trim(poly.to_vec()).len().checked_sub(1)67}6869pub fn leading(poly: &[BigRational]) -> BigRational {70    trim(poly.to_vec())71        .last()72        .cloned()73        .unwrap_or_else(|| rational(0))74}7576pub fn fit(code: Code, dimension: usize, base: usize) -> (Vec<Poly>, bool) {77    let polys: Vec<Poly> = (0..base)78        .map(|class| {79            let points: Vec<(i128, u128)> = (0..)80                .map(|j| class + j * base)81                .filter(|&n| n >= 1)82                .take(dimension + 1 + EXTRA)83                .map(|n| {84                    let value = fill(code, n, dimension, 1, base).expect("the fill counts");85                    (n as i128, value)86                })87                .collect();88            let poly = interpolate(&points[..=dimension]);89            for &(x, y) in &points {90                assert!(91                    evaluate(&poly, x) == rational(y as i128),92                    "the class polynomial misses the fill"93                );94            }95            poly96        })97        .collect();98    let collapses = polys.iter().all(|poly| *poly == polys[0]);99    (polys, collapses)100}101102pub fn fraction(value: &BigRational) -> String {103    if value.denom() == &BigInt::from(1) {104        value.numer().to_string()105    } else {106        format!("{}/{}", value.numer(), value.denom())107    }108}109110fn term(magnitude: &BigRational, degree: usize) -> String {111    if degree == 0 {112        return fraction(magnitude);113    }114    let monomial = if degree == 1 {115        String::from("n")116    } else {117        format!("n**{degree}")118    };119    let one = BigInt::from(1);120    match (magnitude.numer() == &one, magnitude.denom() == &one) {121        (true, true) => monomial,122        (true, false) => format!("{monomial}/{}", magnitude.denom()),123        (false, true) => format!("{}*{monomial}", magnitude.numer()),124        (false, false) => format!("{}*{monomial}/{}", magnitude.numer(), magnitude.denom()),125    }126}127128pub fn text(poly: &[BigRational]) -> String {129    let trimmed = trim(poly.to_vec());130    if trimmed.is_empty() {131        return String::from("0");132    }133    let mut out = String::new();134    for degree in (0..trimmed.len()).rev() {135        let coefficient = &trimmed[degree];136        if is_zero(coefficient) {137            continue;138        }139        let negative = coefficient.numer().sign() == Sign::Minus;140        let magnitude = if negative {141            -coefficient.clone()142        } else {143            coefficient.clone()144        };145        if out.is_empty() {146            if negative {147                out.push('-');148            }149        } else {150            out.push_str(if negative { " - " } else { " + " });151        }152        out.push_str(&term(&magnitude, degree));153    }154    out155}