1. Formal Invariant — Non-Trivial Square Roots of Unity & Birthday Collisions
Invariant 1: Strong Pseudoprimality (Miller-Rabin Witness)
For an odd integer $N > 2$, write $N - 1 = 2^s d$ with $d$ odd. If $N$ is prime and $\gcd(a, N) = 1$, then either:
$$a^d \equiv 1 \pmod N \quad \text{or} \quad \exists\, r \in [0, s-1] \text{ such that } a^{2^r d} \equiv -1 \pmod N$$
Why? Because if $x = a^{2^r d}$ satisfies $x^2 \equiv 1 \pmod N$ for a prime $N$, then $N \mid (x - 1)(x + 1)$, which forces $x \equiv \pm 1 \pmod N$. Any base $a$ violating this condition is an irrefutable witness that $N$ is composite (and overcomes Carmichael numbers like $561 = 3 \times 11 \times 17$). For every odd composite $N$, at least $3/4$ of bases $a \in [1, N-1]$ are witnesses (Rabin 1980).
Invariant 2: Pollard's $\rho$ Implicit Modulo-$p$ Cycle
Let $p \le \sqrt{N}$ be a prime factor of composite $N$, and iterate $f(x) = (x^2 + c) \bmod N$. Because $(x^2 + c \bmod N) \bmod p = (x^2 + c) \bmod p$, the sequence implicitly performs a pseudo-random walk in $\mathbb{Z}/p\mathbb{Z}$! By the Birthday Paradox, two states $x, y$ satisfy $x \equiv y \pmod p$ after expected $O(\sqrt{p}) \le O(N^{1/4})$ steps, at which point $\gcd(|x - y|, N) > 1$ extracts a non-trivial divisor of $N$.
2. Complete 64-Bit Factorization Pipeline
with Batched GCD to find divisor d"] RunRho --> Recurse["Recurse Factor(d) and Factor(N / d)"]
3. Hardware & Algorithmic Optimizations: __int128 & Brent's Batched GCD
Avoiding Overflow (__int128_t) & Amortizing gcd via Richard Brent's Batching
Multiplying two 64-bit integers $a \cdot b \pmod N$ when $N \sim 10^{18}$ overflows uint64_t! Using GCC/Clang's hardware-backed static_cast<__uint128_t>(a) * b % N compiles directly to a single x86-64 mulq + divq instruction pair. Furthermore, calling std::gcd(|x - y|, N) on every step of Pollard's $\rho$ adds an $O(\log N)$ factor. **Brent's optimization** accumulates the product $P = \prod_{j=1}^{128} |x_j - y| \bmod N$ across $128$ steps and computes std::gcd(P, N) only once every 128 iterations, accelerating Pollard's $\rho$ by $6\times$!
4. Complexity Matrix
Step through Extended GCD Bézout coefficients ax + by = gcd(a,b), Miller-Rabin modular squaring witness chains, and Pollard Rho cycle detection gcd(|x-y|, n) > 1.
Invariant & State
| Algorithm | Time Complexity | Space Complexity | 64-Bit ($N \le 2^{64}-1$) Guarantee |
|---|---|---|---|
| Trial Division | O(sqrt(N)) | O(1) | Up to $4 \times 10^9$ ops (Too slow for $10^{18}$) |
| Deterministic 64-bit Miller-Rabin (7 Bases) | O(7 * log^3 N) | O(1) | Exact 100% deterministic for all uint64_t |
| Pollard's $\rho$ (Brent + Batched GCD) | O(N^{1/4}) expected | O(1) | $\sim 10^4$ steps max for $N \le 10^{18}$ (< 2 ms) |
5. Production-Grade C++20 Implementation
#include <cstdint>
#include <vector>
#include <numeric>
#include <algorithm>
using u64 = uint64_t;
using u128 = __uint128_t;
inline u64 mul_mod(u64 a, u64 b, u64 m) {
return static_cast<u64>((static_cast<u128>(a) * b) % m);
}
u64 power_mod(u64 base, u64 exp, u64 mod) {
u64 res = 1;
base %= mod;
while (exp > 0) {
if (exp & 1) res = mul_mod(res, base, mod);
base = mul_mod(base, base, mod);
exp >>= 1;
}
return res;
}
// 100% Deterministic Miller-Rabin for all 64-bit unsigned integers
bool is_prime_u64(u64 n) {
if (n < 2) return false;
for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (n == p) return true;
if (n % p == 0) return false;
}
u64 d = n - 1;
int s = __builtin_ctzll(d);
d >>= s;
// Sinclair's 7 bases sufficient for all n < 2^64
for (u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
if (a % n == 0) continue;
u64 x = power_mod(a, d, n);
if (x == 1 || x == n - 1) continue;
bool composite = true;
for (int r = 1; r < s; ++r) {
x = mul_mod(x, x, n);
if (x == n - 1) {
composite = false;
break;
}
}
if (composite) return false;
}
return true;
}
// Pollard's Rho with Brent's cycle detection and batched GCD
u64 pollard_rho(u64 n) {
if (n % 2 == 0) return 2;
auto f = [n](u64 x, u64 c) { return (mul_mod(x, x, n) + c) % n; };
for (u64 c = 1; ; ++c) {
u64 x = 2, y = 2, g = 1, q = 1, xs = 2;
int m = 128;
for (int r = 1; g == 1; r <<= 1) {
x = y;
for (int i = 0; i < r; ++i) y = f(y, c);
for (int k = 0; k < r && g == 1; k += m) {
xs = y;
int lim = std::min(m, r - k);
for (int i = 0; i < lim; ++i) {
y = f(y, c);
u64 diff = x > y ? x - y : y - x;
q = mul_mod(q, diff, n);
}
g = std::gcd(q, n);
}
}
if (g == n) { // Overshot: backtrack step-by-step from saved xs
do {
xs = f(xs, c);
u64 diff = x > xs ? x - xs : xs - x;
g = std::gcd(diff, n);
} while (g == 1);
}
if (g < n) return g;
}
}
void factorize_rec(u64 n, std::vector<u64>& factors) {
if (n == 1) return;
if (is_prime_u64(n)) {
factors.push_back(n);
return;
}
u64 d = pollard_rho(n);
factorize_rec(d, factors);
factorize_rec(n / d, factors);
}Production Pitfall — Choosing $c = 0$ or $c = -2$ in Pollard's $\rho$ Polynomial
In Pollard's $\rho$, the iteration polynomial $f(x) = (x^2 + c) \bmod N$ must NEVER use $c = 0$ (which simply powers $x^{2^k}$ and fails to mix residues) or $c = -2$ (which maps $x = 2 \cos \theta$ to $2 \cos 2\theta$, trapping the walk on a tiny algebraic subset). Start at $c = 1$ and increment $c$ whenever a full cycle lands on the trivial divisor $g = N$.