diff --git a/DIRECTORY.md b/DIRECTORY.md index 746db33f9b8..bd3c98b36f8 100644 --- a/DIRECTORY.md +++ b/DIRECTORY.md @@ -67,13 +67,16 @@ * [Lib](https://github.com/TheAlgorithms/Rust/blob/master/src/lib.rs) * Math + * [Amicable Numbers Below N](https://github.com/TheAlgorithms/Rust/blob/master/src/math/amicable_numbers.rs) * [Baby-Step Giant-Step Algorithm](https://github.com/TheAlgorithms/Rust/blob/master/src/math/baby_step_giant_step.rs) * [Extended Euclidean Algorithm](https://github.com/TheAlgorithms/Rust/blob/master/src/math/extended_euclidean_algorithm.rs) + * [Fast Inverse Square Root 'Quake' Algorithm](https://github.com/TheAlgorithms/Rust/blob/master/src/math/square_root.rs) * [Greatest Common Divisor](https://github.com/TheAlgorithms/Rust/blob/master/src/math/greatest_common_divisor.rs) * [Pascal Triangle](https://github.com/TheAlgorithms/Rust/blob/master/src/math/pascal_triangle.rs) * [Perfect Numbers](https://github.com/TheAlgorithms/Rust/blob/master/src/math/perfect_numbers.rs) * [Prime Check](https://github.com/TheAlgorithms/Rust/blob/master/src/math/prime_check.rs) * [Prime Numbers](https://github.com/TheAlgorithms/Rust/blob/master/src/math/prime_numbers.rs) + * [Square root with Newton's method](https://github.com/TheAlgorithms/Rust/blob/master/src/math/square_root.rs) * [Trial Division](https://github.com/TheAlgorithms/Rust/blob/master/src/math/trial_division.rs) * [Miller Rabin](https://github.com/TheAlgorithms/Rust/blob/master/src/math/miller_rabin.rs) * [Linear Sieve](https://github.com/TheAlgorithms/Rust/blob/master/src/math/linear_sieve.rs) diff --git a/README.md b/README.md index 29bb483c7ca..7c0e1d6f50a 100644 --- a/README.md +++ b/README.md @@ -50,8 +50,10 @@ These are for demonstration purposes only. ## [Math](./src/math) +- [x] [Amicable numbers below N](./src/math/amicable_numbers.rs) - [x] [Baby-Step Giant-Step Algorithm](./src/math/baby_step_giant_step.rs) - [x] [Extended euclidean algorithm](./src/math/extended_euclidean_algorithm.rs) +- [x] [Fast Inverse Square Root 'Quake' Algorithm](./src/math/square_root.rs) - [x] [Gaussian Elimination](./src/math/gaussian_elimination.rs) - [x] [Greatest common divisor](./src/math/greatest_common_divisor.rs) - [x] [Greatest common divisor of n numbers](./src/math/gcd_of_n_numbers.rs) diff --git a/src/math/amicable_numbers.rs b/src/math/amicable_numbers.rs new file mode 100644 index 00000000000..4d466cd4134 --- /dev/null +++ b/src/math/amicable_numbers.rs @@ -0,0 +1,69 @@ +// Operations based around amicable numbers +// Suports u32 but should be interchangable with other types +// Wikipedia reference: https://en.wikipedia.org/wiki/Amicable_numbers + +// Returns vec of amicable pairs below N +// N must be positive +pub fn amicable_pairs_under_n(n: u32) -> Option> { + let mut factor_sums = vec![0; n as usize]; + + // Make a list of the sum of the factors of each number below N + for i in 1..n { + for j in (i * 2..n).step_by(i as usize) { + factor_sums[j as usize] += i; + } + } + + // Default value of (0, 0) if no pairs are found + let mut out = vec![(0, 0)]; + // Check if numbers are amicable then append + for (i, x) in factor_sums.iter().enumerate() { + if (*x != i as u32) && (*x < n) && (factor_sums[*x as usize] == i as u32) && (*x > i as u32) + { + out.push((i as u32, *x)); + } + } + + // Check if anything was added to the vec, if so remove the (0, 0) and return + if out.len() == 1 { + None + } else { + out.remove(0); + Some(out) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + pub fn test_amicable_numbers_below_n() { + // First 10 amicable numbers, sorted (low, high) + let expected_result = vec![ + (220, 284), + (1184, 1210), + (2620, 2924), + (5020, 5564), + (6232, 6368), + (10744, 10856), + (12285, 14595), + (17296, 18416), + (63020, 76084), + (66928, 66992), + ]; + + // Generate pairs under 100,000 + let mut result = amicable_pairs_under_n(100_000).unwrap(); + + // There should be 13 pairs under 100,000 + assert_eq!(result.len(), 13); + + // Check the first 10 against known values + result = result[..10].to_vec(); + assert_eq!(result, expected_result); + + // N that does not have any amicable pairs below it, the result should be None + assert_eq!(amicable_pairs_under_n(100), None); + } +} diff --git a/src/math/mod.rs b/src/math/mod.rs index 821034673e1..efcf841cf14 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -1,3 +1,4 @@ +mod amicable_numbers; mod armstrong_number; mod baby_step_giant_step; mod extended_euclidean_algorithm; @@ -29,6 +30,7 @@ mod square_root; mod trial_division; mod zellers_congruence_algorithm; +pub use self::amicable_numbers::amicable_pairs_under_n; pub use self::armstrong_number::is_armstrong_number; pub use self::baby_step_giant_step::baby_step_giant_step; pub use self::extended_euclidean_algorithm::extended_euclidean_algorithm; @@ -63,6 +65,6 @@ pub use self::quadratic_residue::cipolla; pub use self::random::PCG32; pub use self::sieve_of_eratosthenes::sieve_of_eratosthenes; pub use self::simpson_integration::simpson_integration; -pub use self::square_root::square_root; +pub use self::square_root::{fast_inv_sqrt, square_root}; pub use self::trial_division::trial_division; pub use self::zellers_congruence_algorithm::zellers_congruence_algorithm; diff --git a/src/math/square_root.rs b/src/math/square_root.rs index 908f8d55987..42d87e395ef 100644 --- a/src/math/square_root.rs +++ b/src/math/square_root.rs @@ -14,12 +14,43 @@ pub fn square_root(num: f64) -> f64 { root } +// fast_inv_sqrt returns an approximation of the inverse square root +// This algorithm was fist used in Quake and has been reimplimented in a few other languages +// This crate impliments it more throughly: https://docs.rs/quake-inverse-sqrt/latest/quake_inverse_sqrt/ +pub fn fast_inv_sqrt(num: f32) -> f32 { + // If you are confident in your input this can be removed for speed + if num < 0.0f32 { + return f32::NAN; + } + + let i = num.to_bits(); + let i = 0x5f3759df - (i >> 1); + let y = f32::from_bits(i); + + println!("num: {:?}, out: {:?}", num, y * (1.5 - 0.5 * num * y * y)); + // First iteration of newton approximation + y * (1.5 - 0.5 * num * y * y) + // The above can be repeated again for more precision +} + #[cfg(test)] mod tests { use super::*; #[test] - fn test() { + fn test_fast_inv_sqrt() { + // Negatives don't have square roots: + assert!(fast_inv_sqrt(-1.0f32).is_nan()); + + // Test a few cases, expect less than 1% error: + let test_pairs = [(4.0, 0.5), (16.0, 0.25), (25.0, 0.2)]; + for pair in test_pairs { + assert!((fast_inv_sqrt(pair.0) - pair.1).abs() <= (0.01 * pair.0)); + } + } + + #[test] + fn test_sqare_root() { assert!((square_root(4.0_f64) - 2.0_f64).abs() <= 1e-10_f64); assert!(square_root(-4.0_f64).is_nan()); }