Add linear sieve (rust-lang#288)

Co-authored-by: Andrii Siriak <siryaka@gmail.com>

Author: erfan-khadem

DIRECTORY.md

@@ -61,6 +61,7 @@
     * [Prime Numbers](https://github.com/TheAlgorithms/Rust/blob/master/src/math/prime_numbers.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)
   * Searching
     * [Binary Search](https://github.com/TheAlgorithms/Rust/blob/master/src/searching/binary_search.rs)
     * [Binary Search Recursive](https://github.com/TheAlgorithms/Rust/blob/master/src/searching/binary_search_recursive.rs)

README.md

@@ -42,6 +42,7 @@ RESTART BUILD
 - [x] [Pascal's triangle](./src/math/pascal_triangle.rs)
 - [x] [Square root with Newton's method](./src/math/square_root.rs)
 - [x] [Fast power algorithm](./src/math/fast_power.rs)
+- [x] [Linear Sieve](./src/math/linear_sieve.rs)
 
 ## [Dynamic Programming](./src/dynamic_programming)
 

src/math/linear_sieve.rs

@@ -0,0 +1,128 @@
+/*
+Linear Sieve algorithm:
+Time complexity is indeed O(n) with O(n) memory, but the sieve generally
+runs slower than a well implemented sieve of Eratosthenes. Some use cases are:
+- factorizing any number k in the sieve in O(log(k))
+- calculating arbitrary multiplicative functions on sieve numbers
+  without increasing the time complexity
+- As a by product, all prime numbers less than `max_number` are stored
+  in `primes` vector.
+ */
+pub struct LinearSieve {
+    max_number: usize,
+    primes: Vec<usize>,
+    minimum_prime_factor: Vec<usize>,
+}
+
+impl LinearSieve {
+    pub const fn new() -> Self {
+        LinearSieve {
+            max_number: 0,
+            primes: vec![],
+            minimum_prime_factor: vec![],
+        }
+    }
+
+    pub fn prepare(&mut self, max_number: usize) -> Result<(), &'static str> {
+        if max_number <= 1 {
+            return Err("Sieve size should be more than 1");
+        }
+        if self.max_number > 0 {
+            return Err("Sieve already initialized");
+        }
+        self.max_number = max_number;
+        self.minimum_prime_factor.resize(max_number + 1, 0);
+        for i in 2..=max_number {
+            if self.minimum_prime_factor[i] == 0 {
+                self.minimum_prime_factor[i] = i;
+                self.primes.push(i);
+                /*
+                   if needed, a multiplicative function can be
+                   calculated for this prime number here:
+                   function[i] = base_case(i);
+                */
+            }
+            for p in self.primes.iter() {
+                let mlt = (*p) * i;
+                if *p > i || mlt > max_number {
+                    break;
+                }
+                self.minimum_prime_factor[mlt] = *p;
+                /*
+                   multiplicative function for mlt can be calculated here:
+                   if i % p:
+                       function[mlt] = add_to_prime_exponent(function[i], i, p);
+                   else:
+                       function[mlt] = function[i] * function[p]
+                */
+            }
+        }
+        Ok(())
+    }
+
+    pub fn factorize(&self, mut number: usize) -> Result<Vec<usize>, &'static str> {
+        if number > self.max_number {
+            return Err("Number is too big, its minimum_prime_factor was not calculated");
+        }
+        if number == 0 {
+            return Err("Number is zero");
+        }
+        let mut result: Vec<usize> = Vec::new();
+        while number > 1 {
+            result.push(self.minimum_prime_factor[number]);
+            number /= self.minimum_prime_factor[number];
+        }
+        Ok(result)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::LinearSieve;
+
+    #[test]
+    fn small_primes_list() {
+        let mut ls = LinearSieve::new();
+        ls.prepare(25).unwrap();
+        assert_eq!(ls.primes, vec![2, 3, 5, 7, 11, 13, 17, 19, 23]);
+    }
+
+    #[test]
+    fn divisible_by_mpf() {
+        let mut ls = LinearSieve::new();
+        ls.prepare(1000).unwrap();
+        for i in 2..=1000 {
+            let div = i / ls.minimum_prime_factor[i];
+            assert_eq!(i % ls.minimum_prime_factor[i], 0);
+            if div == 1 {
+                // Number must be prime
+                assert!(ls.primes.binary_search(&i).is_ok());
+            }
+        }
+    }
+
+    #[test]
+    fn check_factorization() {
+        let mut ls = LinearSieve::new();
+        ls.prepare(1000).unwrap();
+        for i in 1..=1000 {
+            let factorization = ls.factorize(i).unwrap();
+            let mut product = 1usize;
+            for (idx, p) in factorization.iter().enumerate() {
+                assert!(ls.primes.binary_search(&p).is_ok());
+                product *= *p;
+                if idx > 0 {
+                    assert!(*p >= factorization[idx - 1]);
+                }
+            }
+            assert_eq!(product, i);
+        }
+    }
+
+    #[test]
+    fn check_number_of_primes() {
+        let mut ls = LinearSieve::new();
+        ls.prepare(100_000).unwrap();
+        assert_eq!(ls.primes.len(), 9592);
+    }
+}

src/math/mod.rs

@@ -1,6 +1,7 @@
 mod extended_euclidean_algorithm;
 mod fast_power;
 mod greatest_common_divisor;
+mod linear_sieve;
 mod miller_rabin;
 mod pascal_triangle;
 mod perfect_numbers;
@@ -14,6 +15,7 @@ pub use self::fast_power::fast_power;
 pub use self::greatest_common_divisor::{
     greatest_common_divisor_iterative, greatest_common_divisor_recursive,
 };
+pub use self::linear_sieve::LinearSieve;
 pub use self::miller_rabin::miller_rabin;
 pub use self::pascal_triangle::pascal_triangle;
 pub use self::perfect_numbers::perfect_numbers;