primality_test/
lib.rs

1//! Provide a method to determine whether an unsigned integer is a prime number through the `IsPrime` Trait.
2//!
3//! In the current implementation, this crate uses the Miller-Rabin primality test.  
4//! The Miller-Rabin primality test is known to have witnesses that can conclusively determine unsigned integers of at most 64 bits.  
5//! In this crate, the following information is used to select the witnesses.  
6//!
7//! [Deterministic variants of the Miller-Rabin primality test](https://miller-rabin.appspot.com/)
8
9mod constfn;
10mod montgomery;
11mod sieve;
12
13pub use constfn::is_prime;
14use montgomery::Montgomery;
15pub use sieve::LinearSieve;
16
17const SMALL_PRIMES_MEMO: LinearSieve<255> = LinearSieve::new();
18
19pub trait IsPrime {
20    fn is_prime(&self) -> bool;
21}
22
23macro_rules! impl_is_prime {
24    ( $t:ty, $witness_ty:ty, [ $( $witness:expr ),* ] ) => {
25        impl IsPrime for $t {
26            /// Determine an unsigned integer is prime number or not.
27            ///
28            /// # Examples
29            /// ```rust
30            /// use primality_test::IsPrime;
31            ///
32            /// assert!(998244353u32.is_prime());
33            /// assert!(!561u16.is_prime());
34            ///
35            /// let primes = (1..20u16).filter(IsPrime::is_prime).collect::<Vec<_>>();
36            /// assert_eq!(primes, vec![2, 3, 5, 7, 11, 13, 17, 19]);
37            /// ```
38            fn is_prime(&self) -> bool {
39                let p = *self as $witness_ty;
40
41                if p == 1 || p & 1 == 0 {
42                    return p == 2;
43                }
44
45                if <$witness_ty>::BITS == 64 && (p as u64) < 1795265022 {
46                    return (p as u32).is_prime();
47                }
48
49                if p < SMALL_PRIMES_MEMO.len() as $witness_ty {
50                    return SMALL_PRIMES_MEMO.is_prime(p as usize);
51                }
52
53                let mont = Montgomery::<$witness_ty>::new(p);
54
55                let s = (p - 1).trailing_zeros();
56                let t = (p - 1) >> s;
57
58                [$( $witness ),*]
59                    .iter()
60                    // These two lines are necessary if any of the witnesses may be greater than or equal to `p`.
61                    // The maximum witness for upto 32-bit integers is 61, and when `p` <= 61, this line is unreachable because they dictionary SMALL_PRIMES_MEMO.
62                    // The maximum witness for 64-bit integers is 1795265022, but since numbers less than this can be expressed in 32 bits, they are cast to u32 and judged again, so this point is unreachable.
63                    // .map(|&a| a % p)
64                    // .filter(|&a| a != 0)
65                    .all(|&a| {
66                        let a = mont.convert(a);
67                        let at = mont.pow(a, t);
68                        // a^t = 1 (mod p) or a^t = -1 (mod p)
69                        if at == mont.r || at == p - mont.r {
70                            return true;
71                        }
72
73                        // found i satisfying a^((2^i)*t) = -1 (mod p)
74                        (1..s)
75                            .scan(at, |at, _| {
76                                *at = mont.multiply(*at, *at);
77                                Some(*at)
78                            })
79                            .any(|at| at == p - mont.r)
80                    })
81            }
82        }
83    };
84}
85
86impl_is_prime!(u8, u8, [2, 7, 61]);
87impl_is_prime!(u16, u16, [2, 7, 61]);
88impl_is_prime!(u32, u32, [2, 7, 61]);
89impl_is_prime!(u64, u64, [2, 325, 9375, 28178, 450775, 9780504, 1795265022]);
90#[cfg(target_pointer_width = "64")]
91impl_is_prime!(
92    usize,
93    u64,
94    [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
95);
96#[cfg(target_pointer_width = "32")]
97impl_is_prime!(usize, u32, [2, 7, 61]);
98
99#[cfg(test)]
100mod tests {
101    use super::*;
102
103    #[rustfmt::skip]
104    const PRIME: &[u64] = &[
105        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 47, 53, 59, 61, 67, 73, 79, 83, 89, 97,
106        998244353, 1000000007, 67280421310721, 999999999999999989
107    ];
108
109    const NO_PRIME: &[u64] = &[
110        1,
111        57,
112        5329,
113        49141,
114        4759123141,
115        1122004669633,
116        21652684502221,
117        31858317218647,
118        47636622961201,
119        55245642489451,
120        3071837692357849,
121        3770579582154547,
122        7999252175582851,
123        585226005592931977,
124    ];
125
126    #[rustfmt::skip]
127    const CARMICHAEL: &[u64] = &[
128        561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745,
129        63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461, 252601, 278545, 294409,
130        314821, 334153, 340561, 399001, 410041, 449065, 488881, 512461,
131    ];
132
133    #[test]
134    fn primality_test() {
135        for &p in PRIME {
136            if p <= u8::MAX as u64 {
137                assert!((p as u8).is_prime());
138            }
139            if p <= u16::MAX as u64 {
140                assert!((p as u16).is_prime());
141            }
142            if p <= u32::MAX as u64 {
143                assert!((p as u32).is_prime());
144            }
145            assert!(p.is_prime());
146        }
147    }
148
149    #[test]
150    fn non_primality_test() {
151        for &p in NO_PRIME {
152            if p <= u8::MAX as u64 {
153                assert!(!(p as u8).is_prime());
154            }
155            if p <= u16::MAX as u64 {
156                assert!(!(p as u16).is_prime());
157            }
158            if p <= u32::MAX as u64 {
159                assert!(!(p as u32).is_prime());
160            }
161            assert!(!p.is_prime());
162        }
163
164        for &p in CARMICHAEL {
165            if p <= u8::MAX as u64 {
166                assert!(!(p as u8).is_prime());
167            }
168            if p <= u16::MAX as u64 {
169                assert!(!(p as u16).is_prime());
170            }
171            if p <= u32::MAX as u64 {
172                assert!(!(p as u32).is_prime());
173            }
174            assert!(!p.is_prime());
175        }
176    }
177}
primality_test/lib.rs

primality_test/
lib.rs
