diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4390bdb..49351d8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,18 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [0.7.0] - in development
+
+### 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
+
+
 ## [0.7.0-pre.0] - 2025-02-22
 
 ### Changed
diff --git a/Cargo.toml b/Cargo.toml
index 9438273..26a51cb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "crypto-primes"
-version = "0.7.0-pre.0"
+version = "0.7.0-dev"
 edition = "2021"
 license = "Apache-2.0 OR MIT"
 description = "Random prime number generation and primality checking library"
@@ -30,6 +30,7 @@ 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 }
diff --git a/benches/bench.rs b/benches/bench.rs
index 14ee204..f1139a4 100644
--- a/benches/bench.rs
+++ b/benches/bench.rs
@@ -12,10 +12,10 @@ use rug::{integer::Order, Integer as GmpInteger};
 use openssl::bn::BigNum;
 
 use crypto_primes::{
-    generate_prime_with_rng, generate_safe_prime_with_rng,
+    fips_is_prime_with_rng, generate_prime_with_rng, generate_safe_prime_with_rng,
     hazmat::{
-        lucas_test, random_odd_integer, AStarBase, BruteForceBase, LucasCheck, MillerRabin, SelfridgeBase, SetBits,
-        SmallPrimesSieve,
+        lucas_test, minimum_mr_iterations, random_odd_integer, AStarBase, BruteForceBase, LucasCheck, MillerRabin,
+        SelfridgeBase, SetBits, SmallPrimesSieve,
     },
     is_prime_with_rng, is_safe_prime_with_rng,
 };
@@ -221,6 +221,23 @@ fn bench_presets(c: &mut Criterion) {
         )
     });
 
