polynominal_interpolation/
newton.rs

1use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
2
3fn get_newton_basis(x: f64, xs: &[f64]) -> f64 {
4    xs.par_iter().map(|&x_i| x - x_i).product()
5}
6
7pub fn newton_interpolation(xs: Vec<f64>, ys: Vec<f64>) -> impl Fn(f64) -> f64 {
8    let n = xs.len();
9    assert!(n != 0);
10
11    move |x| {
12        if n == 1 {
13            return ys[0];
14        }
15
16        let prev_xs = xs[..n - 1].to_vec();
17        let prev_ys = ys[..n - 1].to_vec();
18
19        let x_n = xs[n - 1];
20        let y_n = ys[n - 1];
21
22        let basis = get_newton_basis(x, &prev_xs);
23        let x_n_basis = get_newton_basis(x_n, &prev_xs);
24
25        let prev_polynominal = newton_interpolation(prev_xs, prev_ys);
26
27        let basis_coefficient = (y_n - prev_polynominal(x_n)) / x_n_basis;
28
29        prev_polynominal(x) + basis_coefficient * basis
30    }
31}