1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

//! Performs Pollard's rho algorithm.
//!
//! # Example
//!
//! ```
//! assert_eq!([1_162_193, 1_347_377], *rho::factorize(1_565_912_117_761));
//! ```

/// Performs Pollard's rho algorithm.
pub fn factorize(n: u64) -> Vec<u64> {
    if n <= 1 {
        return vec![];
    }

    if miller_rabin::is_prime(n) {
        return vec![n];
    }

    let factor = rho(n);
    let mut ret = factorize(factor);
    ret.extend(factorize(n / factor));
    ret.sort_unstable();
    ret
}

#[allow(clippy::many_single_char_names)]
fn rho(n: u64) -> u64 {
    debug_assert!(!miller_rabin::is_prime(n));

    if n % 2 == 0 {
        return 2;
    }

    for c in 1.. {
        let add = |lhs: u64, rhs: u64| -> _ {
            let mut ret = lhs + rhs;
            if ret >= n {
                ret -= n;
            }
            ret
        };

        let sub = |lhs: u64, rhs: u64| -> _ {
            if lhs < rhs {
                n + lhs - rhs
            } else {
                lhs - rhs
            }
        };

        let mul = |lhs: u64, rhs: u64| -> _ {
            if let Some(mul) = lhs.checked_mul(rhs) {
                mul % n
            } else {
                (u128::from(lhs) * u128::from(rhs) % u128::from(n)) as _
            }
        };

        let g = |x: u64| add(mul(x, x), c);

        let mut x = 2;
        let mut y = 2;
        let d = loop {
            x = g(x);
            y = g(g(y));
            let d = gcd::gcd(sub(x, y), n);
            if d > 1 {
                break d;
            }
        };
        if d < n {
            return d;
        }
    }
    unreachable!();
}

#[cfg(test)]
mod tests {
    #[test]
    fn test() {
        test(0, &[]);
        test(1, &[]);
        test(2, &[2]);
        test(3, &[3]);
        test(4, &[2, 2]);
        test(5, &[5]);
        test(6, &[2, 3]);
        test(7, &[7]);
        test(8, &[2, 2, 2]);
        test(9, &[3, 3]);
        test(10, &[2, 5]);

        fn test(n: u64, expected: &[u64]) {
            assert_eq!(*expected, *super::factorize(n));
        }
    }
}