Algorithmic Number Theory: Extended Euclidean Algorithm, Modular Inverse, CRT & Sieve of Eratosthenes

1. Formal Invariant — Bézout's Identity, CRT & Linear Sieve Uniqueness

2. Chinese Remainder Theorem Congruence Merging Pipeline

3. Linear-Time Modular Inverse Table in $O(N)$

When precomputing factorials and inverse factorials modulo a prime $M$ up to $N$, we can compute all modular inverses $\text{inv}[1..N]$ in $O(N)$ total time without a single log M binary exponentiation! Writing $M = q \cdot i + r$ where $q = \lfloor M / i \rfloor$ and $r = M \bmod i < i$, we have $q \cdot i + r \equiv 0 \pmod M$. Multiplying by $i^{-1} r^{-1}$ gives the one-line recurrence:

$$\text{inv}[i] = \left( M - \lfloor M / i \rfloor \right) \cdot \text{inv}[M \bmod i] \bmod M$$

4. Hardware & Cache Reality

5. Complexity Matrix

6. Production-Grade C++20 Implementation

#include <vector>
#include <cstdint>
#include <utility>

// Iterative Extended GCD: returns g = gcd(a, b) and sets x, y such that a*x + b*y = g
int64_t ext_gcd(int64_t a, int64_t b, int64_t& x, int64_t& y) {
    x = 1; y = 0;
    int64_t x1 = 0, y1 = 1, a1 = a, b1 = b;
    while (b1 != 0) {
        int64_t q = a1 / b1;
        std::tie(x, x1) = std::make_tuple(x1, x - q * x1);
        std::tie(y, y1) = std::make_tuple(y1, y - q * y1);
        std::tie(a1, b1) = std::make_tuple(b1, a1 - q * b1);
    }
    return a1;
}

// Modular inverse of a modulo m (works for ANY m where gcd(a, m) == 1)
int64_t mod_inv(int64_t a, int64_t m) {
    int64_t x, y;
    int64_t g = ext_gcd(a, m, x, y);
    if (g != 1) return -1; // Inverse does not exist
    return (x % m + m) % m;
}

// Generalized Chinese Remainder Theorem for potentially non-coprime moduli:
// Returns {r, lcm} such that X = r (mod lcm), or {-1, -1} if no solution exists
std::pair<int64_t, int64_t> crt_merge(int64_t r1, int64_t m1, int64_t r2, int64_t m2) {
    int64_t p, q;
    int64_t g = ext_gcd(m1, m2, p, q);
    if ((r2 - r1) % g != 0) return {-1, -1};

    int64_t step = m2 / g;
    __int128 diff = (r2 - r1) / g;
    int64_t k = static_cast<int64_t>((diff % step * (p % step)) % step);
    if (k < 0) k += step;

    int64_t lcm = m1 * step;
    int64_t r = static_cast<int64_t>((r1 + static_cast<__int128>(m1) * k) % lcm);
    return {r, lcm};
}

// Euler's O(N) Linear Sieve computing smallest prime factor spf[x] and totient phi[x]
struct LinearSieve {
    std::vector<int> spf;
    std::vector<int> phi;
    std::vector<int> primes;

    explicit LinearSieve(int n) : spf(n + 1, 0), phi(n + 1, 0) {
        phi[1] = 1;
        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                phi[i] = i - 1;
                primes.push_back(i);
            }
            for (int p : primes) {
                if (p > spf[i] || 1LL * i * p > n) break;
                spf[i * p] = p;
                if (i % p == 0) {
                    phi[i * p] = phi[i] * p; // p^2 divides i*p
                    break;
                } else {
                    phi[i * p] = phi[i] * (p - 1); // gcd(i, p) == 1
                }
            }
        }
    }
};