1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

//! Basic polynomials in a Newton basis.

use num_traits::One;
use std::iter;
use std::ops::Add;
use std::{
    marker::PhantomData,
    ops::{Mul, Sub},
};

/// A Newton Polynomial that owns its basis and coeffs in form of a `Vec`.
pub type OwnedNewtonPolynomial<C, X> = NewtonPolynomial<C, X, Vec<C>, Vec<X>>;

#[cfg(feature = "generic_const_exprs")]
/// A Newton Polynomial of fixed Degree owning its basis and coeffs as static arrays.
pub type StaticNewtonPolynomial<C, X, const DEGREE: usize> =
    NewtonPolynomial<C, X, [C; DEGREE + 1], [X; DEGREE]>;

/// A Newton Polynomial with borrowed basis and coeffs.
pub type RefNewtonPolynomial<'a, C, X> = NewtonPolynomial<C, X, &'a [C], &'a [X]>;

#[derive(Debug, Clone, Copy, Eq, PartialEq)]
pub struct NewtonPolynomial<C, X, DataC, DataX>
where
    DataC: AsRef<[C]>,
    DataX: AsRef<[X]>,
{
    coeffs: DataC,
    basis_elems: DataX,
    _phantom: PhantomData<(C, X)>,
}

impl<C, X, DataC, DataX> NewtonPolynomial<C, X, DataC, DataX>
where
    DataC: AsRef<[C]>,
    DataX: AsRef<[X]>,
{
    /// Construct the polynomial
    ///     coeffs\[0\] + coeffs\[1\] (x - basis_elems\[0\]) + coeffs\[2\] (x - basis_elems\[0\]) * (x - basis_elems\[1\]) + ...
    /// with `coeffs.len()` degrees of freedom / of degree `basis_elems.len()`.
    pub fn new(coeffs: DataC, basis_elems: DataX) -> Self {
        assert_eq!(coeffs.as_ref().len(), basis_elems.as_ref().len() + 1);
        assert!(!coeffs.as_ref().is_empty());
        Self {
            coeffs,
            basis_elems,
            _phantom: PhantomData,
        }
    }

    /// Evaluate the polynomial "from the left" as 1ₓα₀ + (x-x₁)(α₁+(x-x₂)(α₂ + (x-x₃)(...)))
    pub fn left_eval<Y>(&self, x: X) -> Y
    where
        C: Clone,
        X: Clone + One + Sub<Output = X> + Mul<C, Output = Y> + Mul<Y, Output = Y>,
        Y: Add<Output = Y>,
    {
        let mut it = self
            .basis_elems
            .as_ref()
            .iter()
            .map(|x_i| x.clone() - x_i.clone())
            .rev()
            .chain(iter::once(X::one()))
            .zip(self.coeffs.as_ref().iter().rev());
        let init = {
            let (x, alpha) = it.next().unwrap();
            x * alpha.clone()
        };
        it.fold(init, |acc, (x_i, alpha_i)| {
            x_i * (X::one() * alpha_i.clone() + acc)
        })
    }

    /// Turn a polynomial into it's constituent coefficients and basis elements
    pub fn into_raw(self) -> (DataC, DataX) {
        let Self {
            coeffs,
            basis_elems,
            ..
        } = self;
        (coeffs, basis_elems)
    }
}

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

    #[test]
    fn left_eval() {
        let poly = NewtonPolynomial::new(vec![-1, 2, 3], vec![10, 20]);
        assert_eq!(poly.left_eval(10), -1);
        assert_eq!(poly.left_eval(20), 19);
        assert_eq!(poly.left_eval(15), -66);
        assert_eq!(poly.left_eval(2), 415);
        assert_eq!(poly.left_eval(5), 214);
    }
}