crypto_primes/hazmat/
miller_rabin.rs

1//! Miller-Rabin primality test.
2
3use crypto_bigint::{Integer, Limb, Monty, NonZero as CTNonZero, Odd, PowBoundedExp, RandomMod, Square};
4use rand_core::CryptoRngCore;
5
6use super::Primality;
7
8/// Precomputed data used to perform Miller-Rabin primality test[^Pomerance1980].
9///
10/// The numbers that pass it are commonly called "strong probable primes"
11/// (or "strong pseudoprimes" if they are, in fact, composite).
12///
13/// [^Pomerance1980]:
14///   C. Pomerance, J. L. Selfridge, S. S. Wagstaff "The Pseudoprimes to 25*10^9",
15///   Math. Comp. 35 1003-1026 (1980),
16///   DOI: [10.2307/2006210](https://dx.doi.org/10.2307/2006210)
17#[derive(Clone, Debug, PartialEq, Eq)]
18pub struct MillerRabin<T: Integer> {
19    // The odd number that may or may not be a prime.
20    candidate: T,
21    /// The number of bits necessesary to represent the candidate. Note: this is not the number of
22    /// bits used by a `T` in memory.
23    bits: u32,
24    /// Pre-computed parameters for the Montgomery form of `T`.
25    montgomery_params: <<T as Integer>::Monty as Monty>::Params,
26    /// The number 1 in Montgomery form.
27    one: <T as Integer>::Monty,
28    /// The number -1 in Montgomery form.
29    minus_one: <T as Integer>::Monty,
30    /// The `s` exponent in the Miller-Rabin test, that finds `s` and `d` odd s.t. `candidate - 1 ==
31    /// 2^s * d` (the pair `s` and `d` is unique).
32    s: u32,
33    /// The `d` factor in the Miller-Rabin test, that finds `s` and `d` odd s.t. `candidate -
34    /// 1 == 2^s * d` (the pair `s` and `d` is unique).
35    d: T,
36}
37
38impl<T: Integer + RandomMod> MillerRabin<T> {
39    /// Initializes a Miller-Rabin test for `candidate`.
40    pub fn new(candidate: Odd<T>) -> Self {
41        let params = <T as Integer>::Monty::new_params_vartime(candidate.clone());
42        let m_one = <T as Integer>::Monty::one(params.clone());
43        let m_minus_one = -m_one.clone();
44
45        let one = T::one_like(candidate.as_ref());
46
47        // Find `s` and odd `d` such that `candidate - 1 == 2^s * d`.
48        let (s, d) = if candidate.as_ref() == &one {
49            (0, one)
50        } else {
51            let candidate_minus_one = candidate.wrapping_sub(&one);
52            let s = candidate_minus_one.trailing_zeros_vartime();
53            // Will not overflow because `candidate` is odd and greater than 1.
54            let d = candidate_minus_one
55                .overflowing_shr_vartime(s)
56                .expect("shifting by `s` is within range by construction: `candidate` is odd and greater than 1");
57            (s, d)
58        };
59
60        Self {
61            bits: candidate.bits_vartime(),
62            candidate: candidate.get(),
63            montgomery_params: params,
64            one: m_one,
65            minus_one: m_minus_one,
66            s,
67            d,
68        }
69    }
70
71    /// Perform a Miller-Rabin check with a given base.
72    pub fn test(&self, base: &T) -> Primality {
73        // One could check here if `gcd(base, candidate) == 1` and return `Composite` otherwise.
74        // In practice it doesn't make any performance difference in normal operation.
75
76        let base = <T as Integer>::Monty::new(base.clone(), self.montgomery_params.clone());
77
78        // Implementation detail: bounded exp gets faster every time we decrease the bound
79        // by the window length it uses, which is currently 4 bits.
80        // So even when the bound isn't low enough that the number can fit
81        // in a smaller number of limbs, there is still a performance gain
82        // from specifying the bound.
83        let mut test = base.pow_bounded_exp(&self.d, self.bits);
84
85        if test == self.one || test == self.minus_one {
86            return Primality::ProbablyPrime;
87        }
88        for _ in 1..self.s {
89            test = test.square();
90            if test == self.one {
91                return Primality::Composite;
92            } else if test == self.minus_one {
93                return Primality::ProbablyPrime;
94            }
95        }
96        Primality::Composite
97    }
98
99    /// Perform a Miller-Rabin check with base 2.
100    pub fn test_base_two(&self) -> Primality {
101        self.test(&T::from_limb_like(Limb::from(2u32), &self.candidate))
102    }
103
104    /// Perform a Miller-Rabin check with a random base (in the range `[3, candidate-2]`)
105    /// drawn using the provided RNG.
106    ///
107    /// Note: panics if `candidate == 3` (so the range above is empty).
108    pub fn test_random_base(&self, rng: &mut (impl CryptoRngCore + ?Sized)) -> Primality {
109        // We sample a random base from the range `[3, candidate-2]`:
110        // - we have a separate method for base 2;
111        // - the test holds trivially for bases 1 or `candidate-1`.
112        if self.candidate.bits() < 3 {
113            panic!("No suitable random base possible when `candidate == 3`; use the base 2 test.")
114        }
115
116        let range = self.candidate.wrapping_sub(&T::from(4u32));
117        // Can unwrap here since `candidate` is odd, and `candidate >= 4` (as checked above)
118        let range_nonzero = CTNonZero::new(range).expect("the range should be non-zero by construction");
119        // This should not overflow as long as `random_mod()` behaves according to the contract
120        // (that is, returns a number within the given range).
121        let random = T::random_mod(rng, &range_nonzero)
122            .checked_add(&T::from(3u32))
123            .expect("addition should not overflow by construction");
124        self.test(&random)
125    }
126
127    /// Returns the number of bits necessary to represent the candidate.
128    /// NOTE: This is different than the number of bits of *storage* the integer takes up.
129    ///
130    /// For example, a U512 type occupies 8 64-bit words, but the number `7` contained in such a type
131    /// has a bit length of 3 because 7 is `b111`.
132    pub(crate) fn bits(&self) -> u32 {
133        self.bits
134    }
135}
136
137#[cfg(test)]
138mod tests {
139    use alloc::format;
140    use core::num::NonZero;
141
142    use crypto_bigint::{Integer, Odd, RandomMod, Uint, U1024, U128, U1536, U64};
143    use rand_chacha::ChaCha8Rng;
144    use rand_core::{CryptoRngCore, OsRng, SeedableRng};
145
146    #[cfg(feature = "tests-exhaustive")]
147    use num_prime::nt_funcs::is_prime64;
148
149    use super::MillerRabin;
150    use crate::hazmat::{primes, pseudoprimes, random_odd_integer, SetBits, SmallPrimesSieve};
151
152    #[test]
153    fn miller_rabin_derived_traits() {
154        let mr = MillerRabin::new(Odd::new(U64::ONE).unwrap());
155        assert!(format!("{mr:?}").starts_with("MillerRabin"));
156        assert_eq!(mr.clone(), mr);
157    }
158
159    #[test]
160    #[should_panic(expected = "No suitable random base possible when `candidate == 3`; use the base 2 test.")]
161    fn random_base_range_check() {
162        let mr = MillerRabin::new(Odd::new(U64::from(3u32)).unwrap());
163        mr.test_random_base(&mut OsRng);
164    }
165
166    fn is_spsp(num: u32) -> bool {
167        pseudoprimes::STRONG_BASE_2.iter().any(|x| *x == num)
168    }
169
170    fn random_checks<T: Integer + RandomMod>(
171        rng: &mut (impl CryptoRngCore + ?Sized),
172        mr: &MillerRabin<T>,
173        count: usize,
174    ) -> usize {
175        (0..count)
176            .map(|_| -> usize { mr.test_random_base(rng).is_probably_prime().into() })
177            .sum()
178    }
179
180    fn test_composites(numbers: &[u32], expected_result: bool) {
181        let mut rng = ChaCha8Rng::from_seed(*b"01234567890123456789012345678901");
182        for num in numbers.iter() {
183            let base_2_false_positive = is_spsp(*num);
184            let actual_expected_result = if base_2_false_positive { true } else { expected_result };
185
186            // A random base MR test is expected to report a composite as a prime
187            // with about 1/4 probability. So we're expecting less than
188            // 35 out of 100 false positives, seems to work.
189
190            let mr = MillerRabin::new(Odd::new(U64::from(*num)).unwrap());
191            assert_eq!(mr.test_base_two().is_probably_prime(), actual_expected_result);
192            let reported_prime = random_checks(&mut rng, &mr, 100);
193            assert!(
194                reported_prime < 35,
195                "reported as prime in {reported_prime} out of 100 tests",
196            );
197        }
198    }
199
200    #[test]
201    fn trivial() {
202        let mut rng = ChaCha8Rng::from_seed(*b"01234567890123456789012345678901");
203        let start = random_odd_integer::<U1024>(&mut rng, NonZero::new(1024).unwrap(), SetBits::Msb).unwrap();
204        for num in SmallPrimesSieve::new(start.get(), NonZero::new(1024).unwrap(), false).take(10) {
205            let mr = MillerRabin::new(Odd::new(num).unwrap());
206
207            // Trivial tests, must always be true.
208            assert!(mr.test(&1u32.into()).is_probably_prime());
209            assert!(mr.test(&num.wrapping_sub(&1u32.into())).is_probably_prime());
210        }
211    }
212
213    #[test]
214    fn mersenne_prime() {
215        let mut rng = ChaCha8Rng::from_seed(*b"01234567890123456789012345678901");
216
217        // Mersenne prime 2^127-1
218        let num = Odd::new(U128::from_be_hex("7fffffffffffffffffffffffffffffff")).unwrap();
219
220        let mr = MillerRabin::new(num);
221        assert!(mr.test_base_two().is_probably_prime());
222        for _ in 0..10 {
223            assert!(mr.test_random_base(&mut rng).is_probably_prime());
224        }
225    }
226
227    #[test]
228    fn strong_fibonacci_pseudoprimes() {
229        let mut rng = ChaCha8Rng::from_seed(*b"01234567890123456789012345678901");
230
231        for num in pseudoprimes::STRONG_FIBONACCI.iter() {
232            let mr = MillerRabin::new(Odd::new(*num).unwrap());
233            assert!(!mr.test_base_two().is_probably_prime());
234            for _ in 0..1000 {
235                assert!(!mr.test_random_base(&mut rng).is_probably_prime());
236            }
237        }
238    }
239
240    #[test]
241    fn strong_pseudoprimes_base_2() {
242        // These known exceptions for the base 2 MR test.
243        test_composites(pseudoprimes::STRONG_BASE_2, true);
244    }
245
246    #[test]
247    fn lucas_pseudoprimes() {
248        // Cross-test against the pseudoprimes that circumvent the Lucas test.
249        // We expect the MR test to correctly classify them as composites (most of the time).
250        test_composites(pseudoprimes::EXTRA_STRONG_LUCAS, false);
251        test_composites(pseudoprimes::STRONG_LUCAS, false);
252        test_composites(pseudoprimes::ALMOST_EXTRA_STRONG_LUCAS, false);
253        test_composites(pseudoprimes::FIBONACCI, false);
254        test_composites(pseudoprimes::BRUCKMAN_LUCAS, false);
255        test_composites(pseudoprimes::LUCAS, false);
256    }
257
258    #[test]
259    fn large_carmichael_number() {
260        let mr = MillerRabin::new(Odd::new(pseudoprimes::LARGE_CARMICHAEL_NUMBER).unwrap());
261
262        // It is known to pass MR tests for all prime bases <307
263        assert!(mr.test_base_two().is_probably_prime());
264        assert!(mr.test(&U1536::from(293u64)).is_probably_prime());
265
266        // A test with base 307 correctly reports the number as composite.
267        assert!(!mr.test(&U1536::from(307u64)).is_probably_prime());
268    }
269
270    fn test_large_primes<const L: usize>(nums: &[Uint<L>]) {
271        let mut rng = ChaCha8Rng::from_seed(*b"01234567890123456789012345678901");
272        for num in nums {
273            let mr = MillerRabin::new(Odd::new(*num).unwrap());
274            assert!(mr.test_base_two().is_probably_prime());
275            for _ in 0..10 {
276                assert!(mr.test_random_base(&mut rng).is_probably_prime());
277            }
278        }
279    }
280
281    #[test]
282    fn large_primes() {
283        test_large_primes(primes::PRIMES_128);
284        test_large_primes(primes::PRIMES_256);
285        test_large_primes(primes::PRIMES_384);
286        test_large_primes(primes::PRIMES_512);
287        test_large_primes(primes::PRIMES_1024);
288    }
289
290    #[cfg(feature = "tests-exhaustive")]
291    #[test]
292    fn exhaustive() {
293        // Test all the odd numbers up to the limit where we know the false positives,
294        // and compare the results with the reference.
295        for num in (3..pseudoprimes::EXHAUSTIVE_TEST_LIMIT).step_by(2) {
296            let res_ref = is_prime64(num.into());
297
298            let spsp = is_spsp(num);
299
300            let mr = MillerRabin::new(Odd::new(U64::from(num)).unwrap());
301            let res = mr.test_base_two().is_probably_prime();
302            let expected = spsp || res_ref;
303            assert_eq!(
304                res, expected,
305                "Miller-Rabin: n={num}, expected={expected}, actual={res}",
306            );
307        }
308    }
309}
crypto_primes/hazmat/miller_rabin.rs

crypto_primes/hazmat/
miller_rabin.rs
