Pollard's Rho Algorithm
Pollard rho algorithm is a randomized algorithm which factorizes a number in an average time complexity of \(O(n^{1/4})\).
x.factorize()
factorizes x
and returns a vector with the factors. The order of factors in the vector is undefined.
Miller-Rabin primality test should be with this snippet in the code.
Example
use factorization::*; fn main() { let mut rng = RNG::new(15163487); let a: u64 = 484387724796727379; let mut factors = a.factorize(&mut rng); factors.sort_unstable(); println!("{:?}", factors); // [165551, 2925912406429] println!("{}", factors.iter().product::<u64>()); // 484387724796727379 use millerrabin::Primality; println!("{}", factors.iter().all(|&x| x.is_prime())); } mod factorization { use super::millerrabin::Primality; use std::ops::*; pub trait PollardRho: Primality + From<u8> + PartialOrd + ShrAssign + BitAnd<Output = Self> + Clone { fn rho(self, arr: &mut Vec<Self>, rng: &mut rng::RNG); fn factorize(mut self, rng: &mut rng::RNG) -> Vec<Self> { let mut arr: Vec<Self> = Vec::new(); if self <= 1.into() { return arr; } while self.clone() & 1.into() == 0.into() { self >>= 1.into(); arr.push(2.into()); } self.rho(&mut arr, rng); arr } } macro_rules! impl_pollardrho { ($t:ty, $u:ty, $reset:expr) => { impl PollardRho for $t { fn rho(self, arr: &mut Vec<Self>, rng: &mut rng::RNG) { if self <= 1 { return; } else if self.is_prime() { arr.push(self); return; } let mut i: u64 = 0; let mut x: $t = (rng.next_u64() % self as u64) as $t; let mut y: $t = x; let mut k: u64 = 2; let mut d: $t; let mut reset_limit: u64 = $reset; loop { i += 1; x = (((x as $u * x as $u % self as $u) + (self - 1) as $u) % self as $u) as $t; d = gcd(y.abs_diff(x), self); if d == self || i >= reset_limit { // Reset reset_limit = reset_limit * 3 / 2; i = 0; x = (rng.next_u64() % self as u64) as $t; y = x; } if d != 1 { break; } if i == k { y = x; k <<= 1; } } if d != self { d.rho(arr, rng); (self / d).rho(arr, rng); return; } let mut i = 3; while i * i <= self { if self % i == 0 { i.rho(arr, rng); (d / i).rho(arr, rng); return; } i += 2; } } } }; } impl_pollardrho!(u8, u16, 100000); impl_pollardrho!(u16, u32, 100000); impl_pollardrho!(u32, u64, 100000); impl_pollardrho!(u64, u128, 100000); pub fn gcd<T>(x: T, y: T) -> T where T: Copy + PartialEq + PartialOrd + core::ops::Rem<Output = T> + From<u8>, { if y == 0.into() { x } else { let v = x % y; gcd(y, v) } } pub mod rng { pub struct RNG { val: u64, } impl RNG { pub fn new(seed: u64) -> Self { Self { val: seed } } pub fn next_u64(&mut self) -> u64 { let mut x = self.val; x ^= x << 13; x ^= x >> 7; x ^= x << 17; self.val = x; x } } } pub use rng::*; } mod millerrabin { pub trait Primality { fn is_prime(self) -> bool; } macro_rules! impl_primality { ($t:ty, $u:ty, $thres:expr, $bcnt:expr, $($basis:expr),+) => { impl Primality for $t { fn is_prime(self) -> bool { if self <= 1 { return false; } else if self & 1 == 0 { return self == 2; } const THRES: $t = $thres; const TEST: [$t; $bcnt] = [$($basis,)+]; if self <= THRES { for p in (2..).take_while(|&p| p * p <= self) { if self % p == 0 { return false; } } return true; } let pow = |base: $t, mut exp: $t| -> $t { let mut base = base as $u; let mut ret = 1 as $u; while exp != 0 { if exp & 1 != 0 { ret = (ret * base) % self as $u; } exp >>= 1; base = (base * base) % self as $u; } ret as $t }; let s = (self - 1).trailing_zeros(); let d = (self - 1) >> s; for &a in TEST.iter().take_while(|&&a| a < self - 1) { let mut x = pow(a, d); for _ in 0..s { let y = ((x as $u).pow(2) % self as $u) as $t; if y == 1 && x != 1 && x != self - 1 { return false; } x = y; } if x != 1 { return false; } } true } } }; } impl_primality!(u8, u16, 255, 1, 2); impl_primality!(u16, u32, 2000, 2, 2, 3); impl_primality!(u32, u64, 7000, 3, 2, 7, 61); impl_primality!(u64, u128, 300000, 7, 2, 325, 9375, 28178, 450775, 9780504, 1795265022); }
Code
mod factorization {
use super::millerrabin::Primality;
use std::ops::*;
pub trait PollardRho: Primality + From<u8> + PartialOrd + ShrAssign + BitAnd<Output = Self> + Clone {
fn rho(self, arr: &mut Vec<Self>, rng: &mut rng::RNG);
fn factorize(mut self, rng: &mut rng::RNG) -> Vec<Self> {
let mut arr: Vec<Self> = Vec::new();
if self <= 1.into() {
return arr;
}
while self.clone() & 1.into() == 0.into() {
self >>= 1.into();
arr.push(2.into());
}
self.rho(&mut arr, rng);
arr
}
}
macro_rules! impl_pollardrho {
($t:ty, $u:ty, $reset:expr) => {
impl PollardRho for $t {
fn rho(self, arr: &mut Vec<Self>, rng: &mut rng::RNG) {
if self <= 1 {
return;
} else if self.is_prime() {
arr.push(self);
return;
}
let mut i: u64 = 0;
let mut x: $t = (rng.next_u64() % self as u64) as $t;
let mut y: $t = x;
let mut k: u64 = 2;
let mut d: $t;
let mut reset_limit: u64 = $reset;
loop {
i += 1;
x = (((x as $u * x as $u % self as $u) + (self - 1) as $u) % self as $u) as $t;
d = gcd(y.abs_diff(x), self);
if d == self || i >= reset_limit {
// Reset
reset_limit = reset_limit * 3 / 2;
i = 0;
x = (rng.next_u64() % self as u64) as $t;
y = x;
}
if d != 1 {
break;
}
if i == k {
y = x;
k <<= 1;
}
}
if d != self {
d.rho(arr, rng);
(self / d).rho(arr, rng);
return;
}
let mut i = 3;
while i * i <= self {
if self % i == 0 {
i.rho(arr, rng);
(d / i).rho(arr, rng);
return;
}
i += 2;
}
}
}
};
}
impl_pollardrho!(u8, u16, 100000);
impl_pollardrho!(u16, u32, 100000);
impl_pollardrho!(u32, u64, 100000);
impl_pollardrho!(u64, u128, 100000);
pub fn gcd<T>(x: T, y: T) -> T
where
T: Copy + PartialEq + PartialOrd + core::ops::Rem<Output = T> + From<u8>,
{
if y == 0.into() {
x
} else {
let v = x % y;
gcd(y, v)
}
}
pub mod rng {
pub struct RNG {
val: u64,
}
impl RNG {
pub fn new(seed: u64) -> Self {
Self { val: seed }
}
pub fn next_u64(&mut self) -> u64 {
let mut x = self.val;
x ^= x << 13;
x ^= x >> 7;
x ^= x << 17;
self.val = x;
x
}
}
}
pub use self::rng::*;
}
Last modified on 231008.