1use crate::{linalg, Mat};
2
3#[must_use]
4pub fn fit(x: &[f64], y: &[f64], deg: usize) -> Option<Vec<f64>> {
5 let n = x.len().min(y.len());
6 let (lo, hi) = x[..n].iter().fold((f64::MAX, f64::MIN), |(l, h), &v| (l.min(v), h.max(v)));
7 let (mid, half) = ((lo + hi) / 2.0, ((hi - lo) / 2.0).max(f64::MIN_POSITIVE));
8 let mut a = Mat::zeros(n, deg + 1);
9 for i in 0..n {
10 let t = (x[i] - mid) / half;
11 let mut p = 1.0;
12 for d in 0..=deg {
13 a[(i, d)] = p;
14 p *= t;
15 }
16 }
17 let c = linalg::lstsq(&a, &y[..n], 1e-12)?;
18 let mut out = vec![0.0; deg + 1];
19 let mut basis = vec![1.0];
20 for &cd in &c {
21 for (o, &b) in out.iter_mut().zip(&basis) {
22 *o += cd * b;
23 }
24 let mut next = vec![0.0; basis.len() + 1];
25 for (k, &b) in basis.iter().enumerate() {
26 next[k + 1] += b / half;
27 next[k] -= b * mid / half;
28 }
29 basis = next;
30 }
31 Some(out)
32}
33
34#[must_use]
35pub fn eval(c: &[f64], x: f64) -> f64 {
36 c.iter().rev().fold(0.0, |acc, &k| acc * x + k)
37}
38
39#[must_use]
40pub fn integral(c: &[f64]) -> Vec<f64> {
41 std::iter::once(0.0).chain(c.iter().enumerate().map(|(d, &k)| k / (d + 1) as f64)).collect()
42}
43
44#[must_use]
45pub fn integrate(c: &[f64], a: f64, b: f64) -> f64 {
46 let p = integral(c);
47 eval(&p, b) - eval(&p, a)
48}
49
50#[cfg(test)]
51mod tests {
52 use super::*;
53
54 #[test]
55 fn fits_and_integrates_a_cubic() {
56 let want = [2.0, -1.0, 0.5, 0.25];
57 let x: Vec<f64> = (0..20).map(|i| 30.0 + f64::from(i) * 0.7).collect();
58 let y: Vec<f64> = x.iter().map(|&v| eval(&want, v)).collect();
59 let c = fit(&x, &y, 3).unwrap();
60 for &v in &x {
61 assert!((eval(&c, v) - eval(&want, v)).abs() < 1e-6 * eval(&want, v).abs().max(1.0));
62 }
63 assert!((integrate(&[0.0, 0.0, 3.0], 0.0, 2.0) - 8.0).abs() < 1e-12);
64 }
65}