haar_lib/math/prime_test/
miller_rabin.rs

1//! Miller-Rabin素数判定法
2pub use crate::math::prime_test::CheckPrime;
3
4const B: u32 = 64;
5const R: u128 = 1 << B;
6const MASK: u128 = R - 1;
7
8type Montgomery = (u64, u128, u64);
9
10fn montgomery(modulo: u64) -> Montgomery {
11    assert!(modulo % 2 != 0);
12    assert!(modulo > 0);
13
14    let r = R % modulo as u128;
15    let r2 = r * r % modulo as u128;
16    let m = {
17        let mut ret: u64 = 0;
18        let mut r = R;
19        let mut i = 1;
20        let mut t = 0;
21        while r > 1 {
22            if t % 2 == 0 {
23                t += modulo;
24                ret += i;
25            }
26            t >>= 1;
27            r >>= 1;
28            i <<= 1;
29        }
30        ret
31    };
32
33    (modulo, r2, m)
34}
35
36fn reduce(value: u128, modulo: u64, m: u64) -> u64 {
37    let mut ret = (((((value & MASK) * m as u128) & MASK) * modulo as u128 + value) >> B) as u64;
38    if ret >= modulo {
39        ret -= modulo;
40    }
41    ret
42}
43
44fn pow(mut a: u64, mut p: u64, mg: Montgomery) -> u64 {
45    let (modulo, r2, m) = mg;
46
47    let mut value = reduce(r2, modulo, m);
48
49    while p > 0 {
50        if (p & 1) != 0 {
51            value = reduce(value as u128 * a as u128, modulo, m);
52        }
53        a = reduce(a as u128 * a as u128, modulo, m);
54        p >>= 1;
55    }
56
57    value
58}
59
60fn is_composite(a: u64, s: u32, d: u64, mg: Montgomery) -> bool {
61    let (p, r2, m) = mg;
62    let a = reduce(a as u128 * r2, p, m);
63    let pp = reduce((p as u128 - 1) * r2, p, m);
64    let mut x = pow(a, d, mg);
65
66    if reduce(x as u128, p, m) == 1 {
67        false
68    } else {
69        for _ in 0..s {
70            if x == pp {
71                return false;
72            }
73            x = reduce(x as u128 * x as u128, p, m);
74        }
75
76        true
77    }
78}
79
80/// Miller-Rabin素数判定法
81pub struct MillerRabin;
82
83impl CheckPrime<u64> for MillerRabin {
84    fn is_prime(&self, n: u64) -> bool {
85        if n <= 1 {
86            false
87        } else if n == 2 {
88            true
89        } else if n % 2 == 0 {
90            false
91        } else {
92            let s = (n - 1).trailing_zeros();
93            let d = (n - 1) >> s;
94
95            let mg = montgomery(n);
96
97            if n < 4_759_123_141 {
98                ![2, 7, 61]
99                    .into_iter()
100                    .any(|a| a < n && is_composite(a, s, d, mg))
101            } else {
102                ![2, 325, 9375, 28178, 450775, 9780504, 1795265022]
103                    .into_iter()
104                    .any(|a| a < n && is_composite(a, s, d, mg))
105            }
106        }
107    }
108}
haar_lib/math/prime_test/miller_rabin.rs

haar_lib/math/prime_test/
miller_rabin.rs
