mirror of
https://github.com/RGBCube/uutils-coreutils
synced 2025-07-30 12:37:49 +00:00
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.
This commit is contained in:
parent
e3ecc81d97
commit
6b9585b1dc
1 changed files with 38 additions and 35 deletions
|
@ -1,52 +1,55 @@
|
|||
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;
|
||||
}
|
||||
loop {
|
||||
if a == 1 {
|
||||
return true;
|
||||
for a in BASIS.iter() {
|
||||
let mut x = a % n;
|
||||
if x == 0 {
|
||||
break;
|
||||
}
|
||||
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 {
|
||||
return false;
|
||||
}
|
||||
if num % 2 == 0 {
|
||||
return num == 2;
|
||||
}
|
||||
let mut exponent = num - 1;
|
||||
while exponent & 1 == 0 {
|
||||
exponent >>= 1;
|
||||
loop {
|
||||
let y = mul(x, x, n);
|
||||
if y == 1 {
|
||||
return false;
|
||||
}
|
||||
if y == n - 1 {
|
||||
// This basis element is not a witness of `n` being composite.
|
||||
// Keep looking.
|
||||
break;
|
||||
}
|
||||
x = y;
|
||||
}
|
||||
}
|
||||
|
||||
// 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
|
||||
}
|
||||
|
|
Loading…
Add table
Add a link
Reference in a new issue