+    group.bench_function("(U1024) Prime test", |b| {
+        b.iter_batched(
+            || random_odd_uint::<U1024, _>(&mut OsRng.unwrap_err(), 1024),
+            |num| is_prime_with_rng(&mut OsRng.unwrap_err(), num.as_ref()),
+            BatchSize::SmallInput,
+        )
+    });
+
+    let iters = minimum_mr_iterations(1024, 128).unwrap();
+    group.bench_function("(U1024) Prime test (FIPS, 1/2^128 failure bound)", |b| {
+        b.iter_batched(
+            || random_odd_uint::<U1024, _>(&mut OsRng.unwrap_err(), 1024),
+            |num| fips_is_prime_with_rng(&mut OsRng.unwrap_err(), num.as_ref(), iters, false),
+            BatchSize::SmallInput,
+        )
+    });
+
     group.bench_function("(U128) Safe prime test", |b| {
         b.iter_batched(
             || random_odd_uint::<U128, _>(&mut OsRng.unwrap_err(), 128),
diff --git a/src/hazmat.rs b/src/hazmat.rs
index 9e397d4..53f88f3 100644
--- a/src/hazmat.rs
+++ b/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.
diff --git a/src/hazmat/float.rs b/src/hazmat/float.rs
new file mode 100644
index 0000000..799d1a6
--- /dev/null
+++ b/src/hazmat/float.rs
@@ -0,0 +1,163 @@
+//! 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};
+
+    #[test]
+    fn sqrt_corner_cases() {
+        assert_eq!(floor_sqrt(0), 0);
+        assert_eq!(floor_sqrt(1), 1);
+        assert_eq!(floor_sqrt(2), 1);
+    }
+
+    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);
+        }
+    }
+}
diff --git a/src/hazmat/lucas.rs b/src/hazmat/lucas.rs
index 15784a5..dbeed29 100644
--- a/src/hazmat/lucas.rs
+++ b/src/hazmat/lucas.rs
@@ -208,6 +208,24 @@ fn decompose<T: Integer>(n: &Odd<T>) -> (u32, Odd<T>) {
 /// the checks are defined as follows:
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub enum LucasCheck {
+    /// Introduced by Baillie & Wagstaff[^Baillie1980].
+    /// If `U(n - (D/n)) == 0`, report the number as prime.
+    ///
+    /// This is the Lucas test prescribed by the FIPS-186.5[^FIPS] standard (Section B.3.3).
+    ///
+    /// If the base is [`SelfridgeBase`], known false positives constitute OEIS:A217120[^A217120].
+    ///
+    /// [^Baillie1980]:
+    ///   R. Baillie, S. S. Wagstaff, "Lucas pseudoprimes",
+    ///   Math. Comp. 35 1391-1417 (1980),
+    ///   DOI: [10.2307/2006406](https://dx.doi.org/10.2307/2006406),
+    ///   <http://mpqs.free.fr/LucasPseudoprimes.pdf>
+    ///
+    /// [^FIPS]: <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
+    ///
+    /// [^A217120]: <https://oeis.org/A217120>
+    Regular,
+
     /// Introduced by Baillie & Wagstaff[^Baillie1980].
     /// If either of the following is true:
     /// - any of `V(d*2^r) == 0` for `0 <= r < s`,
@@ -285,8 +303,14 @@ pub enum LucasCheck {
 }
 
 /// Performs the primality test based on Lucas sequence.
+///
 /// See [`LucasCheck`] for possible checks, and the implementors of [`LucasBase`]
 /// for the corresponding bases.
+///
+/// When used with [`SelfridgeBase`] and [`LucasCheck::Regular`], implements the algorithm
+/// prescribed by the FIPS.186-5 standard[^FIPS].
+///
+/// [^FIPS]: FIPS-186.5 standard, <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
 pub fn lucas_test<T: Integer>(candidate: Odd<T>, base: impl LucasBase, check: LucasCheck) -> Primality {
     // The comments in this function use references in `LucasCheck`, plus this one:
     //
@@ -402,7 +426,7 @@ pub fn lucas_test<T: Integer>(candidate: Odd<T>, base: impl LucasBase, check: Lu
         // k' = k * 2
 
         let u_2k = uk * &vk;
-        let v_2k = vk.square() - &(qk.clone() + &qk);
+        let v_2k = vk.square() - &qk.double();
         let q_2k = qk.square();
 
         uk = u_2k;
@@ -430,57 +454,71 @@ pub fn lucas_test<T: Integer>(candidate: Odd<T>, base: impl LucasBase, check: Lu
 
     // Now k=d, so vk = V_d and uk = U_d.
 
-    // Check for the first sufficient condition in various strong checks.
+    // The `U_d == 0` criterion.
+    let ud_equals_zero = uk == zero;
 
-    if check == LucasCheck::Strong && uk == zero {
-        // Strong check: `U_d == 0 mod n`.
+    // The `V_d == ±2 mod n` criterion.
+    //
+    // Note that this criterion only applies if `Q = 1`, since it is a consequence
+    // of a property of Lucas series: `V_k^2 - 4 Q^k = D U_k^2 mod n`.
+    // If `Q = 1` we can easily decompose the left side of the equation
+    // leading to the check above.
+    //
+    // If `Q != 1` we just consider it passed (we don't have a corresponding
+    // pseudoprime list anyway).
+    let vk_equals_two = !q_is_one || (vk == two || vk == minus_two);
+
+    // Early exit for some of the checks.
+    // If the conditions are not satisfied, we have to continue propagating the Lucas sequence.
+
+    if check == LucasCheck::Strong && ud_equals_zero {
         return Primality::ProbablyPrime;
-    } else if check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong {
-        // Extra strong check (from [^Mo1993]): `V_d == ±2 mod n` and `U_d == 0 mod n`.
-        //
-        // Note that the first identity only applies if `Q = 1`, since it is a consequence
-        // of a property of Lucas series: `V_k^2 - 4 Q^k = D U_k^2 mod n`.
-        // If `Q = 1` we can easily decompose the left side of the equation
-        // leading to the check above.
-        //
-        // If `Q != 1` we just consider it passed (we don't have a corresponding
-        // pseudoprime list anyway).
-
-        let vk_equals_two = !q_is_one || (vk == two || vk == minus_two);
-
-        if check == LucasCheck::ExtraStrong && uk == zero && vk_equals_two {
-            return Primality::ProbablyPrime;
-        }
+    }
 
-        // "Almost extra strong" check skips the `U_d` check.
-        // Since we have `U_d` anyway, it does not improve performance,
-        // so it is only here for testing purposes, since we have a corresponding pseudoprime list.
-        if check == LucasCheck::AlmostExtraStrong && vk_equals_two {
-            return Primality::ProbablyPrime;
-        }
+    if check == LucasCheck::ExtraStrong && ud_equals_zero && vk_equals_two {
+        return Primality::ProbablyPrime;
     }
 
-    // Second sufficient condition requires further propagating `V_k` up to `V_{n+1}`.
+    // "Almost extra strong" check skips the `U_d` check.
+    // Since we have `U_d` anyway, it does not improve performance,
+    // so it is only here for testing purposes, since we have a corresponding pseudoprime list.
+    if check == LucasCheck::AlmostExtraStrong && vk_equals_two {
+        return Primality::ProbablyPrime;
+    }
+
+    // Propagate the Lucas sequence from `d` to `n+1`, exiting early for some of the checks.
 
-    // Check if V_{2^t d} == 0 mod n for some 0 <= t < s.
-    // (unless we're in Lucas-V mode, then we just propagate V_k)
+    let mut one_of_vk_equals_zero = vk == zero;
 
-    if check != LucasCheck::LucasV && vk == zero {
+    if (check == LucasCheck::Strong || check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong)
+        && one_of_vk_equals_zero
+    {
         return Primality::ProbablyPrime;
     }
 
     for _ in 1..s {
         // Optimization: V_k = ±2 is a fixed point for V_k' = V_k^2 - 2 Q^k with Q = 1,
         // so if V_k = ±2, we can stop: we will never find a future V_k == 0.
-        if check != LucasCheck::LucasV && q_is_one && (vk == two || vk == minus_two) {
+        if (check == LucasCheck::Strong || check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong)
+            && q_is_one
+            && (vk == two || vk == minus_two)
+        {
             return Primality::Composite;
         }
 
+        if check == LucasCheck::Regular {
+            uk *= &vk;
+        }
+
         // k' = 2k
         // V_{k'} = V_k^2 - 2 Q^k
-        vk = vk.square() - &qk - &qk;
+        vk = vk.square() - &qk.double();
+
+        one_of_vk_equals_zero |= vk == zero;
 
-        if check != LucasCheck::LucasV && vk == zero {
+        if (check == LucasCheck::Strong || check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong)
+            && one_of_vk_equals_zero
+        {
             return Primality::ProbablyPrime;
         }
 
@@ -489,13 +527,34 @@ pub fn lucas_test<T: Integer>(candidate: Odd<T>, base: impl LucasBase, check: Lu
         }
     }
 
+    if check == LucasCheck::Strong || check == LucasCheck::ExtraStrong || check == LucasCheck::AlmostExtraStrong {
+        return Primality::Composite;
+    }
+
+    // At this point:
+    //   vk = V_{d * 2^(s-1)} == V_{(n + 1) / 2}.
+    //   qk = Q^{(n + 1) / 2}
+    // In case of `check == Regular`, also
+    //   uk = U_{(n + 1) / 2}
+
+    if check == LucasCheck::Regular {
+        // Double the index again:
+        uk *= &vk; // now `uk = U_{d * 2^s} = U_{n+1}`
+        if uk == zero {
+            return Primality::ProbablyPrime;
+        } else {
+            return Primality::Composite;
+        }
+    }
+
+    // Lucas-V check[^Baillie2021]: if `V_{n+1} != 2 Q`, report `n` as composite.
     if check == LucasCheck::LucasV {
-        // At this point vk = V_{d * 2^(s-1)}.
         // Double the index again:
-        vk = vk.square() - &qk - &qk; // now vk = V_{d * 2^s} = V_{n+1}
+        vk = vk.square() - &qk.double(); // now `vk = V_{d * 2^s} = V_{n+1}`
 
-        // Lucas-V check[^Baillie2021]: if V_{n+1} == 2 Q, report `n` as prime.
-        if vk == q.clone() + &q {
+        if vk != q.double() {
+            return Primality::Composite;
+        } else {
             return Primality::ProbablyPrime;
         }
     }
@@ -586,96 +645,71 @@ mod tests {
         );
     }
 
-    fn is_slpsp(num: u32) -> bool {
-        pseudoprimes::STRONG_LUCAS.iter().any(|x| *x == num)
+    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
+    enum BaseType {
+        Selfridge,
+        BruteForce,
+        AStar,
     }
 
-    fn is_aeslpsp(num: u32) -> bool {
-        pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS.iter().any(|x| *x == num)
+    trait HasBaseType: LucasBase + Copy {
+        const BASE_TYPE: BaseType;
     }
 
-    fn is_eslpsp(num: u32) -> bool {
-        pseudoprimes::EXTRA_STRONG_LUCAS.iter().any(|x| *x == num)
+    impl HasBaseType for SelfridgeBase {
+        const BASE_TYPE: BaseType = BaseType::Selfridge;
     }
 
-    fn is_vpsp(num: u32) -> bool {
-        pseudoprimes::LUCAS_V.iter().any(|x| *x == num)
+    impl HasBaseType for BruteForceBase {
+        const BASE_TYPE: BaseType = BaseType::BruteForce;
     }
 
-    fn test_composites_selfridge(numbers: &[u32], expected_result: bool) {
-        for num in numbers.iter() {
-            let false_positive = is_slpsp(*num);
-            let actual_expected_result = if false_positive { true } else { expected_result };
-
-            // Test both single-limb and multi-limb, just in case.
-            assert_eq!(
-                lucas_test(
-                    Odd::new(Uint::<1>::from(*num)).unwrap(),
-                    SelfridgeBase,
-                    LucasCheck::Strong
-                )
-                .is_probably_prime(),
-                actual_expected_result
-            );
-            assert_eq!(
-                lucas_test(
-                    Odd::new(Uint::<2>::from(*num)).unwrap(),
-                    SelfridgeBase,
-                    LucasCheck::Strong
-                )
-                .is_probably_prime(),
-                actual_expected_result
-            );
-        }
+    impl HasBaseType for AStarBase {
+        const BASE_TYPE: BaseType = BaseType::AStar;
     }
 
-    fn test_composites_a_star(numbers: &[u32], expected_result: bool) {
-        for num in numbers.iter() {
-            let false_positive = is_vpsp(*num);
-            let actual_expected_result = if false_positive { true } else { expected_result };
+    // Returns `true` if `num` is a composite that is known to be reported as a prime
+    // by a Lucas test with the given base `T` and check type `check`.
+    // Returns `false` is `num` is not in the list of such composites for the given base and check type.
+    //
+    // Panics if there is no data for the given base and check type.
+    fn is_pseudoprime<T: HasBaseType>(num: u32, check: LucasCheck) -> bool {
+        let pseudoprimes = match (T::BASE_TYPE, check) {
+            (BaseType::Selfridge, LucasCheck::Regular) => &pseudoprimes::LUCAS,
+            (BaseType::Selfridge, LucasCheck::Strong) => &pseudoprimes::STRONG_LUCAS,
+            (BaseType::BruteForce, LucasCheck::AlmostExtraStrong) => &pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS,
+            (BaseType::BruteForce, LucasCheck::ExtraStrong) => &pseudoprimes::EXTRA_STRONG_LUCAS,
+            (BaseType::AStar, LucasCheck::LucasV) => &pseudoprimes::LUCAS_V,
+            _ => panic!("We do not have pseudoprimes listed for this combination of base and check"),
+        };
 
-            // Test both single-limb and multi-limb, just in case.
-            assert_eq!(
-                lucas_test(Odd::new(Uint::<1>::from(*num)).unwrap(), AStarBase, LucasCheck::LucasV).is_probably_prime(),
-                actual_expected_result
-            );
-            assert_eq!(
-                lucas_test(Odd::new(Uint::<2>::from(*num)).unwrap(), AStarBase, LucasCheck::LucasV).is_probably_prime(),
-                actual_expected_result
-            );
-        }
+        pseudoprimes.iter().any(|x| *x == num)
     }
 
-    fn test_composites_brute_force(numbers: &[u32], almost_extra: bool, expected_result: bool) {
+    fn test_pseudoprimes<T: HasBaseType>(numbers: &[u32], base: T, check: LucasCheck, expected_result: bool) {
         for num in numbers.iter() {
-            let false_positive = if almost_extra {
-                is_aeslpsp(*num)
-            } else {
-                is_eslpsp(*num)
-            };
+            let false_positive = is_pseudoprime::<T>(*num, check);
             let actual_expected_result = if false_positive { true } else { expected_result };
-            let check = if almost_extra {
-                LucasCheck::AlmostExtraStrong
-            } else {
-                LucasCheck::ExtraStrong
-            };
 
             // Test both single-limb and multi-limb, just in case.
+
+            let is_prime = lucas_test(Odd::new(Uint::<1>::from(*num)).unwrap(), base, check).is_probably_prime();
             assert_eq!(
-                lucas_test(Odd::new(Uint::<1>::from(*num)).unwrap(), BruteForceBase, check).is_probably_prime(),
-                actual_expected_result,
-                "Brute force base, n = {num}, almost_extra = {almost_extra}",
+                is_prime, actual_expected_result,
+                "{num} reported as prime={is_prime}, expected prime={actual_expected_result}"
             );
+
+            let is_prime = lucas_test(Odd::new(Uint::<2>::from(*num)).unwrap(), base, check).is_probably_prime();
             assert_eq!(
-                lucas_test(Odd::new(Uint::<2>::from(*num)).unwrap(), BruteForceBase, check).is_probably_prime(),
-                actual_expected_result
+                is_prime, actual_expected_result,
+                "{num} reported as prime={is_prime}, expected prime={actual_expected_result}"
             );
         }
     }
 
     #[test]
     fn strong_fibonacci_pseudoprimes() {
-        // Can't use `test_composites()` since `STRONG_FIBONACCI` is `U64`.
+        // Can't use `test_pseudoprimes()` since `STRONG_FIBONACCI` is `U64`.
         // Good thing we don't need to test for intersection
         // with `EXTRA_STRONG_LUCAS` or `STRONG_LUCAS` - there's none.
         for num in pseudoprimes::STRONG_FIBONACCI.iter() {
@@ -686,72 +720,84 @@ mod tests {
 
     #[test]
     fn fibonacci_pseudoprimes() {
-        test_composites_selfridge(pseudoprimes::FIBONACCI, false);
-        test_composites_a_star(pseudoprimes::FIBONACCI, false);
-        test_composites_brute_force(pseudoprimes::FIBONACCI, false, false);
-        test_composites_brute_force(pseudoprimes::FIBONACCI, true, false);
+        let nums = pseudoprimes::FIBONACCI;
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::AlmostExtraStrong, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, false);
     }
 
     #[test]
     fn bruckman_lucas_pseudoprimes() {
-        test_composites_selfridge(pseudoprimes::BRUCKMAN_LUCAS, false);
-        test_composites_a_star(pseudoprimes::BRUCKMAN_LUCAS, false);
-        test_composites_brute_force(pseudoprimes::BRUCKMAN_LUCAS, false, false);
-        test_composites_brute_force(pseudoprimes::BRUCKMAN_LUCAS, true, false);
+        let nums = pseudoprimes::BRUCKMAN_LUCAS;
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::AlmostExtraStrong, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, false);
     }
 
     #[test]
     fn almost_extra_strong_lucas_pseudoprimes() {
-        test_composites_selfridge(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false);
-        test_composites_a_star(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false);
+        let nums = pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS;
+
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
 
         // Check for the difference between the almost extra strong and extra strong tests.
-        test_composites_brute_force(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false, false);
-        test_composites_brute_force(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, true, true);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::AlmostExtraStrong, true);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, false);
     }
 
     #[test]
     fn extra_strong_lucas_pseudoprimes() {
-        test_composites_selfridge(pseudoprimes::EXTRA_STRONG_LUCAS, false);
-        test_composites_a_star(pseudoprimes::EXTRA_STRONG_LUCAS, false);
+        let nums = pseudoprimes::EXTRA_STRONG_LUCAS;
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
 
         // These are the known false positives for the extra strong test
         // with brute force base selection.
-        test_composites_brute_force(pseudoprimes::EXTRA_STRONG_LUCAS, false, true);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, true);
     }
 
     #[test]
     fn lucas_pseudoprimes() {
-        test_composites_selfridge(pseudoprimes::LUCAS, false);
-        test_composites_a_star(pseudoprimes::LUCAS, false);
-        test_composites_brute_force(pseudoprimes::LUCAS, false, false);
-        test_composites_brute_force(pseudoprimes::LUCAS, true, false);
+        let nums = pseudoprimes::LUCAS;
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Regular, true);
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::AlmostExtraStrong, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, false);
     }
 
     #[test]
     fn strong_lucas_pseudoprimes() {
+        let nums = pseudoprimes::STRONG_LUCAS;
+
         // These are the known false positives for the strong test
         // with Selfridge base selection.
-        test_composites_selfridge(pseudoprimes::STRONG_LUCAS, true);
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, true);
 
-        test_composites_a_star(pseudoprimes::STRONG_LUCAS, false);
-        test_composites_brute_force(pseudoprimes::STRONG_LUCAS, false, false);
-        test_composites_brute_force(pseudoprimes::STRONG_LUCAS, true, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::AlmostExtraStrong, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, false);
     }
 
     #[test]
     fn strong_pseudoprimes_base_2() {
         // Cross-test against the pseudoprimes that circumvent the MR test base 2.
         // We expect the Lucas test to correctly classify them as composites.
-        test_composites_selfridge(pseudoprimes::STRONG_BASE_2, false);
-        test_composites_a_star(pseudoprimes::STRONG_BASE_2, false);
-        test_composites_brute_force(pseudoprimes::STRONG_BASE_2, false, false);
-        test_composites_brute_force(pseudoprimes::STRONG_BASE_2, true, false);
+
+        let nums = pseudoprimes::STRONG_BASE_2;
+        test_pseudoprimes(nums, SelfridgeBase, LucasCheck::Strong, false);
+        test_pseudoprimes(nums, AStarBase, LucasCheck::LucasV, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::AlmostExtraStrong, false);
+        test_pseudoprimes(nums, BruteForceBase, LucasCheck::ExtraStrong, false);
     }
 
     #[test]
     fn large_carmichael_number() {
         let p = Odd::new(pseudoprimes::LARGE_CARMICHAEL_NUMBER).unwrap();
+        assert!(!lucas_test(p, SelfridgeBase, LucasCheck::Regular).is_probably_prime());
         assert!(!lucas_test(p, SelfridgeBase, LucasCheck::Strong).is_probably_prime());
         assert!(!lucas_test(p, AStarBase, LucasCheck::LucasV).is_probably_prime());
         assert!(!lucas_test(p, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime());
@@ -761,6 +807,7 @@ mod tests {
     fn test_large_primes<const L: usize>(nums: &[Uint<L>]) {
         for num in nums {
             let num = Odd::new(*num).unwrap();
+            assert!(lucas_test(num, SelfridgeBase, LucasCheck::Regular).is_probably_prime());
             assert!(lucas_test(num, SelfridgeBase, LucasCheck::Strong).is_probably_prime());
             assert!(lucas_test(num, AStarBase, LucasCheck::LucasV).is_probably_prime());
             assert!(lucas_test(num, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime());
@@ -810,13 +857,17 @@ mod tests {
         for num in (3..pseudoprimes::EXHAUSTIVE_TEST_LIMIT).step_by(2) {
             let res_ref = is_prime64(num.into());
 
-            let eslpsp = is_eslpsp(num);
-            let aeslpsp = is_aeslpsp(num);
-            let slpsp = is_slpsp(num);
-            let vpsp = is_vpsp(num);
-
             let odd_num = Odd::new(Uint::<1>::from(num)).unwrap();
 
+            let lpsp = is_pseudoprime::<SelfridgeBase>(num, LucasCheck::Regular);
+            let res = lucas_test(odd_num, SelfridgeBase, LucasCheck::Regular).is_probably_prime();
+            let expected = lpsp || res_ref;
+            assert_eq!(
+                res, expected,
+                "Selfridge base, regular: n={num}, expected={expected}, actual={res}",
+            );
+
+            let aeslpsp = is_pseudoprime::<BruteForceBase>(num, LucasCheck::AlmostExtraStrong);
             let res = lucas_test(odd_num, BruteForceBase, LucasCheck::AlmostExtraStrong).is_probably_prime();
             let expected = aeslpsp || res_ref;
             assert_eq!(
@@ -824,6 +875,7 @@ mod tests {
                 "Brute force base, almost extra strong: n={num}, expected={expected}, actual={res}",
             );
 
+            let eslpsp = is_pseudoprime::<BruteForceBase>(num, LucasCheck::ExtraStrong);
             let res = lucas_test(odd_num, BruteForceBase, LucasCheck::ExtraStrong).is_probably_prime();
             let expected = eslpsp || res_ref;
             assert_eq!(
@@ -831,6 +883,7 @@ mod tests {
                 "Brute force base: n={num}, expected={expected}, actual={res}",
             );
 
+            let slpsp = is_pseudoprime::<SelfridgeBase>(num, LucasCheck::Strong);
             let res = lucas_test(odd_num, SelfridgeBase, LucasCheck::Strong).is_probably_prime();
             let expected = slpsp || res_ref;
             assert_eq!(
@@ -838,6 +891,7 @@ mod tests {
                 "Selfridge base: n={num}, expected={expected}, actual={res}",
             );
 
+            let vpsp = is_pseudoprime::<AStarBase>(num, LucasCheck::LucasV);
             let res = lucas_test(odd_num, AStarBase, LucasCheck::LucasV).is_probably_prime();
             let expected = vpsp || res_ref;
 
diff --git a/src/hazmat/miller_rabin.rs b/src/hazmat/miller_rabin.rs
index 78dba30..da6add1 100644
--- a/src/hazmat/miller_rabin.rs
+++ b/src/hazmat/miller_rabin.rs
@@ -3,17 +3,23 @@
 use crypto_bigint::{Integer, Limb, Monty, NonZero as CTNonZero, Odd, PowBoundedExp, RandomMod, Square};
 use rand_core::CryptoRng;
 
-use super::Primality;
+use super::{
+    float::{floor_sqrt, two_powf_upper_bound, two_powi},
+    Primality,
+};
 
 /// Precomputed data used to perform Miller-Rabin primality test[^Pomerance1980].
 ///
 /// The numbers that pass it are commonly called "strong probable primes"
 /// (or "strong pseudoprimes" if they are, in fact, composite).
 ///
+/// The implementation satisfies the FIPS.186-5 standard[^FIPS].
+///
 /// [^Pomerance1980]:
 ///   C. Pomerance, J. L. Selfridge, S. S. Wagstaff "The Pseudoprimes to 25*10^9",
 ///   Math. Comp. 35 1003-1026 (1980),
 ///   DOI: [10.2307/2006210](https://dx.doi.org/10.2307/2006210)
+/// [^FIPS]: FIPS-186.5 standard, <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
 #[derive(Clone, Debug, PartialEq, Eq)]
 pub struct MillerRabin<T: Integer> {
     // The odd number that may or may not be a prime.
@@ -123,15 +129,98 @@ impl<T: Integer + RandomMod> MillerRabin<T> {
             .expect("addition should not overflow by construction");
         self.test(&random)
     }
+}
 
-    /// Returns the number of bits necessary to represent the candidate.
-    /// NOTE: This is different than the number of bits of *storage* the integer takes up.
-    ///
-    /// For example, a U512 type occupies 8 64-bit words, but the number `7` contained in such a type
-    /// has a bit length of 3 because 7 is `b111`.
-    pub(crate) fn bits(&self) -> u32 {
-        self.bits
+/**
+Returns the probability `p_{k,t}` of an odd `k`-bit integer passing `t` rounds of MR testing with random bases
+is actually composite.
+
+Taken from FIPS-186.5[^FIPS], Section C.1, Eq. (2); FIPS in turn quotes Damgård et al[^Damgard].
+There it can be found as Eq. (4.1), with the factors bounded by Proposition 1 and Proposition 2 in Section 4.
+
+**Warning**:[^FIPS]
+Care must be taken when using the phrase "error probability" in connection with
+the recommended number of rounds of MR testing. The probability that a composite number of bit length `k`
+survives `t` rounds of Miller-Rabin testing is not the same as `p_{k,t}`, which is the probability that a
+number surviving `t` rounds of Miller-Rabin testing is composite. Ordinarily, the latter probability
+is the one that should be of most interest to a party responsible for generating primes, while the
+former may be more important to a party responsible for validating the primality of a number
+generated by someone else. However, for sufficiently large `k` (e.g., `k ≥ 51`), it can be shown that
+`p_{k,t} ≤ 4^{–t}` under the same assumptions concerning the selection of candidates as those made to
+obtain formula for `p_{k,t}` (see Damgård et al[^Damgard]).
+In such cases, `t = ceil(–log2(p_target)/2)` rounds of Miller-Rabin testing
+can be used to both generate and validate primes with `p_target` serving as an upper bound on both
+the probability that the generation process yields a composite number and the probability that a
+composite number would survive an attempt to validate its primality.
+
+[^FIPS]: FIPS-186.5 standard, <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
+[^Damgard]:
+    Damgård I, Landrock P, Pomerance C (1993) Average Case Error Estimates for the Strong
+    Probable Prime Test. Mathematics of Computation 61(203):177-194 (1993),
+    <https://www.ams.org/journals/mcom/1993-61-203/S0025-5718-1993-1189518-9/S0025-5718-1993-1189518-9.pdf>
+*/
+const fn pseudoprime_probability(k: u32, t: u32, cap_m: u32) -> f64 {
+    // Simplify Eq. (2) from [^FIPS] by factoring out `2^{k-2}`, which makes it
+    // p_{k,t} = 2.000743 ln(2) k 2^{-2} (
+    //     2^{-Mt} +
+    //     8 (pi^2 - 6) / 3 * sum(m=3..M) sum(j=2..m) 2^{m - (m-1)t - j - (k-1)/j})
+
+    // `2.00743 ln(2) k` comes from Proposition 2 in [^Damgard], from the estimate for `(pi(2^k) - pi(2^{k-1}))^{-1}`.
+
+    // Calculate the powers of 2 under the summations.
+    //
+    // Technically, only a few terms really contribute to the result, and the rest are extremely small in comparison.
+    // But finding out which ones specifically is a little messy.
+    // Can be done if more performance is desired.
+    let mut s = 0.;
+    let mut m = 3i32;
+    while m <= cap_m as i32 {
+        let mut j = 2i32;
+        while j <= m {
+            // Note that we are using `two_powf_upper_bound()` which means we are getting a slightly larger result,
+            // and therefore very slightly overestimating the resulting probability.
+            // This means safety is not compromised.
+            s += two_powf_upper_bound((m - (m - 1) * (t as i32) - j) as f64 - (k - 1) as f64 / j as f64);
+            j += 1;
+        }
+        m += 1;
+    }
+
+    const PI: f64 = core::f64::consts::PI;
+
+    // `2.00743 * ln(2) * 2^(-2)`
+    const COEFF: f64 = 0.3478611111678627;
+
+    COEFF * k as f64 * (1. / two_powi(cap_m * t) + 8. * (PI * PI - 6.) / 3. * s)
+}
+
+/// For a candidate of size `bit_length`, returns the minimum number of Miller-Rabin tests with random bases
+/// required for the probability of hitting a pseudoprime (that is, a composite for which all M-R tests pass)
+/// to be smaller than `2^{-log2_target}`.
+///
+/// Returns `None` if the number of iterations could not be found for the given bounds.
+///
+/// This function implements the formula prescribed by the FIPS.186-5 standard[^FIPS].
+///
+/// [^FIPS]: FIPS-186.5 standard, <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
+pub const fn minimum_mr_iterations(bit_length: u32, log2_target: u32) -> Option<usize> {
+    let mut t = 1;
+    while t <= (log2_target + 1) / 2 {
+        let cap_m_limit = floor_sqrt(4 * (bit_length - 1)) - 1;
+
+        let mut cap_m = 3;
+        while cap_m <= cap_m_limit {
+            let p = pseudoprime_probability(bit_length, t, cap_m);
+            if p < 1. / two_powi(log2_target) {
+                return Some(t as usize);
+            }
+            cap_m += 1;
+        }
+
+        t += 1;
     }
+
+    None
 }
 
 #[cfg(test)]
@@ -146,7 +235,7 @@ mod tests {
     #[cfg(feature = "tests-exhaustive")]
     use num_prime::nt_funcs::is_prime64;
 
-    use super::MillerRabin;
+    use super::{minimum_mr_iterations, MillerRabin};
     use crate::hazmat::{primes, pseudoprimes, random_odd_integer, SetBits, SmallPrimesSieve};
 
     #[test]
@@ -287,6 +376,29 @@ mod tests {
         test_large_primes(primes::PRIMES_1024);
     }
 
+    #[test]
+    fn mr_iterations() {
+        // The test values are taken from Table B.1 in FIPS-186.5 standard.
+
+        assert_eq!(minimum_mr_iterations(140, 100).unwrap(), 32);
+        assert_eq!(minimum_mr_iterations(140, 112).unwrap(), 38);
+        assert_eq!(minimum_mr_iterations(1024, 100).unwrap(), 4);
+        assert_eq!(minimum_mr_iterations(1024, 112).unwrap(), 5);
+
+        assert_eq!(minimum_mr_iterations(170, 100).unwrap(), 27);
+        assert_eq!(minimum_mr_iterations(170, 128).unwrap(), 41);
+        assert_eq!(minimum_mr_iterations(1536, 100).unwrap(), 3);
+        assert_eq!(minimum_mr_iterations(1536, 128).unwrap(), 4);
+
+        assert_eq!(minimum_mr_iterations(200, 100).unwrap(), 22);
+        assert_eq!(minimum_mr_iterations(200, 144).unwrap(), 44);
+        assert_eq!(minimum_mr_iterations(2048, 100).unwrap(), 2);
+        assert_eq!(minimum_mr_iterations(2048, 144).unwrap(), 4);
+
+        // Test the case where the requested error probability cannot be reached.
+        assert!(minimum_mr_iterations(128, 1024).is_none());
+    }
+
     #[cfg(feature = "tests-exhaustive")]
     #[test]
     fn exhaustive() {
diff --git a/src/hazmat/pseudoprimes.rs b/src/hazmat/pseudoprimes.rs
index 1503cc7..47cb3e4 100644
--- a/src/hazmat/pseudoprimes.rs
+++ b/src/hazmat/pseudoprimes.rs
@@ -8,7 +8,7 @@ use crypto_bigint::{U1536, U64};
 pub(crate) const EXHAUSTIVE_TEST_LIMIT: u32 = 500000;
 
 /// Extra strong Lucas pseudoprimes (OEIS:A217719) under `EXHAUSTIVE_TEST_LIMIT`.
-/// Should pass the strong Lucas test with Selfridge base (Baillie method A).
+/// Should pass the extra strong Lucas test with brute force base (Baillie method C).
 pub(crate) const EXTRA_STRONG_LUCAS: &[u32] = &[
     989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059, 72389, 73919, 75077, 100127, 113573, 125249, 137549,
     137801, 153931, 155819, 161027, 162133, 189419, 218321, 231703, 249331, 370229, 429479, 430127, 459191, 473891,
@@ -16,7 +16,7 @@ pub(crate) const EXTRA_STRONG_LUCAS: &[u32] = &[
 ];
 
 /// Strong Lucas pseudoprimes (OEIS:A217255) under `EXHAUSTIVE_TEST_LIMIT`.
-/// Should pass the extra strong Lucas test with brute force base (Baillie method C).
+/// Should pass the strong Lucas test with Selfridge base (Baillie method A).
 pub(crate) const STRONG_LUCAS: &[u32] = &[
     5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309, 58519, 75077, 97439, 100127, 113573, 115639, 130139,
     155819, 158399, 161027, 162133, 176399, 176471, 189419, 192509, 197801, 224369, 230691, 231703, 243629, 253259,
@@ -24,6 +24,7 @@ pub(crate) const STRONG_LUCAS: &[u32] = &[
 ];
 
 /// Almost extra strong Lucas pseudoprimes under `EXHAUSTIVE_TEST_LIMIT`.
+/// Should pass the almost extra strong Lucas test with brute force base (Baillie method C).
 /// Taken from D. Jacobsen, "Pseudoprime Statistics, Tables, and Data",
 /// <http://ntheory.org/pseudoprimes.html>
 pub(crate) const ALMOST_EXTRA_STRONG_LUCAS: &[u32] = &[
@@ -36,6 +37,23 @@ pub(crate) const ALMOST_EXTRA_STRONG_LUCAS: &[u32] = &[
 /// under `EXHAUSTIVE_TEST_LIMIT`.
 pub(crate) const LUCAS_V: &[u32] = &[913];
 
+/// Lucas pseudoprimes (OEIS:A217120) under `EXHAUSTIVE_TEST_LIMIT`.
+/// Should pass the regular Lucas test with Selfridge base (Baillie method A).
+/// Taken from D. Jacobsen, "Pseudoprime Statistics, Tables, and Data",
+/// <http://ntheory.org/pseudoprimes.html>
+pub(crate) const LUCAS: &[u32] = &[
+    323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179, 10877, 11419, 11663, 13919, 14839, 16109, 16211, 18407, 18971,
+    19043, 22499, 23407, 24569, 25199, 25877, 26069, 27323, 32759, 34943, 35207, 39059, 39203, 39689, 40309, 44099,
+    46979, 47879, 50183, 51983, 53663, 56279, 58519, 60377, 63881, 69509, 72389, 73919, 75077, 77219, 79547, 79799,
+    82983, 84419, 86063, 90287, 94667, 97019, 97439, 100127, 101919, 103739, 104663, 113573, 113849, 115439, 115639,
+    120581, 121103, 121393, 130139, 142883, 150079, 155819, 157079, 158399, 158717, 161027, 162133, 162719, 164699,
+    167969, 176399, 176471, 178949, 182513, 184781, 189419, 192509, 195227, 197801, 200147, 201871, 203699, 216659,
+    218129, 223901, 224369, 226529, 230159, 230691, 231703, 238999, 242079, 243629, 250277, 253259, 256409, 265481,
+    268349, 271991, 275099, 277399, 283373, 284171, 288919, 294527, 306287, 308699, 309959, 313499, 317249, 324899,
+    324911, 327359, 345913, 353219, 364229, 366799, 368351, 380393, 381923, 383921, 385307, 391169, 391859, 396899,
+    427349, 429263, 430127, 436409, 436589, 441599, 454607, 455519, 475799, 480689, 487199,
+];
+
 /// First 5 pseudoprimes for Lucas-V test[^Baillie2021].
 ///
 /// [^Baillie2021]: R. Baillie, A. Fiori, S. S. Wagstaff,
@@ -77,13 +95,6 @@ pub(crate) const BRUCKMAN_LUCAS: &[u32] = &[
     64681, 67861, 68251, 75077, 80189, 90061, 96049, 97921, 100065, 100127, 105281, 113573, 118441, 146611, 161027,
 ];
 
-/// Lucas pseudoprimes (OEIS:A217120).
-pub(crate) const LUCAS: &[u32] = &[
-    323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179, 10877, 11419, 11663, 13919, 14839, 16109, 16211, 18407, 18971,
-    19043, 22499, 23407, 24569, 25199, 25877, 26069, 27323, 32759, 34943, 35207, 39059, 39203, 39689, 40309, 44099,
-    46979, 47879,
-];
-
 /// A large Carmichael number.
 /// Source: F. Arnault, "Constructing Carmichael Numbers Which Are Strong
 /// Pseudoprimes to Several Bases". Journal of Symbolic Computation 20(2) 151–161 (1995),
diff --git a/src/lib.rs b/src/lib.rs
index 19ff3ec..f763e85 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -22,7 +22,10 @@ mod presets;
 mod traits;
 
 pub use generic::{sieve_and_find, SieveIterator};
-pub use presets::{generate_prime_with_rng, generate_safe_prime_with_rng, is_prime_with_rng, is_safe_prime_with_rng};
+pub use presets::{
+    fips_is_prime_with_rng, fips_is_safe_prime_with_rng, generate_prime_with_rng, generate_safe_prime_with_rng,
+    is_prime_with_rng, is_safe_prime_with_rng,
+};
 pub use traits::{RandomPrimeWithRng, SieveFactory};
 
 #[cfg(feature = "default-rng")]
diff --git a/src/presets.rs b/src/presets.rs
index c08c98f..eb1f96b 100644
--- a/src/presets.rs
+++ b/src/presets.rs
@@ -1,4 +1,4 @@
-use crypto_bigint::{Integer, Limb, Odd, RandomBits, RandomMod};
+use crypto_bigint::{Integer, Odd, RandomBits, RandomMod, Word};
 use rand_core::CryptoRng;
 
 #[cfg(feature = "default-rng")]
@@ -6,7 +6,9 @@ use rand_core::{OsRng, TryRngCore};
 
 use crate::{
     generic::sieve_and_find,
-    hazmat::{lucas_test, AStarBase, LucasCheck, MillerRabin, Primality, SetBits, SmallPrimesSieveFactory},
+    hazmat::{
+        lucas_test, AStarBase, LucasCheck, MillerRabin, Primality, SelfridgeBase, SetBits, SmallPrimesSieveFactory,
+    },
 };
 
 #[cfg(feature = "multicore")]
@@ -156,6 +158,10 @@ where
     .expect("will produce a result eventually")
 }
 
+fn equals_primitive<T: Integer>(num: &T, primitive: Word) -> bool {
+    num.bits_vartime() <= u16::BITS && num.as_ref()[0].0 == primitive
+}
+
 /// Probabilistically checks if the given number is prime using the provided RNG.
 ///
 /// Performed checks:
@@ -181,83 +187,169 @@ where
 ///       "Strengthening the Baillie-PSW primality test",
 ///       Math. Comp. 90 1931-1955 (2021),
 ///       DOI: [10.1090/mcom/3616](https://doi.org/10.1090/mcom/3616)
-pub fn is_prime_with_rng<T: Integer + RandomMod, R: CryptoRng + ?Sized>(rng: &mut R, num: &T) -> bool {
-    if num == &T::from_limb_like(Limb::from(2u32), num) {
+pub fn is_prime_with_rng<T: Integer + RandomMod>(rng: &mut (impl CryptoRng + ?Sized), candidate: &T) -> bool {
+    if equals_primitive(candidate, 1) {
+        return false;
+    }
+
+    if equals_primitive(candidate, 2) {
         return true;
     }
 
-    match Odd::new(num.clone()).into() {
-        Some(x) => _is_prime_with_rng(rng, x),
-        None => false,
+    let odd_candidate: Odd<T> = match Odd::new(candidate.clone()).into() {
+        Some(x) => x,
+        None => return false,
+    };
+
+    let mr = MillerRabin::new(odd_candidate.clone());
+
+    if !mr.test_base_two().is_probably_prime() {
+        return false;
     }
+
+    match lucas_test(odd_candidate, AStarBase, LucasCheck::Strong) {
+        Primality::Composite => return false,
+        Primality::Prime => return true,
+        _ => {}
+    }
+
+    // The random base test only makes sense when `num > 3`.
+    if !equals_primitive(candidate, 3) && !mr.test_random_base(rng).is_probably_prime() {
+        return false;
+    }
+
+    true
 }
 
 /// Probabilistically checks if the given number is a safe prime using the provided RNG.
 ///
 /// See [`is_prime_with_rng`] for details about the performed checks.
-pub fn is_safe_prime_with_rng<T: Integer + RandomMod, R: CryptoRng + ?Sized>(rng: &mut R, num: &T) -> bool {
-    // Since, by the definition of safe prime, `(num - 1) / 2` must also be prime,
-    // and therefore odd, `num` has to be equal to 3 modulo 4.
+pub fn is_safe_prime_with_rng<T: Integer + RandomMod>(rng: &mut (impl CryptoRng + ?Sized), candidate: &T) -> bool {
+    // Since, by the definition of safe prime, `(candidate - 1) / 2` must also be prime,
+    // and therefore odd, `candidate` has to be equal to 3 modulo 4.
     // 5 is the only exception, so we check for it.
-    if num == &T::from_limb_like(Limb::from(5u32), num) {
+    if equals_primitive(candidate, 5) {
         return true;
     }
 
     // Safe primes are always of the form 4k + 3 (i.e. n ≡ 3 mod 4)
     // The last two digits of a binary number give you its value modulo 4.
-    // Primes p=4n+3 will always end in 11​ in binary because p ≡ 3 mod 4.
-    if num.as_ref()[0].0 & 3 != 3 {
+    // Primes p=4n+3 will always end in 11 in binary because p ≡ 3 mod 4.
+    if candidate.as_ref()[0].0 & 3 != 3 {
         return false;
     }
 
-    // These are ensured to be odd by the check above.
-    let odd_num = Odd::new(num.clone()).expect("`num` is odd here given the checks above");
-    let odd_half_num = Odd::new(num.wrapping_shr_vartime(1)).expect("The binary rep of `num` ends in `11`, so shifting right by one is guaranteed leave a `1` at the end, so it's odd");
-
-    _is_prime_with_rng(rng, odd_num) && _is_prime_with_rng(rng, odd_half_num)
+    is_prime_with_rng(rng, candidate) && is_prime_with_rng(rng, &candidate.wrapping_shr_vartime(1))
 }
 
-/// Checks for primality.
-/// First run a Miller-Rabin test with base 2
-/// If the outcome of M-R is "probably prime", then run a Lucas test
-/// If the Lucas test is inconclusive, run a Miller-Rabin with random base and unless this second
-/// M-R test finds it's composite, then conclude that it's prime.
-fn _is_prime_with_rng<T: Integer + RandomMod, R: CryptoRng + ?Sized>(rng: &mut R, candidate: Odd<T>) -> bool {
-    let mr = MillerRabin::new(candidate.clone());
-
-    if !mr.test_base_two().is_probably_prime() {
+/// Probabilistically checks if the given number is prime using the provided RNG
+/// according to FIPS-186.5[^FIPS] standard.
+///
+/// Performed checks:
+/// - `mr_iterations` of Miller-Rabin check with random bases;
+/// - Regular Lucas check with Selfridge base (see [`SelfridgeBase`] for details), if `add_lucas_test` is `true`.
+///
+/// See [`MillerRabin`] and [`lucas_test`] for more details about the checks;
+/// use [`minimum_mr_iterations`](`crate::hazmat::minimum_mr_iterations`)
+/// to calculate the number of required iterations.
+///
+/// [^FIPS]: FIPS-186.5 standard, <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
+pub fn fips_is_prime_with_rng<T: Integer + RandomMod>(
+    rng: &mut (impl CryptoRng + ?Sized),
+    candidate: &T,
+    mr_iterations: usize,
+    add_lucas_test: bool,
+) -> bool {
+    if equals_primitive(candidate, 1) {
         return false;
     }
 
-    match lucas_test(candidate, AStarBase, LucasCheck::Strong) {
-        Primality::Composite => return false,
-        Primality::Prime => return true,
-        _ => {}
+    if equals_primitive(candidate, 2) {
+        return true;
     }
 
-    // The random base test only makes sense when `num > 3`.
-    if mr.bits() > 2 && !mr.test_random_base(rng).is_probably_prime() {
-        return false;
+    let odd_candidate: Odd<T> = match Odd::new(candidate.clone()).into() {
+        Some(x) => x,
+        None => return false,
+    };
+
+    // The random base test only makes sense when `candidate > 3`.
+    if !equals_primitive(candidate, 3) {
+        let mr = MillerRabin::new(odd_candidate.clone());
+        for _ in 0..mr_iterations {
+            if !mr.test_random_base(rng).is_probably_prime() {
+                return false;
+            }
+        }
+    }
+
+    if add_lucas_test {
+        match lucas_test(odd_candidate, SelfridgeBase, LucasCheck::Strong) {
+            Primality::Composite => return false,
+            Primality::Prime => return true,
+            _ => {}
+        }
     }
 
     true
 }
 
+/// Probabilistically checks if the given number is a safe prime using the provided RNG
+/// according to FIPS-186.5[^FIPS] standard.
+///
+/// See [`fips_is_prime_with_rng`] for details about the performed checks.
+///
+/// [^FIPS]: FIPS-186.5 standard, <https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-5.pdf>
+pub fn fips_is_safe_prime_with_rng<T: Integer + RandomMod>(
+    rng: &mut (impl CryptoRng + ?Sized),
+    candidate: &T,
+    mr_iterations: usize,
+    add_lucas_test: bool,
+) -> bool {
+    // Since, by the definition of safe prime, `(candidate - 1) / 2` must also be prime,
+    // and therefore odd, `candidate` has to be equal to 3 modulo 4.
+    // 5 is the only exception, so we check for it.
+    if equals_primitive(candidate, 5) {
+        return true;
+    }
+
+    // Safe primes are always of the form 4k + 3 (i.e. n ≡ 3 mod 4)
+    // The last two digits of a binary number give you its value modulo 4.
+    // Primes p=4n+3 will always end in 11 in binary because p ≡ 3 mod 4.
+    if candidate.as_ref()[0].0 & 3 != 3 {
+        return false;
+    }
+
+    fips_is_prime_with_rng(rng, candidate, mr_iterations, add_lucas_test)
+        && fips_is_prime_with_rng(rng, &candidate.wrapping_shr_vartime(1), mr_iterations, add_lucas_test)
+}
+
 #[cfg(test)]
 mod tests {
-    use crypto_bigint::{nlimbs, BoxedUint, CheckedAdd, Uint, Word, U128, U64};
+    use crypto_bigint::{nlimbs, BoxedUint, CheckedAdd, Integer, RandomMod, Uint, Word, U128, U64};
     use num_prime::nt_funcs::is_prime64;
     use rand_core::{OsRng, TryRngCore};
 
     use super::{
-        generate_prime, generate_prime_with_rng, generate_safe_prime, generate_safe_prime_with_rng, is_prime,
-        is_safe_prime,
+        fips_is_prime_with_rng, fips_is_safe_prime_with_rng, generate_prime, generate_prime_with_rng,
+        generate_safe_prime, generate_safe_prime_with_rng, is_prime, is_safe_prime,
     };
-    use crate::hazmat::{primes, pseudoprimes};
+    use crate::hazmat::{minimum_mr_iterations, primes, pseudoprimes};
+
+    fn fips_is_prime<T: Integer + RandomMod>(num: &T) -> bool {
+        let mr_iterations = minimum_mr_iterations(128, 100).unwrap();
+        fips_is_prime_with_rng(&mut OsRng.unwrap_err(), num, mr_iterations, false)
+    }
+
+    fn fips_is_safe_prime<T: Integer + RandomMod>(num: &T) -> bool {
+        let mr_iterations = minimum_mr_iterations(128, 100).unwrap();
+        fips_is_safe_prime_with_rng(&mut OsRng.unwrap_err(), num, mr_iterations, false)
+    }
 
     fn test_large_primes<const L: usize>(nums: &[Uint<L>]) {
         for num in nums {
             assert!(is_prime(num));
+            assert!(fips_is_prime(num));
         }
     }
 
@@ -273,6 +365,7 @@ mod tests {
     fn test_pseudoprimes(nums: &[u32]) {
         for num in nums {
             assert!(!is_prime(&U64::from(*num)));
+            assert!(!fips_is_prime(&U64::from(*num)));
         }
     }
 
@@ -288,19 +381,23 @@ mod tests {
 
         for num in pseudoprimes::STRONG_FIBONACCI {
             assert!(!is_prime(num));
+            assert!(!fips_is_prime(num));
         }
 
         assert!(!is_prime(&pseudoprimes::LARGE_CARMICHAEL_NUMBER));
+        assert!(!fips_is_prime(&pseudoprimes::LARGE_CARMICHAEL_NUMBER));
     }
 
     fn test_cunningham_chain<const L: usize>(length: usize, num: &Uint<L>) {
         let mut next = *num;
         for i in 0..length {
             assert!(is_prime(&next));
+            assert!(fips_is_prime(&next));
 
             // The start of the chain isn't a safe prime by definition
             if i > 0 {
                 assert!(is_safe_prime(&next));
+                assert!(fips_is_safe_prime(&next));
             }
 
             next = next.wrapping_shl_vartime(1).checked_add(&Uint::<L>::ONE).unwrap();
@@ -308,6 +405,7 @@ mod tests {
 
         // The chain ended.
         assert!(!is_prime(&next));
+        assert!(!fips_is_prime(&next));
     }
 
     #[test]
@@ -326,6 +424,7 @@ mod tests {
             let p: U128 = generate_prime(bit_length);
             assert!(p.bits_vartime() == bit_length);
             assert!(is_prime(&p));
+            assert!(fips_is_prime(&p));
         }
     }
 
@@ -336,6 +435,7 @@ mod tests {
             assert!(p.bits_vartime() == bit_length);
             assert!(p.to_words().len() == nlimbs!(bit_length));
             assert!(is_prime(&p));
+            assert!(fips_is_prime(&p));
         }
     }
 
@@ -345,6 +445,7 @@ mod tests {
             let p: U128 = generate_safe_prime(bit_length);
             assert!(p.bits_vartime() == bit_length);
             assert!(is_safe_prime(&p));
+            assert!(fips_is_safe_prime(&p));
         }
     }
 
@@ -355,6 +456,7 @@ mod tests {
             assert!(p.bits_vartime() == bit_length);
             assert!(p.to_words().len() == nlimbs!(bit_length));
             assert!(is_safe_prime(&p));
+            assert!(fips_is_safe_prime(&p));
         }
     }
 
@@ -365,11 +467,30 @@ mod tests {
             let is_safe_prime_ref = is_prime_ref && is_prime64(num / 2);
 
             let num_uint = U64::from(num);
+
             let is_prime_test = is_prime(&num_uint);
+            assert_eq!(
+                is_prime_ref, is_prime_test,
+                "num={num}, expected={is_prime_ref}, actual={is_prime_test}"
+            );
+
             let is_safe_prime_test = is_safe_prime(&num_uint);
+            assert_eq!(
+                is_safe_prime_ref, is_safe_prime_test,
+                "num={num}, expected={is_safe_prime_ref}, actual={is_safe_prime_test}"
+            );
 
-            assert_eq!(is_prime_ref, is_prime_test);
-            assert_eq!(is_safe_prime_ref, is_safe_prime_test);
+            let is_prime_test = fips_is_prime(&num_uint);
+            assert_eq!(
+                is_prime_ref, is_prime_test,
+                "num={num}, expected={is_prime_ref}, actual={is_prime_test}"
+            );
+
+            let is_safe_prime_test = fips_is_safe_prime(&num_uint);
+            assert_eq!(
+                is_safe_prime_ref, is_safe_prime_test,
+                "num={num}, expected={is_safe_prime_ref}, actual={is_safe_prime_test}"
+            );
         }
     }
 
@@ -380,6 +501,7 @@ mod tests {
         // so the initial sieving cannot tell if it is prime or not,
         // and a full primality test is run.
         assert!(!is_safe_prime(&U64::from(17881u32 * 17891u32)));
+        assert!(!fips_is_safe_prime(&U64::from(17881u32 * 17891u32)));
     }
 
     #[test]
@@ -475,8 +597,8 @@ mod tests_openssl {
     use openssl::bn::{BigNum, BigNumContext};
     use rand_core::{OsRng, TryRngCore};
 
-    use super::{generate_prime, is_prime};
-    use crate::hazmat::{random_odd_integer, SetBits};
+    use super::{fips_is_prime_with_rng, generate_prime, is_prime};
+    use crate::hazmat::{minimum_mr_iterations, random_odd_integer, SetBits};
 
     fn openssl_is_prime(num: &BigNum, ctx: &mut BigNumContext) -> bool {
         num.is_prime(64, ctx).unwrap()
@@ -501,21 +623,33 @@ mod tests_openssl {
             assert!(openssl_is_prime(&p_bn, &mut ctx), "OpenSSL reports {p} as composite",);
         }
 
+        let mr_iterations = minimum_mr_iterations(U128::BITS, 100).unwrap();
+
         // Generate primes with OpenSSL, check them
         let mut p_bn = BigNum::new().unwrap();
         for _ in 0..100 {
             p_bn.generate_prime(128, false, None, None).unwrap();
             let p = from_openssl(&p_bn);
             assert!(is_prime(&p), "we report {p} as composite");
+            assert!(
+                fips_is_prime_with_rng(&mut OsRng.unwrap_err(), &p, mr_iterations, false),
+                "we report {p} as composite"
+            );
         }
 
         // Generate random numbers, check if our test agrees with OpenSSL
         for _ in 0..100 {
-            let p = random_odd_integer::<U128, _>(&mut OsRng.unwrap_err(), NonZero::new(128).unwrap(), SetBits::Msb)
-                .unwrap();
-            let actual = is_prime(p.as_ref());
+            let p = random_odd_integer::<U128, _>(&mut OsRng, NonZero::new(128).unwrap(), SetBits::Msb).unwrap();
             let p_bn = to_openssl(&p);
             let expected = openssl_is_prime(&p_bn, &mut ctx);
+
+            let actual = is_prime(p.as_ref());
+            assert_eq!(
+                actual, expected,
+                "difference between OpenSSL and us: OpenSSL reports {expected}, we report {actual}",
+            );
+
+            let actual = fips_is_prime_with_rng(&mut OsRng.unwrap_err(), p.as_ref(), mr_iterations, false);
             assert_eq!(
                 actual, expected,
                 "difference between OpenSSL and us: OpenSSL reports {expected}, we report {actual}",
@@ -536,8 +670,8 @@ mod tests_gmp {
         Integer,
     };
 
-    use super::{generate_prime, is_prime};
-    use crate::hazmat::{random_odd_integer, SetBits};
+    use super::{fips_is_prime_with_rng, generate_prime, is_prime};
+    use crate::hazmat::{minimum_mr_iterations, random_odd_integer, SetBits};
 
     fn gmp_is_prime(num: &Integer) -> bool {
         matches!(num.is_probably_prime(25), IsPrime::Yes | IsPrime::Probably)
@@ -560,6 +694,8 @@ mod tests_gmp {
             assert!(gmp_is_prime(&p_bn), "GMP reports {p} as composite");
         }
 
+        let mr_iterations = minimum_mr_iterations(U128::BITS, 100).unwrap();
+
         // Generate primes with GMP, check them
         for _ in 0..100 {
             let start =
@@ -569,15 +705,26 @@ mod tests_gmp {
             let p_bn = start_bn.next_prime();
             let p = from_gmp(&p_bn);
             assert!(is_prime(&p), "we report {p} as composite");
+            assert!(
+                fips_is_prime_with_rng(&mut OsRng.unwrap_err(), &p, mr_iterations, false),
+                "we report {p} as composite"
+            );
         }
 
         // Generate random numbers, check if our test agrees with GMP
         for _ in 0..100 {
             let p = random_odd_integer::<U128, _>(&mut OsRng.unwrap_err(), NonZero::new(128).unwrap(), SetBits::Msb)
                 .unwrap();
-            let actual = is_prime(p.as_ref());
             let p_bn = to_gmp(&p);
             let expected = gmp_is_prime(&p_bn);
+
+            let actual = is_prime(p.as_ref());
+            assert_eq!(
+                actual, expected,
+                "difference between GMP and us: GMP reports {expected}, we report {actual}",
+            );
+
+            let actual = fips_is_prime_with_rng(&mut OsRng.unwrap_err(), p.as_ref(), mr_iterations, false);
             assert_eq!(
                 actual, expected,
                 "difference between GMP and us: GMP reports {expected}, we report {actual}",