Machine prime integration (#86)

* Added Machine-Prime, simplified nearest_prime

Added the low-memory variant of Machine-prime (using a modified BPSW test) as an optional dependency. The Criterion benchmark shows that it is approximately 16x faster than the existing version. 

Simplified the next and previous_prime functions

* cargo fmt

* Rename feature to `fast_test`

* Add the `fast_test` feature to docs

* Remove commented code

* Correct placement of whitespace

* Use a typed stride to ensure that `stride` is only ever 1 or `u64::MAX`

* Sentences start with capital letters

* Add changes to log

---------

Co-authored-by: JASory <[email protected]>

Author: JSorngard

CHANGELOG.md

@@ -2,6 +2,11 @@
 
 This file contains the changes to the crate since version 0.4.8.
 
+## 0.9.7
+
+- Added the `fast_test` feature that makes `is_prime` call out to the [`machine-prime`](https://crates.io/crates/machine_prime)
+ crate for a significant speedup.
+
 ## 0.9.6
 
 - Correct function name in README.

Cargo.lock

@@ -1,7 +1,7 @@
 [package]
 name = "const-primes"
 authors = ["Johanna Sörngård <[email protected]>"]
-version = "0.9.6"
+version = "0.9.7"
 edition = "2021"
 license = "MIT OR Apache-2.0"
 keywords = ["const", "primes", "no_std", "prime-numbers"]
@@ -16,6 +16,8 @@ serde = { version = "1.0", default-features = false, features = ["derive"], opti
 serde_arrays = { version = "0.1.0", optional = true }
 zerocopy = { version = "0.8", default-features = false, features = ["derive"], optional = true }
 rkyv = { version = "0.8", default-features = false, optional = true }
+# no-std constant time primality testing, low memory variant
+machine-prime = { version="1.3.*", features = ["small"], optional = true}
 
 [dev-dependencies]
 criterion = { version = "0.5", features = ["html_reports"] }
@@ -33,6 +35,10 @@ zerocopy = ["dep:zerocopy"]
 # Derives the `Serialize`, `Deserialize`, and `Archive` traits from the [`rkyv`](https://crates.io/crates/rkyv) crate for the `Primes` struct.
 rkyv = ["dep:rkyv"]
 
+# Speeds up primality testing significantly by using the [`machine-prime`](https://crates.io/crates/machine-prime) crate.
+fast_test = ["dep:machine-prime"]
+
+
 [package.metadata.docs.rs]
 # Document all features.
 all-features = true

Cargo.toml

@@ -114,6 +114,9 @@ for the `Primes` struct.
 `rkyv`: derives the `Serialize`, `Deserialize`, and `Archive` traits from
 [`rkyv`](https://crates.io/crates/rkyv) for the `Primes` struct.
 
+`fast_test`: speeds up primality testing significantly by using the
+[`machine-prime`](https://crates.io/crates/machine-prime) crate.
+
 <br>
 
 ### License

README.md

@@ -1,5 +1,5 @@
 //! This module contains an implementation of a deterministic Miller-Rabin primality test
-
+#[cfg(not(feature = "fast_test"))]
 use crate::integer_math::{mod_mul, mod_pow};
 
 /// Returns whether `n` is prime.
@@ -17,67 +17,76 @@ use crate::integer_math::{mod_mul, mod_pow};
 /// ```
 #[must_use]
 pub const fn is_prime(n: u64) -> bool {
-    // Since we know the maximum size of the numbers we test against
-    // we can use the fact that there are known perfect bases
-    // in order to make the test both fast and deterministic.
-    // This list of witnesses was taken from
-    // <https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases>.
-    const NUM_BASES: usize = 11;
-    const WITNESSES: [(u64, &[u64]); NUM_BASES] = [
-        (2_046, &[2]),
-        (1_373_652, &[2, 3]),
-        (9_080_190, &[31, 73]),
-        (25_326_000, &[2, 3, 5]),
-        (4_759_123_140, &[2, 7, 61]),
-        (1_112_004_669_632, &[2, 13, 23, 1_662_803]),
-        (2_152_302_898_746, &[2, 3, 5, 7, 11]),
-        (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]),
-        (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]),
-        (3_825_123_056_546_413_050, &[2, 3, 5, 7, 11, 13, 17, 19, 23]),
-        (u64::MAX, &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]),
-    ];
-
-    if n == 2 || n == 3 {
-        return true;
-    } else if n <= 1 || n % 2 == 0 || n % 3 == 0 {
-        return false;
+    #[cfg(feature = "fast_test")]
+    {
+        machine_prime::is_prime(n)
     }
 
-    // Use a small wheel to check up to log2(n).
-    // This keeps the complexity at O(log(n)).
-    let mut candidate_factor = 5;
-    let trial_limit = n.ilog2() as u64;
-    while candidate_factor <= trial_limit {
-        if n % candidate_factor == 0 || n % (candidate_factor + 2) == 0 {
+    #[cfg(not(feature = "fast_test"))]
+    {
+        // Since we know the maximum size of the numbers we test against
+        // we can use the fact that there are known perfect bases
+        // in order to make the test both fast and deterministic.
+        // This list of witnesses was taken from
+        // <https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases>.
+        const NUM_BASES: usize = 11;
+        const WITNESSES: [(u64, &[u64]); NUM_BASES] = [
+            (2_046, &[2]),
+            (1_373_652, &[2, 3]),
+            (9_080_190, &[31, 73]),
+            (25_326_000, &[2, 3, 5]),
+            (4_759_123_140, &[2, 7, 61]),
+            (1_112_004_669_632, &[2, 13, 23, 1_662_803]),
+            (2_152_302_898_746, &[2, 3, 5, 7, 11]),
+            (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]),
+            (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]),
+            (3_825_123_056_546_413_050, &[2, 3, 5, 7, 11, 13, 17, 19, 23]),
+            (u64::MAX, &[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]),
+        ];
+
+        if n == 2 || n == 3 {
+            return true;
+        } else if n <= 1 || n % 2 == 0 || n % 3 == 0 {
             return false;
         }
-        candidate_factor += 6;
-    }
 
-    // Find r such that n = 2^d * r + 1 for some r >= 1
-    let mut d = n - 1;
-    while d % 2 == 0 {
-        d >>= 1;
-    }
+        // Use a small wheel to check up to log2(n).
+        // This keeps the complexity at O(log(n)).
+        let mut candidate_factor = 5;
+        let trial_limit = n.ilog2() as u64;
+        while candidate_factor <= trial_limit {
+            if n % candidate_factor == 0 || n % (candidate_factor + 2) == 0 {
+                return false;
+            }
+            candidate_factor += 6;
+        }
 
-    let mut i = 0;
-    while i < NUM_BASES && WITNESSES[i].0 < n {
-        i += 1;
-    }
-    let witnesses = WITNESSES[i].1;
+        // Find r such that n = 2^d * r + 1 for some r >= 1
+        let mut d = n - 1;
+        while d % 2 == 0 {
+            d >>= 1;
+        }
 
-    let mut i = 0;
-    while i < witnesses.len() && witnesses[i] < n {
-        if !miller_test(d, n, witnesses[i]) {
-            return false;
+        let mut i = 0;
+        while i < NUM_BASES && WITNESSES[i].0 < n {
+            i += 1;
+        }
+        let witnesses = WITNESSES[i].1;
+
+        let mut i = 0;
+        while i < witnesses.len() && witnesses[i] < n {
+            if !miller_test(d, n, witnesses[i]) {
+                return false;
+            }
+            i += 1;
         }
-        i += 1;
-    }
 
-    true
+        true
+    } // end conditional compilation block
 }
 
 /// Performs a Miller-Rabin test with the witness k.
+#[cfg(not(feature = "fast_test"))]
 const fn miller_test(mut d: u64, n: u64, k: u64) -> bool {
     let mut x = mod_pow(k, d, n);
     if x == 1 || x == n - 1 {

src/check.rs

@@ -34,6 +34,7 @@ pub const fn isqrt(n: u64) -> u64 {
 
 /// Calculates (`base` ^ `exp`) mod `modulo` without overflow.
 #[must_use]
+#[cfg(not(feature = "fast_test"))]
 pub const fn mod_pow(mut base: u64, mut exp: u64, modulo: u64) -> u64 {
     let mut res = 1;
 
@@ -52,6 +53,7 @@ pub const fn mod_pow(mut base: u64, mut exp: u64, modulo: u64) -> u64 {
 
 /// Calculates (`a` * `b`) mod `modulo` without overflow.
 #[must_use]
+#[cfg(not(feature = "fast_test"))]
 pub const fn mod_mul(a: u64, b: u64, modulo: u64) -> u64 {
     ((a as u128 * b as u128) % modulo as u128) as u64
 }

src/integer_math.rs

@@ -97,6 +97,8 @@
 //! `zerocopy`: derives the [`IntoBytes`](zerocopy::IntoBytes) trait from the [`zerocopy`] crate for the [`Primes`] struct.
 //!
 //! `rkyv`: derives the [`Serialize`](rkyv::Serialize), [`Deserialize`](rkyv::Deserialize), and [`Archive`](rkyv::Archive) traits from the [`rkyv`] crate for the [`Primes`] struct.
+//!
+//! `fast_test`: speeds up primality testing significantly by using the [`machine-prime`](https://crates.io/crates/machine-prime) crate.
 
 #![forbid(unsafe_code)]
 #![no_std]

src/lib.rs

@@ -2,6 +2,27 @@
 
 use crate::is_prime;
 
+/// Generalised function for nearest search by incrementing/decrementing by 1
+/// Any attempt at optimising this would be largely pointless since the largest prime gap under 2^64 is only 1550
+/// And is_prime's trial division already eliminates most of those
+const fn bounded_search(mut n: u64, stride: Stride) -> Option<u64> {
+    let stride = stride.into_u64();
+
+    loop {
+        // Addition over Z/2^64, aka regular addition under optimisation flags
+        n = n.wrapping_add(stride);
+
+        // If either condition is met then we started either below or above the smallest or largest prime respectively
+        // Any two values from 2^64-58 to 1 would also work
+        if n == 0u64 || n == u64::MAX {
+            return None;
+        }
+
+        if is_prime(n) {
+            return Some(n);
+        }
+    }
+}
 /// Returns the largest prime smaller than `n` if there is one.
 ///
 /// Scans for primes downwards from the input with [`is_prime`].
@@ -17,31 +38,15 @@ use crate::is_prime;
 /// ```
 ///
 /// There's no prime smaller than two:
-/// ```
 ///
+/// ```
 /// # use const_primes::previous_prime;
 /// const NO_SUCH: Option<u64> = previous_prime(2);
 /// assert_eq!(NO_SUCH, None);
 /// ```
 #[must_use = "the function only returns a new value and does not modify its input"]
-pub const fn previous_prime(mut n: u64) -> Option<u64> {
-    if n <= 2 {
-        None
-    } else if n == 3 {
-        Some(2)
-    } else {
-        n -= 1;
-
-        if n % 2 == 0 {
-            n -= 1;
-        }
-
-        while !is_prime(n) {
-            n -= 2;
-        }
-
-        Some(n)
-    }
+pub const fn previous_prime(n: u64) -> Option<u64> {
+    bounded_search(n, Stride::Down)
 }
 
 /// Returns the smallest prime greater than `n` if there is one that
@@ -60,30 +65,28 @@ pub const fn previous_prime(mut n: u64) -> Option<u64> {
 /// ```
 ///
 /// Primes larger than 18446744073709551557 can not be represented by a `u64`:
-/// ```
 ///
+/// ```
 /// # use const_primes::next_prime;
 /// const NO_SUCH: Option<u64> = next_prime(18_446_744_073_709_551_557);
 /// assert_eq!(NO_SUCH, None);
 /// ```
 #[must_use = "the function only returns a new value and does not modify its input"]
-pub const fn next_prime(mut n: u64) -> Option<u64> {
-    // The largest prime smaller than u64::MAX
-    if n >= 18_446_744_073_709_551_557 {
-        None
-    } else if n <= 1 {
-        Some(2)
-    } else {
-        n += 1;
+pub const fn next_prime(n: u64) -> Option<u64> {
+    bounded_search(n, Stride::Up)
+}
 
-        if n % 2 == 0 {
-            n += 1;
-        }
+enum Stride {
+    Up,
+    Down,
+}
 
-        while !is_prime(n) {
-            n += 2;
+impl Stride {
+    const fn into_u64(self) -> u64 {
+        match self {
+            Self::Up => 1,
+            // Adding by 2^64-1 over Z/2^64 is equivalent to subtracting by 1
+            Self::Down => u64::MAX,
         }
-
-        Some(n)
     }
 }