fn gcd(x: u64, y: u64) -> u64 {
if y == 0 {
return x;
}
gcd(y, x % y)
}
fn factorize_inner(n: u64) -> Vec<u64> {
assert!(n != 0);
if n == 1 {
return vec![];
}
if is_prime(n) {
return vec![n];
}
let find_divisor_of_n = |seed: u64| -> Option<u64> {
let f = |x| ((x as u128 * x as u128 + seed as u128) % n as u128) as u64;
let mut x = seed;
let mut y = f(x);
let mut d = 1;
while d == 1 {
x = f(x);
y = f(f(y));
d = gcd(x.abs_diff(y), n);
}
Some(d).filter(|&d| d != n)
};
let divisor = if n % 2 == 0 {
2
} else {
(1..).find_map(find_divisor_of_n).unwrap()
};
let mut res = factorize_inner(divisor);
res.extend(factorize_inner(n / divisor));
res
}
fn factorize(n: u64) -> Vec<u64> {
let mut res = factorize_inner(n);
res.sort();
res
}
fn pow_mod(mut x: u64, mut n: u64, modulus: u64) -> u64 {
let mut res = 1;
while n > 0 {
if n & 1 == 1 {
res = (res as u128 * x as u128 % modulus as u128) as u64;
}
x = (x as u128 * x as u128 % modulus as u128) as u64;
n >>= 1;
}
res
}
fn is_prime(n: u64) -> bool {
if n <= 1 {
return false;
}
let s = (n - 1).trailing_zeros();
let d = (n - 1) >> s;
let is_witness_of_n = |a: u64| -> bool {
if a == 0 {
return false;
}
let mut x = pow_mod(a, d, n);
if x == 1 {
return false;
}
for _ in 0..s {
if x == n - 1 {
return false;
}
x = (x as u128 * x as u128 % n as u128) as u64;
}
true
};
let a = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
a.iter().all(|&a| !is_witness_of_n(a % n))
}