factor::miller_rabbin: Refactor before extracting dividers
Replace iterated division with u64::trailing_zeros, hoist the selection of `mul`
out of the loop, another cool 49.5% runtime improvement.
use crate::numeric::*;
-fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
- if a == 0 {
+// Small set of bases for the Miller-Rabin prime test, valid for all 64b integers;
+// discovered by Jim Sinclair on 2011-04-20, see miller-rabin.appspot.com
+const BASIS: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
+
+// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
+// (some) dividers; it will fail to factor strong pseudoprimes.
+pub(crate) fn is_prime(n: u64) -> bool {
+ if n < 2 {
return false;
}
+ if n % 2 == 0 {
+ return n == 2;
+ }
- let mul = if m < 1 << 63 {
+ let d = (n - 1).trailing_zeros();
+ let r = (n - 1) >> d;
+
+ let mul = if n < 1 << 63 {
sm_mul as fn(u64, u64, u64) -> u64
} else {
big_mul as fn(u64, u64, u64) -> u64
};
- if pow(a, m - 1, m, mul) != 1 {
- return true;
- }
- a = pow(a, exponent, m, mul);
- if a == 1 {
- return false;
+ for a in BASIS.iter() {
+ let mut x = a % n;
+ if x == 0 {
+ break;
}
- loop {
- if a == 1 {
- return true;
- }
- if a == m - 1 {
+
+ if pow(x, n - 1, n, mul) != 1 {
return false;
}
- a = mul(a, a, m);
- }
+ x = pow(x, r, n, mul);
+ if x == 1 || x == n - 1 {
+ break;
}
-// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
-// (some) dividers; it will fail to factor strong pseudoprimes.
-pub(crate) fn is_prime(num: u64) -> bool {
- if num < 2 {
+ loop {
+ let y = mul(x, x, n);
+ if y == 1 {
return false;
}
- if num % 2 == 0 {
- return num == 2;
+ if y == n - 1 {
+ // This basis element is not a witness of `n` being composite.
+ // Keep looking.
+ break;
+ }
+ x = y;
}
- let mut exponent = num - 1;
- while exponent & 1 == 0 {
- exponent >>= 1;
}
- // These witnesses detect all composites up to at least 2^64.
- // Discovered by Jim Sinclair, according to http://miller-rabin.appspot.com
- let witnesses = [2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022];
- !witnesses
- .iter()
- .any(|&wit| witness(wit % num, exponent, num))
+ true
}
