Problem
Extract an arbitrary non-trivial factor of a given integer \(n\).
Prerequisite
Pollard-Rho
Let’s consider a pseudo-random polynomial function \(f_n : Z_n \rightarrow Z_n\), where \(f_n(x) \equiv x^2 + 1\pmod n\).
Consider starting with a random seed \(s_{n, 0}\). We can define the sequence in the form of linear recurrence: \(s_{n, i} = f_n(s_{n, i - 1})\).
The Birthday Paradox states that, given that there are \(d\) days in a year, with \(O(n^{0.5})\) people there is roughly a \(50\%\) chance that there is at least a pair with the same birthday. Since we may assume that \(s_n\) is a random sequence, the birthday paradox enlightens that it is likely there is some pair of elements in the first \(O(n^{0.5})\) elements of the sequence. The distance between them is the period of the sequence. You may think of the Greek letter \(\rho\), where the cycle has length of the distance.
Up to now, we do not see anything that directly serves the purpose, and the complexity (\(O(n^{0.5})\)) is not good.
Let us consider \(n = p\times q\), where \(p > 1 \land q > 1\). Without loss of generality, we may assume that \(p \leq q\).
Instead of looking at sequence \(s_n\), let us see what \(s_p\) can do for us. We may discover that the expected length of cycle of the sequence \(s_p\) is \(O(p^{0.5}) \leq O(n^{0.25})\). Although we have no clue what \(p\) could be at this moment, we are still able to proceed, since we can expect that by \(O(n^{0.25})\) cost we can find some two distinct indices \(i\) and \(j\) that satisfies \(s_{p, i} \equiv s_{p, j} \pmod p\). This implies that \(s_{n, i} \equiv s_{n, j} \pmod p\), where the probability that \(s_{n, i} \equiv s_{n, j} \pmod n\) is \(q^{-1}\) (assuming \(s\) is somewhat uniformly distributed). Hence, there is only \(q^{-1} \leq n^{-0.5}\) chance that the algorithm fails, and otherwise \(gcd(s_{n, i}, s_{n, j})\) is a non-trivial factor of \(n\) that is a multiple of \(p\). There is no worry even in case it failed, as we can start again with a different random seed.
Note that the algorithm will not stop if \(n\) is prime, so we may need primality test algorithm such as Miller-Rabin to support Pollard’s algorithm.
Snippet
struct PrimeFactor {
static long long mul(long long x, long long y, long long mod) {
long long m = x, s = 0;
for (; y; y >>= 1, m <<= 1, m = m >= mod ? m - mod : m)
if (y & 1) s += m, s = s >= mod ? s - mod : s;
return s;
}
static long long powmod(long long x, long long p, long long mod) {
long long s = 1, m = x % mod;
for (; p; m = mul(m, m, mod), p >>= 1)
if (p & 1) s = mul(s, m, mod);
return s;
}
static bool miller_rabin(long long n, int s = 7) {
static const long long wits[7] = {
2, 325, 9375, 28178, 450775, 9780504, 1795265022
};
auto witness = [&](long long a, long long n, long long u, int t) {
long long x = powmod(a, u, n), nx;
for (int i = 0; i < t; ++i, x = nx) {
nx = mul(x, x, n);
if (nx == 1 && x != 1 && x != n - 1) return true;
}
return x != 1;
};
if (n < 2) return 0;
if (~n & 1) return n == 2;
long long u = n - 1, t = 0, a;
for (; ~u & 1; u >>= 1) ++t;
while (s--) if ((a = wits[s] % n) && witness(a, n, u, t)) return false;
return true;
}
static long long pollard_rho(long long n) {
auto f = [](long long x, long long n) { return mul(x, x, n) + 1; };
if (~n & 1) return 2;
while (true) {
long long x = rand() % (n - 1) + 1, y = 2, d = 1;
for (int sz = 2; d == 1; y = x, sz <<= 1)
for (int i = 0; i < sz && d <= 1; ++i)
x = f(x, n), d = __gcd(abs(x - y), n);
if (d && n - d) return d;
}
}
static vector<pair<long long, int>> factorize(long long n) {
set<long long> factors;
function<void(long long)> rec = [&](long long x) {
if (x == 1) return;
if (miller_rabin(x)) return void(factors.emplace(x));
auto u = pollard_rho(x);
rec(u);
for (auto p : factors) while (x % p == 0) x /= p;
rec(x);
};
rec(n);
vector<pair<long long, int>> res;
for (auto p : factors) {
res.emplace_back(p, 0);
for (; n % p == 0; n /= p) ++res.back().second;
}
return res;
}
};