素因数分解(Pollard の rho 法)

GitHub last commit

概要

正の整数 $n$ を受け取り,$n$ の素因数を昇順に列挙する.

制約

計算量

実装

// 👇👇👇👇👇👇👇👇👇👇👇👇 number/factorize_fast 👇👇👇👇👇👇👇👇👇👇👇👇
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
}
// 👆👆👆👆👆👆👆👆👆👆👆👆 number/factorize_fast 👆👆👆👆👆👆👆👆👆👆👆👆

// 👇👇👇👇👇👇👇👇👇👇👇👇 number/is_prime_fast 👇👇👇👇👇👇👇👇👇👇👇👇
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))
}
// 👆👆👆👆👆👆👆👆👆👆👆👆 number/is_prime_fast 👆👆👆👆👆👆👆👆👆👆👆👆

Verified with

参考

戻る