Miller-Rabin Primality Test

The Miller-Rabin algorithm is a probabilistic method used to determine if a large number is prime. It builds upon Fermat's Little Theorem, which states that for a prime $ p $ and integer $ a $ such that $ 1 \le a < p $, the congruence $ a^{p-1} \equiv 1 \pmod p $ holds. However, relying solely on this theorem allows for Carmichael numbers (composite numbers satisfying the congruence for all bases), necessitating a stronger verification method.

The algorithm utilizes the property that if $ p $ is prime, the only square roots of 1 modulo $ p $ are $ 1 $ and $ -1 $. Let $ n-1 $ be decomposed into $ d \cdot 2^s $, where $ d $ is odd. We compute $ x_0 = a^d \pmod n $. If $ x_0 \equiv 1 $ or $ x_0 \equiv -1 $, $ n $ is likely prime. Otherwise, we iteratively square $ x_i $ up to $ s-1 $ times. If any $ x_i \equiv -1 $, the test passes. If the sequence reaches 1 without passing through $ -1 $, or never reaches 1, $ n $ is composite.

bool miller_rabin_check(u64 n, u64 a, u64 d, int s) {
    u64 x = pow_mod(a, d, n);
    if (x == 1 || x == n - 1) return false;
    for (int r = 1; r < s; ++r) {
        x = mul_mod(x, x, n);
        if (x == n - 1) return false;
    }
    return true;
}

bool is_prime(u64 n) {
    if (n < 2) return false;
    if (n == 2 || n == 3) return true;
    if (n % 2 == 0) return false;
    
    u64 d = n - 1;
    int s = 0;
    while (d % 2 == 0) d >>= 1, s++;
    
    std::mt19937_64 rng(std::random_device{}());
    for (int i = 0; i < 10; ++i) {
        u64 a = rng() % (n - 3) + 2;
        if (miller_rabin_check(n, a, d, s)) return false;
    }
    return true;
}

Pollard-Rho Factorization

Pollard's Rho algorithm efficiently factors large composite integers, typically running in $ O(n^{1/4}) $ time complexity. It generates a pseudo-random sequence $ x_{i+1} = f(x_i) \pmod n $ where $ f(x) = x^2 + c $. The sequence eventually enters a cycle. By the birthday paradox, two distinct elements $ x_i $ and $ x_j $ in the sequence will eventually satisfy $ \gcd(|x_i - x_j|, n) > 1 $, yielding a non-trivial factor.

Optimization with Brent's Cycle Detection and Batch GCD

To improve efficiency, we combine Brent's cycle detection (doubling the step length to find the cycle) with a batched Greatest Common Divisor (GCD) computation. Calculating GCD is computationally expensive, so instead of checking every pair, we accumulate the product of differences $ |x_i - y| $ modulo $ n $ over a batch of iterations (e.g., 128 steps). A GCD check is performed on the accumulated product periodically. If $ \gcd(product, n) > 1 $, a factor is found.

u64 pollard_rho(u64 n) {
    if (n % 2 == 0) return 2;
    
    std::mt19937_64 rng(std::random_device{}());
    u64 x = rng() % (n - 2) + 1;
    u64 c = rng() % (n - 1) + 1;
    u64 y = x, d = 1;
    
    // Brent's algorithm with batched multiplication
    for (int len = 1; ; len <<= 1) {
        y = x;
        u64 product = 1;
        
        for (int i = 1; i <= len; ++i) {
            x = mul_mod(x, x, n); 
            x = (x + c) % n;
            u64 diff = (x > y) ? x - y : y - x;
            product = mul_mod(product, diff, n);
            
            // Check GCD every 128 steps to amortize cost
            if ((i & 127) == 0) {
                d = std::gcd(product, n);
                if (d > 1) return d;
            }
        }
        d = std::gcd(product, n);
        if (d > 1) return d;
    }
}

Complete Implementation for Maximum Prime Factor

The following template combines Miller-Rabin and Pollard-Rho to find the maximum prime factor of a given integer. It recursively factors the number $ n $: if $ n $ is prime, it updates the answer; otherwise, it finds a factor $ d $ using Pollard-Rho and recurses on both $ d $ and $ n/d $.

#include <bits/stdc++.h>
using namespace std;
using u64 = uint64_t;
using u128 = __uint128_t;

u128 mul_mod(u64 a, u64 b, u64 mod) { return (u128)a * b % mod; }

u64 pow_mod(u64 base, u64 exp, u64 mod) {
    u64 res = 1;
    while (exp) {
        if (exp & 1) res = mul_mod(res, base, mod);
        base = mul_mod(base, base, mod);
        exp >>= 1;
    }
    return res;
}

bool miller_rabin_check(u64 n, u64 a, u64 d, int s) {
    u64 x = pow_mod(a, d, n);
    if (x == 1 || x == n - 1) return false;
    for (int r = 1; r < s; ++r) {
        x = mul_mod(x, x, n);
        if (x == n - 1) return false;
    }
    return true;
}

bool is_prime(u64 n) {
    if (n < 2) return false;
    if (n == 2 || n == 3) return true;
    if (n % 2 == 0) return false;
    u64 d = n - 1; int s = 0;
    while (d % 2 == 0) d >>= 1, s++;
    std::mt19937_64 rng(random_device{}());
    for (int i = 0; i < 10; ++i) {
        u64 a = rng() % (n - 3) + 2;
        if (miller_rabin_check(n, a, d, s)) return false;
    }
    return true;
}

u64 pollard_rho(u64 n) {
    if (n % 2 == 0) return 2;
    std::mt19937_64 rng(random_device{}());
    u64 x = rng() % (n - 2) + 1, c = rng() % (n - 1) + 1;
    u64 y = x, d = 1;
    for (int len = 1; ; len <<= 1) {
        y = x; u64 product = 1;
        for (int i = 1; i <= len; ++i) {
            x = (mul_mod(x, x, n) + c) % n;
            u64 diff = (x > y) ? x - y : y - x;
            product = mul_mod(product, diff, n);
            if ((i & 127) == 0) {
                d = std::gcd(product, n);
                if (d > 1) return d;
            }
        }
        d = std::gcd(product, n);
        if (d > 1) return d;
    }
}

u64 max_prime_factor = 0;
void factorize(u64 n) {
    if (n <= max_prime_factor) return;
    if (is_prime(n)) {
        max_prime_factor = n;
        return;
    }
    u64 d = pollard_rho(n);
    while (n % d == 0) n /= d;
    factorize(d); factorize(n);
}

int main() {
    ios::sync_with_stdio(false); cin.tie(nullptr);
    int t; cin >> t;
    while (t--) {
        u64 n; cin >> n;
        max_prime_factor = 1;
        factorize(n);
        if (max_prime_factor == n) cout << "Prime\n";
        else cout << max_prime_factor << "\n";
    }
    return 0;
}