qualia_core_db/solvers/interpolation/
lagrange.rs

1//! Polynomial interpolation through `n` points: the Lagrange form (direct evaluation)
2//! and the Newton divided-difference form (build coefficients once, evaluate cheaply).
3
4use super::InterpolationError;
5
6fn validate(xs: &[f64], ys: &[f64]) -> Result<(), InterpolationError> {
7    if xs.is_empty() || xs.len() != ys.len() {
8        return Err(InterpolationError::InsufficientData);
9    }
10    for i in 0..xs.len() {
11        for j in (i + 1)..xs.len() {
12            if xs[i] == xs[j] {
13                return Err(InterpolationError::DuplicateNodes);
14            }
15        }
16    }
17    Ok(())
18}
19
20/// Evaluate the Lagrange interpolating polynomial through `(xs, ys)` at `x`.
21pub fn lagrange_eval(xs: &[f64], ys: &[f64], x: f64) -> Result<f64, InterpolationError> {
22    validate(xs, ys)?;
23    let n = xs.len();
24    let mut sum = 0.0;
25    for i in 0..n {
26        let mut li = 1.0;
27        for j in 0..n {
28            if i != j {
29                li *= (x - xs[j]) / (xs[i] - xs[j]);
30            }
31        }
32        sum += ys[i] * li;
33    }
34    Ok(sum)
35}
36
37/// Newton divided-difference coefficients for `(xs, ys)` (the leading diagonal of the
38/// divided-difference table). Use with [`newton_eval`].
39pub fn newton_coefficients(xs: &[f64], ys: &[f64]) -> Result<Vec<f64>, InterpolationError> {
40    validate(xs, ys)?;
41    let n = xs.len();
42    let mut coef = ys.to_vec();
43    for j in 1..n {
44        for i in (j..n).rev() {
45            coef[i] = (coef[i] - coef[i - 1]) / (xs[i] - xs[i - j]);
46        }
47    }
48    Ok(coef)
49}
50
51/// Evaluate the Newton form (Horner over the nested products) at `x`.
52pub fn newton_eval(xs: &[f64], coef: &[f64], x: f64) -> f64 {
53    let n = coef.len();
54    let mut acc = coef[n - 1];
55    for i in (0..n - 1).rev() {
56        acc = acc * (x - xs[i]) + coef[i];
57    }
58    acc
59}
60
61#[cfg(test)]
62mod tests {
63    use super::*;
64    const EPS: f64 = 1e-9;
65
66    #[test]
67    fn interpolant_passes_through_nodes() {
68        let xs = [0.0, 1.0, 2.0, 3.0];
69        let ys = [1.0, 3.0, 2.0, 5.0];
70        for i in 0..xs.len() {
71            assert!((lagrange_eval(&xs, &ys, xs[i]).unwrap() - ys[i]).abs() < EPS);
72        }
73    }
74
75    #[test]
76    fn reproduces_a_quadratic_exactly() {
77        // f(x) = 2x² − 3x + 1 sampled at 3 points → interpolant equals f everywhere.
78        let f = |x: f64| 2.0 * x * x - 3.0 * x + 1.0;
79        let xs = [-1.0, 0.0, 2.0];
80        let ys = xs.map(f);
81        for &q in &[0.5, 1.7, -0.3, 5.0] {
82            assert!((lagrange_eval(&xs, &ys, q).unwrap() - f(q)).abs() < 1e-8);
83        }
84    }
85
86    #[test]
87    fn newton_matches_lagrange() {
88        let xs = [0.0, 1.0, 2.0, 4.0];
89        let ys = [1.0, 2.0, 0.0, 8.0];
90        let coef = newton_coefficients(&xs, &ys).unwrap();
91        for &q in &[0.3, 1.5, 3.0] {
92            assert!(
93                (newton_eval(&xs, &coef, q) - lagrange_eval(&xs, &ys, q).unwrap()).abs() < 1e-9
94            );
95        }
96    }
97
98    #[test]
99    fn fails_closed() {
100        assert_eq!(
101            lagrange_eval(&[], &[], 0.0).unwrap_err(),
102            InterpolationError::InsufficientData
103        );
104        assert_eq!(
105            lagrange_eval(&[1.0, 1.0], &[2.0, 3.0], 0.0).unwrap_err(),
106            InterpolationError::DuplicateNodes
107        );
108    }
109}
qualia_core_db/solvers/interpolation/lagrange.rs

qualia_core_db/solvers/interpolation/
lagrange.rs
