diff --git a/Cargo.toml b/Cargo.toml index d62cd4affa7..1d6f2b3584d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -10,6 +10,7 @@ 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" @@ -17,4 +18,4 @@ quickcheck_macros = "1.0" [features] default = ["big-math"] -big-math = ["dep:num-bigint", "dep:num-traits"] \ No newline at end of file +big-math = ["dep:num-bigint", "dep:num-traits"] diff --git a/src/math/least_square_approx.rs b/src/math/least_square_approx.rs new file mode 100644 index 00000000000..424bd8c3f07 --- /dev/null +++ b/src/math/least_square_approx.rs @@ -0,0 +1,101 @@ +/// Least Square Approximation

+/// 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 { + 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::::new(); + for i in 1..=(2 * degree + 1) { + sums.push(points.iter().map(|(x, _)| x.powi(i - 1)).sum()); + } + + let mut free_col = Vec::::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 = 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]); + } +} diff --git a/src/math/logarithm.rs b/src/math/logarithm.rs new file mode 100644 index 00000000000..bc2ec73227f --- /dev/null +++ b/src/math/logarithm.rs @@ -0,0 +1,78 @@ +use std::f64::consts::E; + +/// Calculates the **logbase(x)** +/// +/// Parameters: +///

-> base: base of log +///

-> x: value for which log shall be evaluated +///

-> 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); + } +} diff --git a/src/math/mod.rs b/src/math/mod.rs index 5d292c437aa..f27d1a26ca9 100644 --- a/src/math/mod.rs +++ b/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;