use crate::{zz, ZZ, Result};
/// Generated using the following sage code
/// ```sage
/// s = []
/// for j in range(2**10): s.append(int(''.join(['1' if is_prime(2*(j*8+x)+1) else '0' for x in range(8)])[::-1],2))
/// print(s)
/// ```
const LOW_PRIMES: [u8; 1024] = [110, 203, 180, 100, 154, 18, 109, 129, 50, 76, 74, 134, 13, 130, 150, 33, 201, 52, 4, 90, 32, 97, 137, 164, 68, 17, 134, 41, 209, 130, 40, 74, 48, 64, 66, 50, 33, 153, 52, 8, 75, 6, 37, 66, 132, 72, 138, 20, 5, 66, 48, 108, 8, 180, 64, 11, 160, 8, 81, 18, 40, 137, 4, 101, 152, 48, 76, 128, 150, 68, 18, 128, 33, 66, 18, 65, 201, 4, 33, 192, 50, 45, 152, 0, 0, 73, 4, 8, 129, 150, 104, 130, 176, 37, 8, 34, 72, 137, 162, 64, 89, 38, 4, 144, 6, 64, 67, 48, 68, 146, 0, 105, 16, 130, 8, 8, 164, 13, 65, 18, 96, 192, 0, 36, 210, 34, 97, 8, 132, 4, 27, 130, 1, 211, 16, 1, 2, 160, 68, 192, 34, 96, 145, 20, 12, 64, 166, 4, 210, 148, 32, 9, 148, 32, 82, 0, 8, 16, 162, 76, 0, 130, 1, 81, 16, 8, 139, 164, 37, 154, 48, 68, 129, 16, 76, 3, 2, 37, 82, 128, 8, 73, 132, 32, 80, 50, 0, 24, 162, 64, 17, 36, 40, 1, 132, 1, 1, 160, 65, 10, 18, 69, 0, 54, 8, 0, 38, 41, 131, 130, 97, 192, 128, 4, 16, 16, 109, 0, 34, 72, 88, 38, 12, 194, 16, 72, 137, 36, 32, 88, 32, 69, 136, 36, 0, 25, 2, 37, 192, 16, 104, 8, 20, 1, 202, 50, 40, 128, 0, 4, 75, 38, 0, 19, 144, 96, 130, 128, 37, 208, 0, 1, 16, 50, 12, 67, 134, 33, 17, 0, 8, 67, 36, 4, 72, 16, 12, 144, 146, 0, 67, 32, 45, 0, 6, 9, 136, 36, 64, 192, 50, 9, 9, 130, 0, 83, 128, 8, 128, 150, 65, 129, 0, 64, 72, 16, 72, 8, 150, 72, 88, 32, 41, 195, 128, 32, 2, 148, 96, 146, 0, 32, 129, 34, 68, 16, 160, 5, 64, 144, 1, 73, 32, 4, 10, 0, 36, 137, 52, 72, 19, 128, 44, 192, 130, 41, 0, 36, 69, 8, 0, 8, 152, 54, 4, 82, 132, 4, 208, 4, 0, 138, 144, 68, 130, 50, 101, 24, 144, 0, 10, 2, 1, 64, 2, 40, 64, 164, 4, 146, 48, 4, 17, 134, 8, 66, 0, 44, 82, 4, 8, 201, 132, 96, 72, 18, 9, 153, 36, 68, 0, 36, 0, 3, 20, 33, 0, 16, 1, 26, 50, 5, 136, 32, 64, 64, 6, 9, 195, 132, 64, 1, 48, 96, 24, 2, 104, 17, 144, 12, 2, 162, 4, 0, 134, 41, 137, 20, 36, 130, 2, 65, 8, 128, 4, 25, 128, 8, 16, 18, 104, 66, 164, 4, 0, 2, 97, 16, 6, 12, 16, 0, 1, 18, 16, 32, 3, 148, 33, 66, 18, 101, 24, 148, 12, 10, 4, 40, 1, 20, 41, 10, 164, 64, 208, 0, 64, 1, 144, 4, 65, 32, 45, 64, 130, 72, 193, 32, 0, 16, 48, 1, 8, 36, 4, 89, 132, 36, 0, 2, 41, 130, 0, 97, 88, 2, 72, 129, 22, 72, 16, 0, 33, 17, 6, 0, 202, 160, 64, 2, 0, 4, 145, 176, 0, 66, 4, 12, 129, 6, 9, 72, 20, 37, 146, 32, 37, 17, 160, 0, 10, 134, 12, 193, 2, 72, 0, 32, 69, 8, 50, 0, 152, 6, 4, 19, 34, 0, 130, 4, 72, 129, 20, 68, 130, 18, 36, 24, 16, 64, 67, 128, 40, 208, 4, 32, 129, 36, 100, 216, 0, 44, 9, 18, 8, 65, 162, 0, 0, 2, 65, 202, 32, 65, 192, 16, 1, 24, 164, 4, 24, 164, 32, 18, 148, 32, 131, 160, 64, 2, 50, 68, 128, 4, 0, 24, 0, 12, 64, 134, 96, 138, 0, 100, 136, 18, 5, 1, 130, 0, 74, 162, 1, 193, 16, 97, 9, 4, 1, 136, 0, 96, 1, 180, 64, 8, 6, 1, 3, 128, 8, 64, 148, 4, 138, 32, 41, 128, 2, 12, 82, 2, 1, 66, 132, 0, 128, 132, 100, 2, 50, 72, 0, 48, 68, 64, 34, 33, 0, 2, 8, 195, 160, 4, 208, 32, 64, 24, 22, 64, 64, 0, 40, 82, 144, 8, 130, 20, 1, 24, 16, 8, 9, 130, 64, 10, 160, 32, 147, 128, 8, 192, 0, 32, 82, 0, 5, 1, 16, 64, 17, 6, 12, 130, 0, 0, 75, 144, 68, 154, 0, 40, 128, 144, 4, 74, 6, 9, 67, 2, 40, 0, 52, 1, 24, 0, 101, 9, 128, 68, 3, 0, 36, 2, 130, 97, 72, 20, 65, 0, 18, 40, 0, 52, 8, 81, 4, 5, 18, 144, 40, 137, 132, 96, 18, 16, 73, 16, 38, 64, 73, 130, 0, 145, 16, 1, 10, 36, 64, 136, 16, 76, 16, 4, 0, 80, 162, 44, 64, 144, 72, 10, 176, 1, 80, 18, 8, 0, 164, 4, 9, 160, 40, 146, 2, 0, 67, 16, 33, 2, 32, 65, 129, 50, 0, 8, 4, 12, 82, 0, 33, 73, 132, 32, 16, 2, 1, 129, 16, 72, 64, 34, 1, 1, 132, 105, 193, 48, 1, 200, 2, 68, 136, 0, 12, 1, 2, 45, 192, 18, 97, 0, 160, 0, 192, 48, 64, 1, 18, 8, 11, 32, 0, 128, 148, 64, 1, 132, 64, 0, 50, 0, 16, 132, 0, 11, 36, 0, 1, 6, 41, 138, 132, 65, 128, 16, 8, 8, 148, 76, 3, 128, 1, 64, 150, 64, 65, 32, 32, 80, 34, 37, 137, 162, 64, 64, 164, 32, 2, 134, 40, 1, 32, 33, 74, 16, 8, 0, 20, 8, 64, 4, 37, 66, 2, 33, 67, 16, 4, 146, 0, 33, 17, 160, 76, 24, 34, 9, 3, 132, 65, 137, 16, 4, 130, 34, 36, 1, 20, 8, 8, 132, 8, 193, 0, 9, 66, 176, 65, 138, 2, 0, 128, 54, 4, 73, 160, 36, 145, 0, 0, 2, 148, 65, 146, 2, 1, 8, 6, 8, 9, 0, 1, 208, 22, 40, 137, 128, 96, 0, 0, 104, 1, 144, 12, 80, 32, 1, 64, 128, 64, 66, 48, 65];
pub fn is_prime(p: impl AsRef<ZZ>) -> bool {
let p = p.as_ref();
if p <= 1 {
return false;
} else if p == 2 {
return true;
} else if p & 1 == 0 {
return false;
}
let t: Result<u64> = p.try_into();
match t {
Ok(p) => {
if p < 16385 {
let mut p: usize = p.try_into().unwrap();
p >>= 1;
let idx = p>>3;
return (LOW_PRIMES[idx]>>(p&7)) & 1 == 1;
}
return primal::is_prime(p);
},
Err(_) => {}
}
miller_rabin(p, 100)
}
pub fn miller_rabin(n: impl AsRef<ZZ>, k: usize) -> bool {
let n = n.as_ref();
if n == 2 || n == 3 || n == 5 {
return true;
}
if n <= 1 || n % 2 == 0 {
return false;
}
let mut d = n - 1;
while &d & 1 == 0 {
d /= 2;
}
for _ in 0..k {
let a = ZZ::rand_range(&zz!(2), n - 4);
let mut x = a.mod_pow(&d, n);
if x == 1 || x == n - 1 {
continue;
}
let mut ok = false;
while d != n - 1 {
x = (&x * &x) % n;
d *= 2;
if x == 1 {
return false;
}
if x == n - 1 {
ok = true;
break;
}
}
if !ok {
return false;
}
}
return true;
}
#[cfg(test)]
mod test {
use super::{miller_rabin, ZZ, is_prime, zz};
#[test]
fn test_miller_rabin() {
let primes = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
89, 97,
];
for p in primes {
assert_eq!(miller_rabin(ZZ::from(p), 4), true);
}
let non_primes = [
1, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34,
35, 36, 38, 39, 40, 42, 44, 45, 46, 48, 49, 50, 51, 52, 54, 55, 56, 57, 58, 60, 62, 63,
64, 65, 66, 68, 69, 70, 72, 74, 75, 76, 77, 78, 80, 81, 82, 84, 85, 86, 87, 88,
];
for p in non_primes {
assert_eq!(miller_rabin(ZZ::from(p), 4), false);
}
}
#[test]
fn test_large() {
let p = zz!(8046047750471549884742864261737478411129806839855555228638804849413140951997507220381046990534016609);
assert_eq!(is_prime(&p), true);
assert_eq!(is_prime(p+2), false);
}
}