1
Fork 0
mirror of https://github.com/RGBCube/uutils-coreutils synced 2025-07-30 12:37:49 +00:00

factor::miller_rabin: Extract dividers from the primality test

Another 36% improvement.
This commit is contained in:
nicoo 2020-05-24 15:41:23 +02:00
parent 6b9585b1dc
commit 8241037690
2 changed files with 42 additions and 11 deletions

View file

@ -4,14 +4,29 @@ use crate::numeric::*;
// discovered by Jim Sinclair on 2011-04-20, see miller-rabin.appspot.com // discovered by Jim Sinclair on 2011-04-20, see miller-rabin.appspot.com
const BASIS: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]; const BASIS: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
#[derive(Eq, PartialEq)]
pub(crate) enum Result {
Prime,
Pseudoprime,
Composite(u64),
}
impl Result {
pub(crate) fn is_prime(&self) -> bool {
return *self == Result::Prime;
}
}
// Deterministic Miller-Rabin primality-checking algorithm, adapted to extract // Deterministic Miller-Rabin primality-checking algorithm, adapted to extract
// (some) dividers; it will fail to factor strong pseudoprimes. // (some) dividers; it will fail to factor strong pseudoprimes.
pub(crate) fn is_prime(n: u64) -> bool { pub(crate) fn test(n: u64) -> Result {
use self::Result::*;
if n < 2 { if n < 2 {
return false; return Pseudoprime;
} }
if n % 2 == 0 { if n % 2 == 0 {
return n == 2; return if n == 2 { Prime } else { Composite(2) };
} }
let d = (n - 1).trailing_zeros(); let d = (n - 1).trailing_zeros();
@ -30,7 +45,7 @@ pub(crate) fn is_prime(n: u64) -> bool {
} }
if pow(x, n - 1, n, mul) != 1 { if pow(x, n - 1, n, mul) != 1 {
return false; return Pseudoprime;
} }
x = pow(x, r, n, mul); x = pow(x, r, n, mul);
if x == 1 || x == n - 1 { if x == 1 || x == n - 1 {
@ -40,7 +55,7 @@ pub(crate) fn is_prime(n: u64) -> bool {
loop { loop {
let y = mul(x, x, n); let y = mul(x, x, n);
if y == 1 { if y == 1 {
return false; return Composite(gcd(x - 1, n));
} }
if y == n - 1 { if y == n - 1 {
// This basis element is not a witness of `n` being composite. // This basis element is not a witness of `n` being composite.
@ -51,5 +66,11 @@ pub(crate) fn is_prime(n: u64) -> bool {
} }
} }
true Prime
}
// Used by build.rs' tests
#[allow(dead_code)]
pub(crate) fn is_prime(n: u64) -> bool {
test(n).is_prime()
} }

View file

@ -1,3 +1,4 @@
use crate::miller_rabin::Result::*;
use crate::{miller_rabin, Factors}; use crate::{miller_rabin, Factors};
use numeric::*; use numeric::*;
use rand::distributions::{Distribution, Uniform}; use rand::distributions::{Distribution, Uniform};
@ -39,16 +40,25 @@ fn find_divisor(num: u64) -> u64 {
} }
} }
pub(crate) fn factor(num: u64) -> Factors { pub(crate) fn factor(mut num: u64) -> Factors {
let mut factors = Factors::new(); let mut factors = Factors::new();
if num == 1 { if num == 1 {
return factors; return factors;
} }
if miller_rabin::is_prime(num) { match miller_rabin::test(num) {
factors.push(num); Prime => {
return factors; factors.push(num);
} return factors;
}
Composite(d) => {
num = num / d;
factors *= factor(d);
}
Pseudoprime => {}
};
let divisor = find_divisor(num); let divisor = find_divisor(num);
factors *= factor(divisor); factors *= factor(divisor);