numra_interp/
polynomial.rs

1//! Polynomial interpolation via barycentric Lagrange formula.
2//!
3//! Author: Moussa Leblouba
4//! Date: 9 February 2026
5//! Modified: 2 May 2026
6
7use numra_core::Scalar;
8
9use crate::error::InterpError;
10use crate::{validate_data, Interpolant};
11
12/// Barycentric Lagrange interpolant.
13///
14/// Evaluates the unique polynomial of degree `n-1` passing through `n` data points
15/// using the numerically stable barycentric formula. O(n) per evaluation after O(n^2) setup.
16pub struct BarycentricLagrange<S: Scalar> {
17    x: Vec<S>,
18    y: Vec<S>,
19    w: Vec<S>, // barycentric weights
20}
21
22impl<S: Scalar> BarycentricLagrange<S> {
23    /// Create a barycentric Lagrange interpolant from data points.
24    ///
25    /// # Errors
26    ///
27    /// Returns error if fewer than 1 point or data is not sorted.
28    pub fn new(x: &[S], y: &[S]) -> Result<Self, InterpError> {
29        validate_data(x, y, 1)?;
30        let w = compute_barycentric_weights(x);
31        Ok(Self {
32            x: x.to_vec(),
33            y: y.to_vec(),
34            w,
35        })
36    }
37
38    /// Create from Chebyshev nodes on `[a, b]`.
39    ///
40    /// Uses `n` Chebyshev nodes of the first kind, which minimize the
41    /// maximum interpolation error (near-optimal for polynomial interpolation).
42    /// Uses the known closed-form barycentric weights for numerical stability.
43    pub fn chebyshev<F: Fn(S) -> S>(f: F, a: S, b: S, n: usize) -> Self {
44        let mut x = Vec::with_capacity(n);
45        let mut y = Vec::with_capacity(n);
46        let mut orig_weights = Vec::with_capacity(n);
47
48        for k in 0..n {
49            // Chebyshev nodes: x_k = cos((2k+1)/(2n) * pi), mapped to [a,b]
50            let theta_f = (2 * k + 1) as f64 / (2 * n) as f64 * core::f64::consts::PI;
51            let theta = S::from_f64(theta_f);
52            let node = (a + b) * S::HALF + (b - a) * S::HALF * theta.cos();
53            x.push(node);
54            y.push(f(node));
55            // Closed-form barycentric weight: (-1)^k * sin((2k+1)*pi/(2n))
56            let sign = if k % 2 == 0 { 1.0 } else { -1.0 };
57            orig_weights.push(sign * theta_f.sin());
58        }
59
60        // Sort by x (Chebyshev nodes come in decreasing order from cos)
61        let mut indices: Vec<usize> = (0..n).collect();
62        indices.sort_by(|&i, &j| x[i].to_f64().partial_cmp(&x[j].to_f64()).unwrap());
63        let x_sorted: Vec<S> = indices.iter().map(|&i| x[i]).collect();
64        let y_sorted: Vec<S> = indices.iter().map(|&i| y[i]).collect();
65        let w: Vec<S> = indices
66            .iter()
67            .map(|&i| S::from_f64(orig_weights[i]))
68            .collect();
69
70        Self {
71            x: x_sorted,
72            y: y_sorted,
73            w,
74        }
75    }
76}
77
78impl<S: Scalar> Interpolant<S> for BarycentricLagrange<S> {
79    fn interpolate(&self, x: S) -> S {
80        let n = self.x.len();
81
82        // Check if x is exactly at a node
83        for i in 0..n {
84            let diff = x - self.x[i];
85            if diff.abs() < S::EPSILON * S::from_f64(100.0) {
86                return self.y[i];
87            }
88        }
89
90        // Barycentric formula: p(x) = (sum w_i * y_i / (x - x_i)) / (sum w_i / (x - x_i))
91        let mut numer = S::ZERO;
92        let mut denom = S::ZERO;
93        for i in 0..n {
94            let t = self.w[i] / (x - self.x[i]);
95            numer += t * self.y[i];
96            denom += t;
97        }
98        numer / denom
99    }
100}
101
102/// Compute barycentric weights.
103/// Computes in f64 and normalizes to avoid cancellation in the barycentric formula.
104fn compute_barycentric_weights<S: Scalar>(x: &[S]) -> Vec<S> {
105    let n = x.len();
106    let xf: Vec<f64> = x.iter().map(|xi| xi.to_f64()).collect();
107    let mut wf = vec![1.0_f64; n];
108    for i in 0..n {
109        for j in 0..n {
110            if i != j {
111                wf[i] /= xf[i] - xf[j];
112            }
113        }
114    }
115    // Normalize to O(1) to reduce cancellation in the barycentric summation
116    let max_abs = wf.iter().map(|w| w.abs()).fold(0.0_f64, f64::max);
117    if max_abs > 0.0 {
118        for w in &mut wf {
119            *w /= max_abs;
120        }
121    }
122    wf.iter().map(|&w| S::from_f64(w)).collect()
123}
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128    use approx::assert_relative_eq;
129
130    #[test]
131    fn test_lagrange_at_knots() {
132        let x = vec![0.0, 1.0, 2.0, 3.0];
133        let y = vec![1.0, 3.0, 2.0, 4.0];
134        let p = BarycentricLagrange::new(&x, &y).unwrap();
135        for (&xi, &yi) in x.iter().zip(y.iter()) {
136            assert_relative_eq!(p.interpolate(xi), yi, epsilon = 1e-10);
137        }
138    }
139
140    #[test]
141    fn test_lagrange_polynomial_exact() {
142        // Polynomial of degree 2: y = x^2
143        let x = vec![0.0, 1.0, 2.0];
144        let y = vec![0.0, 1.0, 4.0];
145        let p = BarycentricLagrange::new(&x, &y).unwrap();
146        // Should reproduce x^2 exactly
147        assert_relative_eq!(p.interpolate(0.5), 0.25, epsilon = 1e-10);
148        assert_relative_eq!(p.interpolate(1.5), 2.25, epsilon = 1e-10);
149    }
150
151    #[test]
152    fn test_lagrange_cubic_exact() {
153        // Polynomial of degree 3: y = x^3 - x
154        let x = vec![0.0, 1.0, 2.0, 3.0];
155        let y: Vec<f64> = x.iter().map(|&xi| xi.powi(3) - xi).collect();
156        let p = BarycentricLagrange::new(&x, &y).unwrap();
157        assert_relative_eq!(p.interpolate(0.5), 0.5_f64.powi(3) - 0.5, epsilon = 1e-10);
158        assert_relative_eq!(p.interpolate(2.5), 2.5_f64.powi(3) - 2.5, epsilon = 1e-10);
159    }
160
161    #[test]
162    fn test_chebyshev_sin() {
163        // Interpolate sin(x) on [0, pi] with 10 Chebyshev nodes
164        let p = BarycentricLagrange::chebyshev(|x: f64| x.sin(), 0.0, core::f64::consts::PI, 10);
165        // Check at several points
166        for k in 1..10 {
167            let x = k as f64 * core::f64::consts::PI / 10.0;
168            let err = (p.interpolate(x) - x.sin()).abs();
169            assert!(err < 1e-7, "Chebyshev error at x={}: {}", x, err);
170        }
171    }
172
173    #[test]
174    fn test_chebyshev_runge() {
175        // Runge's function 1/(1+25x^2) on [-1, 1]
176        // Chebyshev nodes should handle this much better than equidistant.
177        // Need ~50 nodes for high accuracy due to poles at x = ±i/5 in C.
178        let p = BarycentricLagrange::chebyshev(|x: f64| 1.0 / (1.0 + 25.0 * x * x), -1.0, 1.0, 50);
179        // Check at x=0 (should be 1.0)
180        assert_relative_eq!(p.interpolate(0.0), 1.0, epsilon = 1e-3);
181        // Check near the edges (where equidistant nodes fail badly)
182        let y_edge = p.interpolate(0.9);
183        let exact = 1.0 / (1.0 + 25.0 * 0.81);
184        assert!(
185            (y_edge - exact).abs() < 1e-4,
186            "Runge at 0.9: {} vs {}",
187            y_edge,
188            exact
189        );
190    }
191
192    #[test]
193    fn test_lagrange_single_point() {
194        let p = BarycentricLagrange::new(&[1.0], &[5.0]).unwrap();
195        assert_relative_eq!(p.interpolate(1.0), 5.0, epsilon = 1e-14);
196        // Constant function
197        assert_relative_eq!(p.interpolate(2.0), 5.0, epsilon = 1e-14);
198    }
199
200    #[test]
201    fn test_lagrange_f32() {
202        let p = BarycentricLagrange::new(&[0.0f32, 1.0, 2.0], &[0.0, 1.0, 4.0]).unwrap();
203        assert!((p.interpolate(0.5f32) - 0.25).abs() < 1e-4);
204    }
205}
numra_interp/polynomial.rs

numra_interp/
polynomial.rs
