diff --git a/CHANGELOG.md b/CHANGELOG.md index d224398..7556812 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/Cargo.lock b/Cargo.lock index 809e357..f32ca09 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -101,9 +101,10 @@ checksum = "afb84c814227b90d6895e01398aee0d8033c00e7466aca416fb6a8e0eb19d8a7" [[package]] name = "const-primes" -version = "0.9.6" +version = "0.9.7" dependencies = [ "criterion", + "machine-prime", "rand", "rkyv", "serde", @@ -248,6 +249,12 @@ version = "0.4.22" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a7a70ba024b9dc04c27ea2f0c0548feb474ec5c54bba33a7f72f873a39d07b24" +[[package]] +name = "machine-prime" +version = "1.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "31ced57f9093501cfed04a9deca0dcb6f36475fa8f00495a6d99823a75d9f714" + [[package]] name = "memchr" version = "2.7.4" diff --git a/Cargo.toml b/Cargo.toml index b13d5de..5d663bb 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "const-primes" authors = ["Johanna Sörngård "] -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 diff --git a/README.md b/README.md index fca4c45..ec9e097 100644 --- a/README.md +++ b/README.md @@ -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. +
### License diff --git a/src/check.rs b/src/check.rs index 60be5f7..e72b5e2 100644 --- a/src/check.rs +++ b/src/check.rs @@ -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 - // . - 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 + // . + 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 { diff --git a/src/integer_math.rs b/src/integer_math.rs index 665f02a..f4a959a 100644 --- a/src/integer_math.rs +++ b/src/integer_math.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 } diff --git a/src/lib.rs b/src/lib.rs index 0ed70f6..6996fc8 100644 --- a/src/lib.rs +++ b/src/lib.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] diff --git a/src/search.rs b/src/search.rs index 68e0711..cb54006 100644 --- a/src/search.rs +++ b/src/search.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 { + 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 = 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 { - 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 { + 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 { /// ``` /// /// Primes larger than 18446744073709551557 can not be represented by a `u64`: -/// ``` /// +/// ``` /// # use const_primes::next_prime; /// const NO_SUCH: Option = 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 { - // 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 { + 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) } }