lib_modulo/
prime.rs

1use crate::Modulus64;
2
3// from <https://miller-rabin.appspot.com/>
4/// x < 350_269_456_337
5static SET3: [u64; 3] = [
6    4230279247111683200,
7    14694767155120705706,
8    16641139526367750375,
9];
10/// x < 7_999_252_175_582_851
11static SET5: [u64; 5] = [
12    2,
13    4130806001517,
14    149795463772692060,
15    186635894390467037,
16    3967304179347715805,
17];
18/// x < 2^64
19static SET7: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
20
21/// Performs deterministic Miller-Rabin primality test.
22///
23/// # Time complexity
24///
25/// *O*(log *x*)
26pub const fn primality_test(x: u64) -> bool {
27    /// `x` is prime iff `(test >> x) & 1 == 1`
28    const SMALL_PRIME: u64 = {
29        let mut test = 0;
30        let mut x = 64;
31        while x > 0 {
32            x -= 1;
33
34            test = (test << 1) | if primality_test_naive(x) { 1 } else { 0 };
35        }
36        test
37    };
38
39    /// remove multiples of 2, 3 or 5
40    const MAY_BE_PRIME_30: u32 = {
41        let mut table = 0;
42
43        let mut n = 0;
44        while n < 30 {
45            table |= if n % 2 == 0 || n % 3 == 0 || n % 5 == 0 {
46                0 // composite
47            } else {
48                1 // may be prime
49            } << n;
50            n += 1;
51        }
52
53        table
54    };
55
56    if x < 64 {
57        return (SMALL_PRIME >> x) & 1 == 1;
58    } else if (MAY_BE_PRIME_30 >> x % 30) & 1 == 0 || x % 7 == 0 {
59        return false;
60    }
61
62    let (s, d) = {
63        let x = x - 1;
64        let s = x.trailing_zeros();
65        (s - 1, x >> s)
66    };
67
68    let modulus = Modulus64::new(x);
69    let one = modulus.residue(1).x;
70    // (a - a) r = 0 (mod x), r % x != 0
71    let neg_one = x - one;
72
73    let witness = if x < 350_269_456_337 {
74        SET3.as_slice()
75    } else if x < 7_999_252_175_582_851 {
76        SET5.as_slice()
77    } else {
78        SET7.as_slice()
79    };
80
81    let mut i = 0;
82    'test: while i < witness.len() {
83        let mut mint = modulus.residue(witness[i]);
84        i += 1;
85
86        if mint.is_zero() {
87            continue;
88        }
89
90        mint = mint.pow(d);
91        if mint.x == one || mint.x == neg_one {
92            continue;
93        }
94
95        let mut s = s;
96        while s > 0 {
97            s -= 1;
98
99            mint.x = mint.modulus.mul(mint.x, mint.x);
100            if mint.x == neg_one {
101                continue 'test;
102            }
103        }
104
105        return false;
106    }
107
108    true
109}
110
111const fn primality_test_naive(x: u64) -> bool {
112    if x < 2 {
113        return false;
114    }
115
116    let mut d = 1;
117    while d < x.isqrt() {
118        d += 1;
119
120        if x % d == 0 {
121            return false;
122        }
123    }
124
125    true
126}
127
128// #[test]
129// fn f() {
130//     use std::io::Write;
131
132//     let mut f = std::fs::File::create("./src/small_prime_context_u16_raw.rs").unwrap();
133
134//     let _ = f.write(b"[\n");
135//     for n in (3..1 << 16).step_by(2) {
136//         if primality_test(n) {
137//             let modulus = Context64::new(n);
138
139//             let _ = f.write(format!("({}, {}, {}),\n", modulus.n, modulus.inv_n, modulus.r2_mod_n).as_bytes());
140//         }
141//     }
142//     let _ = f.write(b"]");
143// }
144
145#[cfg(test)]
146mod tests {
147    use rand::{rng, Rng};
148
149    use super::*;
150
151    #[test]
152    fn small() {
153        for x in 0..500_000 {
154            assert_eq!(primality_test(x), primality_test_naive(x), "{x}")
155        }
156    }
157
158    #[test]
159    fn intermediate() {
160        let mut rng = rng();
161
162        for x in std::iter::repeat_with(|| rng.random_range(1 << 30..1 << 40)).take(100) {
163            assert_eq!(primality_test(x), primality_test_naive(x), "{x}")
164        }
165    }
166}
lib_modulo/prime.rs

lib_modulo/
prime.rs
