forked from TheAlgorithms/Rust
-
Notifications
You must be signed in to change notification settings - Fork 0
/
miller_rabin.rs
105 lines (96 loc) · 3.69 KB
/
miller_rabin.rs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
fn modulo_power(mut base: u64, mut power: u64, modulo: u64) -> u64 {
base %= modulo;
if base == 0 {
return 0; // return zero if base is divisible by modulo
}
let mut ans: u128 = 1;
let mut bbase: u128 = base as u128;
while power > 0 {
if (power % 2) == 1 {
ans = (ans * bbase) % (modulo as u128);
}
bbase = (bbase * bbase) % (modulo as u128);
power /= 2;
}
ans as u64
}
fn check_prime_base(number: u64, base: u64, two_power: u64, odd_power: u64) -> bool {
// returns false if base is a witness
let mut x: u128 = modulo_power(base, odd_power, number) as u128;
let bnumber: u128 = number as u128;
if x == 1 || x == (bnumber - 1) {
return true;
}
for _ in 1..two_power {
x = (x * x) % bnumber;
if x == (bnumber - 1) {
return true;
}
}
false
}
pub fn miller_rabin(number: u64, bases: &[u64]) -> u64 {
// returns zero on a probable prime, and a witness if number is not prime
// A base set for deterministic performance on 64 bit numbers is:
// [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
// another one for 32 bits:
// [2, 3, 5, 7], with smallest number to fail 3'215'031'751 = 151 * 751 * 28351
// note that all bases should be prime
if number <= 4 {
match number {
0 => {
panic!("0 is invalid input for Miller-Rabin. 0 is not prime by definition, but has no witness");
}
2 => return 0,
3 => return 0,
_ => return number,
}
}
if bases.contains(&number) {
return 0;
}
let two_power: u64 = (number - 1).trailing_zeros() as u64;
let odd_power = (number - 1) >> two_power;
for base in bases {
if !check_prime_base(number, *base, two_power, odd_power) {
return *base;
}
}
0
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn basic() {
let default_bases = vec![2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
// these bases make miller rabin deterministic for any number < 2 ^ 64
// can use smaller number of bases for deterministic performance for numbers < 2 ^ 32
assert_eq!(miller_rabin(3, &default_bases), 0);
assert_eq!(miller_rabin(7, &default_bases), 0);
assert_eq!(miller_rabin(11, &default_bases), 0);
assert_eq!(miller_rabin(2003, &default_bases), 0);
assert_ne!(miller_rabin(1, &default_bases), 0);
assert_ne!(miller_rabin(4, &default_bases), 0);
assert_ne!(miller_rabin(6, &default_bases), 0);
assert_ne!(miller_rabin(21, &default_bases), 0);
assert_ne!(miller_rabin(2004, &default_bases), 0);
// bigger test cases.
// primes are generated using openssl
// non primes are randomly picked and checked using openssl
// primes:
assert_eq!(miller_rabin(3629611793, &default_bases), 0);
assert_eq!(miller_rabin(871594686869, &default_bases), 0);
assert_eq!(miller_rabin(968236663804121, &default_bases), 0);
assert_eq!(miller_rabin(6920153791723773023, &default_bases), 0);
// random non primes:
assert_ne!(miller_rabin(4546167556336341257, &default_bases), 0);
assert_ne!(miller_rabin(4363186415423517377, &default_bases), 0);
assert_ne!(miller_rabin(815479701131020226, &default_bases), 0);
// these two are made of two 31 bit prime factors:
// 1950202127 * 2058609037 = 4014703722618821699
assert_ne!(miller_rabin(4014703722618821699, &default_bases), 0);
// 1679076769 * 2076341633 = 3486337000477823777
assert_ne!(miller_rabin(3486337000477823777, &default_bases), 0);
}
}