polynominal_interpolation/
lagrange.rs

1use itertools::multizip;
2use rayon::iter::{
3    IndexedParallelIterator, IntoParallelIterator, IntoParallelRefIterator, ParallelBridge,
4    ParallelIterator,
5};
6
7fn get_lagrange_barycentric_weights(xs: &[f64]) -> Vec<f64> {
8    xs.par_iter()
9        .enumerate()
10        .map(|(index, &x_index)| {
11            let inversed_weight: f64 = xs
12                .par_iter()
13                .enumerate()
14                .filter(|&(i, _)| i != index)
15                .map(|(_, &x_i)| x_i)
16                .map(|x_i| (x_index - x_i))
17                .product();
18
19            1. / inversed_weight
20        })
21        .collect()
22}
23
24pub fn lagrange_interpolation(xs: Vec<f64>, ys: Vec<f64>) -> impl Fn(f64) -> f64 {
25    let weights = get_lagrange_barycentric_weights(&xs);
26
27    move |x| {
28        if let Some(i) = xs.par_iter().position_first(|&x_i| x == x_i) {
29            return ys[i];
30        }
31
32        let numerator: f64 = multizip((&weights, &xs, &ys))
33            .par_bridge()
34            .map(|(&w_i, &x_i, &y_i)| (w_i * y_i) / (x - x_i))
35            .sum();
36        let denominator: f64 = (&weights, &xs)
37            .into_par_iter()
38            .map(|(&w_i, &x_i)| w_i / (x - x_i))
39            .sum();
40
41        numerator / denominator
42    }
43}