slain-num 0.4.0

use crate::{linalg, Mat};

#[must_use]
pub fn fit(x: &[f64], y: &[f64], deg: usize) -> Option<Vec<f64>> {
    let n = x.len().min(y.len());
    let (lo, hi) = x[..n].iter().fold((f64::MAX, f64::MIN), |(l, h), &v| (l.min(v), h.max(v)));
    let (mid, half) = ((lo + hi) / 2.0, ((hi - lo) / 2.0).max(f64::MIN_POSITIVE));
    let mut a = Mat::zeros(n, deg + 1);
    for i in 0..n {
        let t = (x[i] - mid) / half;
        let mut p = 1.0;
        for d in 0..=deg {
            a[(i, d)] = p;
            p *= t;
        }
    }
    let c = linalg::lstsq(&a, &y[..n], 1e-12)?;
    let mut out = vec![0.0; deg + 1];
    let mut basis = vec![1.0];
    for &cd in &c {
        for (o, &b) in out.iter_mut().zip(&basis) {
            *o += cd * b;
        }
        let mut next = vec![0.0; basis.len() + 1];
        for (k, &b) in basis.iter().enumerate() {
            next[k + 1] += b / half;
            next[k] -= b * mid / half;
        }
        basis = next;
    }
    Some(out)
}

#[must_use]
pub fn eval(c: &[f64], x: f64) -> f64 {
    c.iter().rev().fold(0.0, |acc, &k| acc * x + k)
}

#[must_use]
pub fn integral(c: &[f64]) -> Vec<f64> {
    std::iter::once(0.0).chain(c.iter().enumerate().map(|(d, &k)| k / (d + 1) as f64)).collect()
}

#[must_use]
pub fn integrate(c: &[f64], a: f64, b: f64) -> f64 {
    let p = integral(c);
    eval(&p, b) - eval(&p, a)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn fits_and_integrates_a_cubic() {
        let want = [2.0, -1.0, 0.5, 0.25];
        let x: Vec<f64> = (0..20).map(|i| 30.0 + f64::from(i) * 0.7).collect();
        let y: Vec<f64> = x.iter().map(|&v| eval(&want, v)).collect();
        let c = fit(&x, &y, 3).unwrap();
        for &v in &x {
            assert!((eval(&c, v) - eval(&want, v)).abs() < 1e-6 * eval(&want, v).abs().max(1.0));
        }
        assert!((integrate(&[0.0, 0.0, 3.0], 0.0, 2.0) - 8.0).abs() < 1e-12);
    }
}