polynominal_interpolation/
lagrange.rs1use 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}