FIPS compatibility

CHANGELOG.md

@@ -11,6 +11,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Bumped `rand-core` to 0.9. ([#74])
 
 
+### Added
+
+- `LucasCheck::Regular` (the recommended one from the FIPS-186.5 standard). ([#72])
+- `hamzat::minimum_mr_iterations()` (to calculate the number of MR checks according to the FIPS-186.5 standard). ([#72])
+- Add `fips_is_prime_with_rng()` and `fips_is_safe_prime_with_rng()`. ([#72])
+
+
+[#72]: https://github.com/entropyxyz/crypto-primes/pull/72
 [#74]: https://github.com/entropyxyz/crypto-primes/pull/74
 
 

Cargo.toml

@@ -10,7 +10,7 @@ categories = ["cryptography", "no-std"]
 rust-version = "1.83"
 
 [dependencies]
-crypto-bigint = { version = "0.7.0-pre", default-features = false, features = ["rand_core"] }
+crypto-bigint = { version = "0.7.0-pre.0", default-features = false, features = ["rand_core"] }
 rand_core = { version = "0.9.0", default-features = false }
 
 # Optional dependencies used in tests and benchmarks
@@ -21,7 +21,7 @@ rayon = { version = "1", optional = true }
 
 [dev-dependencies]
 # need `crypto-bigint` with `alloc` to test `BoxedUint`
-crypto-bigint = { version = "0.7.0-pre", default-features = false, features = ["alloc"] }
+crypto-bigint = { version = "0.7.0-pre.0", default-features = false, features = ["alloc"] }
 rand_chacha = "0.9"
 criterion = { version = "0.5", features = ["html_reports"] }
 num-modular = { version = "0.5", features = ["num-bigint"] }
@@ -30,6 +30,11 @@ num-integer = "0.1"
 proptest = "1"
 num-prime = "0.4.3"
 num_cpus = "1.16"
+float-cmp = "0.10"
+
+# Temporary old versions for `glass_pumpkin` tests. Remove when `glass_pumpking` switches to `rand_core=0.9`.
+rand_core_06 = { package = "rand_core", version = "0.6.4", default-features = false }
+rand_chacha_03 = { package = "rand_chacha", version = "0.3", default-features = false }
 
 [features]
 default = ["default-rng"]
@@ -53,13 +58,3 @@ harness = false
 [[bench]]
 name = "cctv"
 harness = false
-
-[patch.crates-io]
-# https://github.com/RustCrypto/crypto-bigint/pull/762
-# https://github.com/RustCrypto/crypto-bigint/pull/765
-crypto-bigint = { git = "https://github.com/RustCrypto/crypto-bigint.git" }
-
-# https://github.com/LF-Decentralized-Trust-labs/agora-glass_pumpkin/pull/26
-glass_pumpkin = { git = "https://github.com/baloo/agora-glass_pumpkin.git", branch = "baloo/rand-core/0.9" }
-# https://github.com/rust-num/num-bigint/pull/317
-num-bigint = { git = "https://github.com/bionicles/num-bigint.git" }

benches/bench.rs

@@ -435,6 +435,10 @@ fn bench_glass_pumpkin(c: &mut Criterion) {
     use crypto_bigint::Limb;
     use crypto_primes::hazmat::{lucas_test, AStarBase, LucasCheck, MillerRabin, Primality};
 
+    fn make_rng_gp() -> rand_chacha_03::ChaCha8Rng {
+        <rand_chacha_03::ChaCha8Rng as rand_core_06::SeedableRng>::from_seed(*b"01234567890123456789012345678901")
+    }
+
     // The `glass-pumpkin` implementation is doing a different number of M-R checks than this crate.
     // For a fair comparison we make a custom implementation here,
     // using the same number of checks that `glass-pumpkin` does.
@@ -527,9 +531,9 @@ fn bench_glass_pumpkin(c: &mut Criterion) {
         b.iter(|| prime_like_gp(1024, &mut rng))
     });
 
-    let mut rng = make_rng();
+    let mut rng_gp = make_rng_gp();
     group.bench_function("(U1024) Random prime", |b| {
-        b.iter(|| glass_pumpkin::prime::from_rng(1024, &mut rng))
+        b.iter(|| glass_pumpkin::prime::from_rng(1024, &mut rng_gp))
     });
 
     group.sample_size(20);
@@ -545,9 +549,9 @@ fn bench_glass_pumpkin(c: &mut Criterion) {
         |b| b.iter(|| safe_prime_like_gp(1024, &mut rng)),
     );
 
-    let mut rng = make_rng();
+    let mut rng_gp = make_rng_gp();
     group.bench_function("(U1024) Random safe prime", |b| {
-        b.iter(|| glass_pumpkin::safe_prime::from_rng(1024, &mut rng))
+        b.iter(|| glass_pumpkin::safe_prime::from_rng(1024, &mut rng_gp))
     });
 }
 

src/hazmat.rs

@@ -1,6 +1,7 @@
 //! Components to build your own primality test.
 //! Handle with care.
 
+mod float;
 mod gcd;
 mod jacobi;
 mod lucas;
@@ -13,7 +14,7 @@ pub(crate) mod pseudoprimes;
 mod sieve;
 
 pub use lucas::{lucas_test, AStarBase, BruteForceBase, LucasBase, LucasCheck, SelfridgeBase};
-pub use miller_rabin::MillerRabin;
+pub use miller_rabin::{minimum_mr_iterations, MillerRabin};
 pub use sieve::{random_odd_integer, SetBits, SmallPrimesSieve, SmallPrimesSieveFactory};
 
 /// Possible results of various primality tests.

src/hazmat/float.rs

@@ -0,0 +1,156 @@
+//! Const-context floating point functions that are currently not present in `core`.
+
+/// Calculates `base^exp`.
+const fn pow(mut base: f64, mut exp: u32) -> f64 {
+    let mut result = 1.;
+    while exp > 0 {
+        if exp & 1 == 1 {
+            result *= base;
+        }
+        base *= base;
+        exp >>= 1;
+    }
+    result
+}
+
+/// Calculates `2^exp`.
+pub(crate) const fn two_powi(exp: u32) -> f64 {
+    pow(2f64, exp.abs_diff(0))
+}
+
+/// Calculates `floor(x)`.
+// Taken from `libm` crate.
+const fn floor(x: f64) -> f64 {
+    const TOINT: f64 = 1. / f64::EPSILON;
+
+    let ui = x.to_bits();
+    let e = ((ui >> 52) & 0x7ff) as i32;
+
+    if (e >= 0x3ff + 52) || (x == 0.) {
+        return x;
+    }
+    /* y = int(x) - x, where int(x) is an integer neighbor of x */
+    let y = if (ui >> 63) != 0 {
+        x - TOINT + TOINT - x
+    } else {
+        x + TOINT - TOINT - x
+    };
+    /* special case because of non-nearest rounding modes */
+    if e < 0x3ff {
+        return if (ui >> 63) != 0 { -1. } else { 0. };
+    }
+    if y > 0. {
+        x + y - 1.
+    } else {
+        x + y
+    }
+}
+
+/// Calculates a lower bound approximation of `2^exp` where `0 <= exp <= 1`.
+const fn two_powf_normalized_lower_bound(exp: f64) -> f64 {
+    debug_assert!(exp >= 0.);
+    debug_assert!(exp <= 1.);
+
+    // Use the first four terms of the Taylor expansion to calculate `2^exp`.
+    // The error is under 2%.
+    //
+    // Since it is a monotonous function, `res <= 2^exp`.
+
+    let exp_2 = exp * exp;
+    let exp_3 = exp_2 * exp;
+
+    // The coefficients are `ln(2)^n / n!`, where `n` is the power of the corresponding term.
+    const LN_2: f64 = core::f64::consts::LN_2;
+    const C1: f64 = LN_2;
+    const C2: f64 = LN_2 * LN_2 / 2.;
+    const C3: f64 = LN_2 * LN_2 * LN_2 / 6.;
+
+    1. + C1 * exp + C2 * exp_2 + C3 * exp_3
+}
+
+/// Calculates an approximation of `2^exp` where `exp < 0`.
+/// The approximation is guaranteed to always be greater than `2^exp`.
+pub(crate) const fn two_powf_upper_bound(exp: f64) -> f64 {
+    debug_assert!(exp < 0.);
+
+    let positive_exp = -exp;
+
+    let int_part = floor(positive_exp);
+    let frac_part = positive_exp - int_part;
+
+    let int_res = two_powi(int_part as u32);
+    let frac_res = two_powf_normalized_lower_bound(frac_part);
+
+    // `int_res * frac_res <= 2^(int_part + frac_part)`,
+    // so when we invert it, we get the upper bound approximation instead.
+    1. / (int_res * frac_res)
+}
+
+/// Calculates `floor(sqrt(x))`.
+pub(crate) const fn floor_sqrt(x: u32) -> u32 {
+    if x < 2 {
+        return x;
+    }
+
+    // Initialize the binary search bounds.
+    let mut low = 1;
+    let mut high = x / 2;
+
+    while low <= high {
+        let mid = (low + high) / 2;
+        let mid_squared = mid * mid;
+
+        // Check if `mid` is the floor of the square root of `x`.
+        if mid_squared <= x && (mid + 1) * (mid + 1) > x {
+            break;
+        } else if mid_squared < x {
+            low = mid + 1;
+        } else {
+            high = mid - 1
+        }
+    }
+
+    (low + high) / 2
+}
+
+#[cfg(test)]
+mod tests {
+    use float_cmp::assert_approx_eq;
+    use proptest::prelude::*;
+
+    use super::{floor, floor_sqrt, pow, two_powf_normalized_lower_bound};
+
+    proptest! {
+        #[test]
+        fn fuzzy_pow(base in 0..100u32, exp in 0..30u32) {
+            let base_f = base as f64 / 100.;
+            let test = pow(base_f, exp);
+            let reference = base_f.powf(exp as f64);
+            assert_approx_eq!(f64, test, reference, ulps = 20);
+        }
+
+        #[test]
+        fn fuzzy_floor(x in proptest::num::f64::NORMAL) {
+            let test = floor(x);
+            let reference = x.floor();
+            assert_approx_eq!(f64, test, reference);
+        }
+
+        #[test]
+        fn fuzzy_two_powf_upper_bound(exp in 0..1000) {
+            let exp_f = exp as f64 / 1000.;
+            let test = two_powf_normalized_lower_bound(exp_f);
+            let reference = 2f64.powf(exp_f);
+            assert!(test <= reference);
+            assert!((reference - test) / reference <= 0.02);
+        }
+
+        #[test]
+        fn fuzzy_floor_sqrt(x in 0..100000u32) {
+            let x_f = x as f64;
+            let test = floor_sqrt(x);
+            let reference = x_f.sqrt().floor() as u32;
+            assert_eq!(test, reference);
+        }
+    }
+}