newton polynomial interpolation

day9/Cargo.toml

@@ -2,7 +2,16 @@
 name = "day9"
 version = "0.1.0"
 edition = "2021"
+default-run = "main"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
+
+[[bin]]
+name = "main"
+path= "src/main.rs"
+
+[[bin]]
+name = "newton"
+path= "src/main_newton.rs"

day9/src/main_newton.rs

@@ -0,0 +1,57 @@
+use std::fs::read_to_string;
+
+fn get_numbers_in_line_iter<T: std::str::FromStr>(str: &str) -> impl Iterator<Item=T> + '_ {
+    str.split_whitespace().filter_map(|substr| substr.parse::<T>().ok())
+}
+
+fn polynomial_value(coeff: &Vec<f64>, x: f64) -> f64 {
+    let mut result: f64 = 0.;
+    let mut x_power: f64 = 1.;
+    for &c in coeff {
+        result += c * x_power;
+        x_power *= x;
+    }
+    result
+}
+
+fn newton_interpolation(x: &Vec<f64>, y: &Vec<f64>) -> Vec<f64> {
+    let degree = std::cmp::min(x.len(), y.len());
+    if degree == 0 {
+        return vec![];
+    }
+    let mut result: Vec<f64> = Vec::with_capacity(degree);
+    let mut w: Vec<f64> = Vec::with_capacity(degree);
+    result.push(y[0]);
+    w.push(x[0]);
+    for i in 1..degree {
+        w.push(1.);
+        for j in 0..i {
+            let xj = x[j];
+            w.push(0.);
+            for k in (0..=j).rev() {
+                w[k + 1] += w[k];
+                w[k] *= -xj;
+            }
+        }
+        let a = (y[i] - polynomial_value(&result, x[i])) / polynomial_value(&w, x[i]);
+        for j in 0..i {
+            result[j] += a * w[j];
+        }
+        result.push(a * w[i]);
+        w.clear();
+    }
+    result
+}
+
+fn main() {let mut sum1 = 0;
+    let mut sum2 = 0;
+    for line in read_to_string("input.txt").unwrap().lines() {
+        let y = get_numbers_in_line_iter::<f64>(line).collect::<Vec<_>>();
+        let x = (0..y.len()).map(|x| x as f64).collect::<Vec<_>>();
+        let coeff = newton_interpolation(&x, &y);
+        sum1 += polynomial_value(&coeff, y.len() as f64).round() as i32;
+        sum2 += polynomial_value(&coeff, -1.).round() as i32;
+    }
+    println!("Sum1: {}", sum1);
+    println!("Sum2: {}", sum2);
+}