Least Square Approximation and Logarithm evaluation

Cargo.toml

@@ -10,11 +10,12 @@ num-bigint = { version = "0.4", optional = true }
 num-traits = { version = "0.2", optional = true }
 rand = "0.8"
 rand_chacha = "0.3"
+nalgebra = "0.32.3"
 
 [dev-dependencies]
 quickcheck = "1.0"
 quickcheck_macros = "1.0"
 
 [features]
 default = ["big-math"]
-big-math = ["dep:num-bigint", "dep:num-traits"]
+big-math = ["dep:num-bigint", "dep:num-traits"]

src/math/least_square_approx.rs

@@ -0,0 +1,101 @@
+/// Least Square Approximation <p>
+/// Function that returns a polynomial which very closely passes through the given points (in 2D)
+///
+/// The result is made of coeficients, in descending order (from x^degree to free term)
+///
+/// Parameters:
+///
+/// points -> coordinates of given points
+///
+/// degree -> degree of the polynomial
+///
+pub fn least_square_approx(points: &[(f64, f64)], degree: i32) -> Vec<f64> {
+    use nalgebra::{DMatrix, DVector};
+
+    /* Used for rounding floating numbers */
+    fn round_to_decimals(value: f64, decimals: u32) -> f64 {
+        let multiplier = 10f64.powi(decimals as i32);
+        (value * multiplier).round() / multiplier
+    }
+
+    /* Computes the sums in the system of equations */
+    let mut sums = Vec::<f64>::new();
+    for i in 1..=(2 * degree + 1) {
+        sums.push(points.iter().map(|(x, _)| x.powi(i - 1)).sum());
+    }
+
+    let mut free_col = Vec::<f64>::new();
+    /* Compute the free terms column vector */
+    for i in 1..=(degree + 1) {
+        free_col.push(points.iter().map(|(x, y)| y * (x.powi(i - 1))).sum());
+    }
+    let b = DVector::from_row_slice(&free_col);
+
+    let size = (degree + 1) as usize;
+    /* Create and fill the system's matrix */
+    let a = DMatrix::from_fn(size, size, |i, j| sums[i + j]);
+
+    /* Solve the system of equations: A * x = b */
+    match a.qr().solve(&b) {
+        Some(x) => {
+            let mut rez: Vec<f64> = x.iter().map(|x| round_to_decimals(*x, 5)).collect();
+            rez.reverse();
+            rez
+        }
+        None => Vec::new(),
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn ten_points_1st_degree() {
+        let points = vec![
+            (5.3, 7.8),
+            (4.9, 8.1),
+            (6.1, 6.9),
+            (4.7, 8.3),
+            (6.5, 7.7),
+            (5.6, 7.0),
+            (5.8, 8.2),
+            (4.5, 8.0),
+            (6.3, 7.2),
+            (5.1, 8.4),
+        ];
+
+        assert_eq!(least_square_approx(&points, 1), [-0.49069, 10.44898]);
+    }
+
+    #[test]
+    fn eight_points_5th_degree() {
+        let points = vec![
+            (4f64, 8f64),
+            (8f64, 2f64),
+            (1f64, 7f64),
+            (10f64, 3f64),
+            (11.0, 0.0),
+            (7.0, 3.0),
+            (10.0, 1.0),
+            (13.0, 13.0),
+        ];
+
+        assert_eq!(
+            least_square_approx(&points, 5),
+            [0.00603, -0.21304, 2.79929, -16.53468, 40.29473, -19.35771]
+        );
+    }
+
+    #[test]
+    fn four_points_2nd_degree() {
+        let points = vec![
+            (2.312, 8.345344),
+            (-2.312, 8.345344),
+            (-0.7051, 3.49716601),
+            (0.7051, 3.49716601),
+        ];
+
+        assert_eq!(least_square_approx(&points, 2), [1.0, 0.0, 3.0]);
+    }
+}

src/math/logarithm.rs

@@ -0,0 +1,78 @@
+use std::f64::consts::E;
+
+/// Calculates the **log<sub>base</sub>(x)**
+///
+/// Parameters:
+///   <p>-> base: base of log
+///   <p>-> x: value for which log shall be evaluated
+///   <p>-> tol: tolerance; the precision of the approximation
+///
+/// Advisable to use **std::f64::consts::*** for specific bases (like 'e')
+pub fn log(base: f64, mut x: f64, tol: f64) -> f64 {
+    let mut rez: f64 = 0f64;
+
+    if x <= 0f64 || base <= 0f64 {
+        println!("Log does not support negative argument or negative base.");
+        f64::NAN
+    } else if x < 1f64 && base == E {
+        /*
+            For x in (0, 1) and base 'e', the function is using MacLaurin Series:
+            ln(|1 + x|) = Σ "(-1)^n-1 * x^n / n", for n = 1..inf
+            Substituting x with x-1 yields:
+            ln(|x|) = Σ "(-1)^n-1 * (x-1)^n / n"
+        */
+        x -= 1f64;
+
+        let mut prev_rez = 1f64;
+        let mut step: i32 = 1;
+
+        while (prev_rez - rez).abs() > tol {
+            prev_rez = rez;
+            rez += (-1f64).powi(step - 1) * x.powi(step) / step as f64;
+            step += 1;
+        }
+
+        rez
+    } else {
+        let ln_x = x.ln();
+        let ln_base = base.ln();
+
+        ln_x / ln_base
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    #[test]
+    fn basic() {
+        assert_eq!(log(E, E, 0.0), 1.0);
+        assert_eq!(log(E, E.powi(100), 0.0), 100.0);
+        assert_eq!(log(10.0, 10000.0, 0.0), 4.0);
+        assert_eq!(log(234501.0, 1.0, 1.0), 0.0);
+    }
+
+    #[test]
+    fn test_log_positive_base() {
+        assert_eq!(log(10.0, 100.0, 0.00001), 2.0);
+        assert_eq!(log(2.0, 8.0, 0.00001), 3.0);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_log_zero_base() {
+        assert_eq!(log(0.0, 100.0, 0.00001), f64::NAN);
+    }
+
+    #[test]
+    #[should_panic] // Should panic because can't compare NAN to NAN
+    fn test_log_negative_base() {
+        assert_eq!(log(-1.0, 100.0, 0.00001), f64::NAN);
+    }
+
+    #[test]
+    fn test_log_tolerance() {
+        assert_eq!(log(10.0, 100.0, 1e-10), 2.0);
+    }
+}

src/math/mod.rs

@@ -37,7 +37,9 @@ mod interquartile_range;
 mod karatsuba_multiplication;
 mod lcm_of_n_numbers;
 mod leaky_relu;
+mod least_square_approx;
 mod linear_sieve;
+mod logarithm;
 mod lucas_series;
 mod matrix_ops;
 mod mersenne_primes;
@@ -118,7 +120,9 @@ pub use self::interquartile_range::interquartile_range;
 pub use self::karatsuba_multiplication::multiply;
 pub use self::lcm_of_n_numbers::lcm;
 pub use self::leaky_relu::leaky_relu;
+pub use self::least_square_approx::least_square_approx;
 pub use self::linear_sieve::LinearSieve;
+pub use self::logarithm::log;
 pub use self::lucas_series::dynamic_lucas_number;
 pub use self::lucas_series::recursive_lucas_number;
 pub use self::matrix_ops::Matrix;