polyfit_residuals/
poly.rs

1//! Basic polynomials in a Newton basis.
2
3use num_traits::One;
4use std::iter;
5use std::ops::Add;
6use std::{
7    marker::PhantomData,
8    ops::{Mul, Sub},
9};
10
11/// A Newton Polynomial that owns its basis and coeffs in form of a `Vec`.
12pub type OwnedNewtonPolynomial<C, X> = NewtonPolynomial<C, X, Vec<C>, Vec<X>>;
13
14#[cfg(feature = "nightly-only")]
15/// A Newton Polynomial of fixed Degree owning its basis and coeffs as static arrays.
16pub type StaticNewtonPolynomial<C, X, const DEGREE: usize> =
17    NewtonPolynomial<C, X, [C; DEGREE + 1], [X; DEGREE]>;
18
19/// A Newton Polynomial with borrowed basis and coeffs.
20pub type RefNewtonPolynomial<'a, C, X> = NewtonPolynomial<C, X, &'a [C], &'a [X]>;
21
22#[derive(Debug, Clone, Copy, Eq, PartialEq)]
23pub struct NewtonPolynomial<C, X, DataC, DataX>
24where
25    DataC: AsRef<[C]>,
26    DataX: AsRef<[X]>,
27{
28    coeffs: DataC,
29    basis_elems: DataX,
30    _phantom: PhantomData<(C, X)>,
31}
32
33impl<C, X, DataC, DataX> NewtonPolynomial<C, X, DataC, DataX>
34where
35    DataC: AsRef<[C]>,
36    DataX: AsRef<[X]>,
37{
38    /// Construct the polynomial
39    ///     coeffs\[0\] + coeffs\[1\] (x - basis_elems\[0\]) + coeffs\[2\] (x - basis_elems\[0\]) * (x - basis_elems\[1\]) + ...
40    /// with `coeffs.len()` degrees of freedom / of degree `basis_elems.len()`.
41    pub fn new(coeffs: DataC, basis_elems: DataX) -> Self {
42        assert_eq!(coeffs.as_ref().len(), basis_elems.as_ref().len() + 1);
43        assert!(!coeffs.as_ref().is_empty());
44        Self {
45            coeffs,
46            basis_elems,
47            _phantom: PhantomData,
48        }
49    }
50
51    /// Evaluate the polynomial "from the left" as 1ₓα₀ + (x-x₁)(α₁+(x-x₂)(α₂ + (x-x₃)(...)))
52    pub fn left_eval<Y>(&self, x: X) -> Y
53    where
54        C: Clone,
55        X: Clone + One + Sub<Output = X> + Mul<C, Output = Y> + Mul<Y, Output = Y>,
56        Y: Add<Output = Y>,
57    {
58        let mut it = self
59            .basis_elems
60            .as_ref()
61            .iter()
62            .map(|x_i| x.clone() - x_i.clone())
63            .rev()
64            .chain(iter::once(X::one()))
65            .zip(self.coeffs.as_ref().iter().rev());
66        let init = {
67            let (x, alpha) = it.next().unwrap();
68            x * alpha.clone()
69        };
70        it.fold(init, |acc, (x_i, alpha_i)| {
71            x_i * (X::one() * alpha_i.clone() + acc)
72        })
73    }
74
75    /// Evaluate the polynomial "from the right" as α₀1ₓ + (α₁ + (α₂ + (...)(x-x₃))(x-x₂))(x-x₁)
76    pub fn right_eval<Y>(&self, x: X) -> Y
77    where
78        C: Clone + Mul<X, Output = Y>,
79        X: Clone + One + Sub<Output = X> + Mul<Y, Output = Y>,
80        Y: Add<Output = Y>,
81    {
82        let mut it = self
83            .basis_elems
84            .as_ref()
85            .iter()
86            .map(|x_i| x.clone() - x_i.clone())
87            .rev()
88            .chain(iter::once(X::one()))
89            .zip(self.coeffs.as_ref().iter().rev());
90        let init = {
91            let (x, alpha) = it.next().unwrap();
92            alpha.clone() * x
93        };
94        it.fold(init, |acc, (x_i, alpha_i)| {
95            x_i * (alpha_i.clone() * X::one() + acc)
96        })
97    }
98
99    /// Evaluate the polynomial at a point
100    pub fn eval<Y>(&self, x: X) -> Y
101    where
102        C: Clone + Mul<X, Output = Y>,
103        X: Clone + One + Sub<Output = X> + Mul<C, Output = Y> + Mul<Y, Output = Y>,
104        Y: Add<Output = Y>,
105    {
106        // debug_assert_eq!(self.left_eval(x.clone()), self.right_eval(x.clone()));
107        self.left_eval(x)
108    }
109
110    /// Turn a polynomial into it's constituent coefficients and basis elements
111    pub fn into_raw(self) -> (DataC, DataX) {
112        let Self {
113            coeffs,
114            basis_elems,
115            ..
116        } = self;
117        (coeffs, basis_elems)
118    }
119}
120
121#[cfg(test)]
122mod tests {
123    use super::NewtonPolynomial;
124
125    #[test]
126    fn left_eval() {
127        let poly = NewtonPolynomial::new(vec![-1, 2, 3], vec![10, 20]);
128        assert_eq!(poly.left_eval(10), -1);
129        assert_eq!(poly.left_eval(20), 19);
130        assert_eq!(poly.left_eval(15), -66);
131        assert_eq!(poly.left_eval(2), 415);
132        assert_eq!(poly.left_eval(5), 214);
133    }
134
135    #[test]
136    fn right_eval() {
137        let poly = NewtonPolynomial::new(vec![-1, 2, 3], vec![10, 20]);
138        assert_eq!(poly.right_eval(10), -1);
139        assert_eq!(poly.right_eval(20), 19);
140        assert_eq!(poly.right_eval(15), -66);
141        assert_eq!(poly.right_eval(2), 415);
142        assert_eq!(poly.right_eval(5), 214);
143    }
144
145    #[test]
146    fn eval() {
147        let poly = NewtonPolynomial::new(vec![-1, 2, 3], vec![10, 20]);
148        assert_eq!(poly.eval(10), -1);
149        assert_eq!(poly.eval(20), 19);
150        assert_eq!(poly.eval(15), -66);
151        assert_eq!(poly.eval(2), 415);
152        assert_eq!(poly.eval(5), 214);
153    }
154}
polyfit_residuals/poly.rs

polyfit_residuals/
poly.rs
