mirror of
https://github.com/RGBCube/uutils-coreutils
synced 2026-01-16 18:21:01 +00:00
6c4e967fc6
This commit builds upon @wikol's Pollard rho implementation. It adds the following: 1. A generator for prime inverse tables. With these, we can do very fast divisibility tests (a single multiply and comparison) for small primes (presently, the first 1000 primes are in the table, which means all numbers of ~26 bits or less can be factored very quickly. 2. Always try prime inverse tables before jumping into Pollard's rho method or using trial division. 3. Since we have eliminated all small factors by the time we're done with the table division, only use slow trial division when the number is big enough to cause overflow issues in Pollard's rho, and jump out of trial division and into Pollard's rho as soon as the number is small enough. 4. Updates the Makefile to regenerate the prime table if it's not up-to-date.
80 lines
1.8 KiB
Rust
80 lines
1.8 KiB
Rust
/*
|
|
* This file is part of the uutils coreutils package.
|
|
*
|
|
* (c) Wiktor Kuropatwa <wiktor.kuropatwa@gmail.com>
|
|
*
|
|
* For the full copyright and license information, please view the LICENSE file
|
|
* that was distributed with this source code.
|
|
*/
|
|
|
|
// computes (a + b) % m using the russian peasant algorithm
|
|
pub fn multiply(mut a: u64, mut b: u64, m: u64) -> u64 {
|
|
let mut result = 0;
|
|
while b > 0 {
|
|
if b & 1 > 0 {
|
|
result = (result + a) % m;
|
|
}
|
|
a = (a << 1) % m;
|
|
b >>= 1;
|
|
}
|
|
result
|
|
}
|
|
|
|
// computes a.pow(b) % m
|
|
fn pow(mut a: u64, mut b: u64, m: u64) -> u64 {
|
|
let mut result = 1;
|
|
while b > 0 {
|
|
if b & 1 > 0 {
|
|
result = multiply(result, a, m);
|
|
}
|
|
a = multiply(a, a, m);
|
|
b >>= 1;
|
|
}
|
|
result
|
|
}
|
|
|
|
fn witness(mut a: u64, exponent: u64, m: u64) -> bool {
|
|
if a == 0 {
|
|
return false;
|
|
}
|
|
if pow(a, m-1, m) != 1 {
|
|
return true;
|
|
}
|
|
a = pow(a, exponent, m);
|
|
if a == 1 {
|
|
return false;
|
|
}
|
|
loop {
|
|
if a == 1 {
|
|
return true;
|
|
}
|
|
if a == m-1 {
|
|
return false;
|
|
}
|
|
a = multiply(a, a, m);
|
|
}
|
|
}
|
|
|
|
// uses deterministic (i.e., fixed witness set) Miller-Rabin test
|
|
pub 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;
|
|
}
|
|
|
|
// 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, 9375, 28178, 450775, 9780504, 1795265022];
|
|
for wit in witnesses.iter() {
|
|
if witness(*wit % num, exponent, num) {
|
|
return false;
|
|
}
|
|
}
|
|
true
|
|
}
|