-
Notifications
You must be signed in to change notification settings - Fork 687
/
Factor.h
71 lines (69 loc) · 3.02 KB
/
Factor.h
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
/**
* Author: chilli, SJTU, pajenegod
* Date: 2020-03-04
* License: CC0
* Source: own
* Description: Pollard-rho randomized factorization algorithm. Returns prime
* factors of a number, in arbitrary order (e.g. 2299 -> \{11, 19, 11\}).
* Time: $O(n^{1/4})$, less for numbers with small factors.
* Status: stress-tested
*
* Details: This implementation uses the improvement described here
* (https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm#Variants), where
* one can accumulate gcd calls by some factor (40 chosen here through
* exhaustive testing). This improves performance by approximately 6-10x
* depending on the inputs and speed of gcd. Benchmark found here:
* (https://ideone.com/nGGD9T)
*
* GCD can be improved by a factor of 1.75x using Binary GCD
* (https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/).
* However, with the gcd accumulation the bottleneck moves from the gcd calls
* to the modmul. As GCD only constitutes ~12% of runtime, speeding it up
* doesn't matter so much.
*
* This code can probably be sped up by using a faster mod mul - potentially
* montgomery reduction on 128 bit integers.
* Alternatively, one can use a quadratic sieve for an asymptotic improvement,
* which starts being faster in practice around 1e13.
*
* Brent's cycle finding algorithm was tested, but doesn't reduce modmul calls
* significantly.
*
* Subtle implementation notes:
* - prd starts off as 2 to handle the case n = 4; it's harmless for other n
* since we're guaranteed that n > 2. (Pollard rho has problems with prime
* powers in general, but all larger ones happen to work.)
* - t starts off as 30 to make the first gcd check come earlier, as an
* optimization for small numbers.
* - we vary f between restarts because the cycle finding algorithm does not
* find the first element in the cycle but rather one at distance k*|cycle|
* from the start, and that can result in continual failures if all cycles
* have the same size for all prime factors. E.g. fixing f(x) = x^2 + 1 would
* loop infinitely for n = 352523 * 352817, where all cycles have size 821.
* - we operate on residues in [i, n + i) which modmul is not designed to
* handle, but specifically modmul(x, x) still turns out to work for small
* enough i. (With reference to the proof in modmul-proof.tex, the argument
* for "S is in [-c, 2c)" goes through unchanged, while S < 2^63 now follows
* from S < 2c and S = x^2 (mod c) together implying S < c + i^2.)
*/
#pragma once
#include "ModMulLL.h"
#include "MillerRabin.h"
ull pollard(ull n) {
ull x = 0, y = 0, t = 30, prd = 2, i = 1, q;
auto f = [&](ull x) { return modmul(x, x, n) + i; };
while (t++ % 40 || __gcd(prd, n) == 1) {
if (x == y) x = ++i, y = f(x);
if ((q = modmul(prd, max(x,y) - min(x,y), n))) prd = q;
x = f(x), y = f(f(y));
}
return __gcd(prd, n);
}
vector<ull> factor(ull n) {
if (n == 1) return {};
if (isPrime(n)) return {n};
ull x = pollard(n);
auto l = factor(x), r = factor(n / x);
l.insert(l.end(), all(r));
return l;
}