Number Theory Algorithms, in depth: gcd, extended Euclid, modular inverses, fast exponentiation, sieves, Euler's phi and the Chinese remainder theorem

Modular arithmetic in one page

Two integers are congruent modulo m, written a ≡ b (mod m), when m divides their difference. Addition, subtraction and multiplication respect congruence, so you can reduce after every operation and keep numbers small: (a * b) mod m equals ((a mod m) * (b mod m)) mod m. Division does not. You cannot divide by 3 modulo 6, because 3 has no partner that multiplies with it to give 1. Whether such a partner, the modular inverse, exists is the question the next two tools answer.

One language detail matters from the start. In Python, -7 % 3 is 2. In C, C++, Java and Go it is -1, because % follows the sign of the dividend. Code that subtracts and then reduces must normalise with ((x % m) + m) % m or Java's Math.floorMod, or negative residues will leak into array indices and comparisons.

Euclid's algorithm and why it is fast

The greatest common divisor of a and b is unchanged if you replace the larger by its remainder modulo the smaller, because any number dividing both also divides the remainder. Repeat until the remainder is zero; the last non-zero value is the gcd. For 240 and 46 the remainders are 10, 6, 4, 2 and then 0, so the gcd is 2 after five division steps.

Each pair of steps at least halves the larger number, so the step count is O(log min(a, b)). Lamé's theorem sharpens this: the worst case is consecutive Fibonacci numbers, and the number of steps never exceeds five times the number of decimal digits of the smaller input. For 64-bit integers that is under a hundred iterations. Complexity notation is explained in Big-O analysis.

Extended Euclid and the modular inverse

Bezout's identity says that for any a and b there are integers x and y with a·x + b·y = gcd(a, b). The extended algorithm finds them by carrying the coefficients along with the remainders. For 240 and 46 it returns (2, -9, 47), and indeed 240 × -9 + 46 × 47 = -2160 + 2162 = 2.

If gcd(a, m) = 1, the identity reads a·x + m·y = 1, and reducing modulo m gives a·x ≡ 1. So x mod m is the inverse of a. If the gcd is larger than 1, no inverse exists, and code must say so instead of returning garbage. This one fact explains why moduli in hashing and combinatorics are chosen prime: modulo a prime, every non-zero value has an inverse.

def ext_gcd(a, b):
    """Return (g, x, y) with a*x + b*y == g == gcd(a, b)."""
    old_r, r = a, b
    old_x, x = 1, 0
    old_y, y = 0, 1
    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_x, x = x, old_x - q * x
        old_y, y = y, old_y - q * y
    return old_r, old_x, old_y

def mod_inverse(a, m):
    g, x, _ = ext_gcd(a % m, m)
    if g != 1:
        raise ValueError(f"{a} has no inverse mod {m} (gcd={g})")
    return x % m

The iterative form matters with big integers. A recursive extended Euclid is shorter and fine for 64-bit inputs, which need fewer than a hundred levels, but on cryptographic-size inputs of thousands of bits the step count reaches four figures, past Python's default recursion limit of 1,000. Python 3.8 and later also expose the inverse directly as pow(a, -1, m), which raises ValueError when none exists.

Fast modular exponentiation

Computing b^e mod m by multiplying e times is hopeless when e has hundreds of digits. Square-and-multiply reads the exponent in binary: it keeps b^(2^k) by repeated squaring and multiplies it into the result whenever bit k of the exponent is set. That is at most 2·log2(e) multiplications, each reduced modulo m so the numbers never grow.

def pow_mod(base, exp, m):
    if m == 1:
        return 0
    result, base = 1, base % m
    while exp > 0:
        if exp & 1:                 # this bit of the exponent is set
            result = result * base % m
        base = base * base % m      # base^(2^k) for the next bit
        exp >>= 1
    return result

The overflow trap is in result * base. In Python integers are unbounded, so the code above is correct for any size. In a language with 64-bit integers, the product of two residues below m is below m², which fits in a signed 64-bit value only when m is below about 3.04 × 10^9. The popular modulus 1,000,000,007 is safe. A modulus near 2^61, common in string hashing, is not: use unsigned __int128 in GCC or Clang, BigInteger.modPow in Java, or Montgomery multiplication when speed matters.

Fermat's little theorem gives a second way to invert. If p is prime and a is not a multiple of p, then a^(p-1) ≡ 1 (mod p), so a^(p-2) is the inverse. For example, 3^5 mod 7 is 5, and 3 × 5 = 15 ≡ 1 (mod 7). This only works for a prime modulus; for a composite modulus use extended Euclid.

Sieves: every prime up to n

To find all primes up to n, the sieve of Eratosthenes takes each unmarked number p in turn and marks its multiples, starting at p² because smaller multiples were already marked by smaller primes. It runs in O(n log log n) time and n bits or bytes of memory. There are 25 primes up to 100 and 168 up to 1,000.

The linear sieve does a little more work per step and gains something valuable. It records the smallest prime factor of every number, marking each composite exactly once, by its smallest prime, so the total work is O(n). With that table, factorising any x ≤ n is a loop of divisions by spf[x], taking O(log x) steps. That turns many queries such as "count the divisors of each of a million numbers" from square-root time each into near-constant time.

def spf_sieve(n):
    """Linear sieve: smallest prime factor of every 2..n, plus the primes."""
    spf = [0] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if spf[i] == 0:             # nobody marked i, so i is prime
            spf[i] = i
            primes.append(i)
        for p in primes:            # mark i*p once, by its smallest prime p
            if p > spf[i] or i * p > n:
                break
            spf[i * p] = p
    return spf, primes

def factorize(x, spf):
    factors = {}
    while x > 1:
        p = spf[x]
        factors[p] = factors.get(p, 0) + 1
        x //= p
    return factors                  # factorize(360, spf) == {2: 3, 3: 2, 5: 1}

Memory is the limit. A sieve to 10^9 needs about a gigabyte as bytes, and an SPF table of 32-bit integers needs four times that. For a range [L, R] far from zero, use a segmented sieve: generate primes up to √R, then cross off their multiples inside one cache-sized window of the range at a time.

Euler's totient

Euler's function φ(n) counts the integers from 1 to n that are coprime to n. For a prime p it is p - 1. In general, for each distinct prime factor p of n, multiply n by (1 - 1/p): 36 = 2² × 3², so φ(36) = 36 × 1/2 × 2/3 = 12. Compute it by trial division up to √n, or from the SPF table when you need it for many values.

Euler's theorem generalises Fermat: if gcd(a, n) = 1 then a^φ(n) ≡ 1 (mod n). Its practical use is exponent reduction. To compute a power tower or a huge exponent modulo n, reduce the exponent modulo φ(n) first, but only when the base is coprime to n.

The Chinese remainder theorem

Given x ≡ r1 (mod m1) and x ≡ r2 (mod m2), the theorem says there is exactly one solution modulo lcm(m1, m2) if r1 ≡ r2 (mod gcd(m1, m2)), and none otherwise. Most textbook code assumes coprime moduli and silently returns a wrong answer when they are not. The version below handles the general case and reports contradictions.

def crt(r1, m1, r2, m2):
    """Solve x = r1 (mod m1), x = r2 (mod m2). Return (x, lcm) or None."""
    g, p, _ = ext_gcd(m1, m2)
    if (r2 - r1) % g:
        return None                 # the congruences contradict each other
    lcm = m1 // g * m2
    k = (r2 - r1) // g * p % (m2 // g)
    return (r1 + m1 * k) % lcm, lcm

Fold it over a list to combine many congruences. With x ≡ 2 (mod 3), x ≡ 3 (mod 5) and x ≡ 2 (mod 7), the first two give 8 modulo 15, and adding the third gives 23 modulo 105. With non-coprime moduli, x ≡ 3 (mod 4) and x ≡ 5 (mod 6) give 11 modulo 12, while x ≡ 3 (mod 4) and x ≡ 4 (mod 6) have no solution, since one demands an odd number and the other an even one.

In practice the theorem appears in scheduling, where jobs with periods 4 and 6 and offsets 3 and 5 first coincide at time 11 and every 12 ticks after, in exact big-integer arithmetic done as several small modular computations, and in speeding up RSA decryption by working modulo p and q separately.

Worked example: a toy RSA key

Take primes p = 61 and q = 53. The modulus is n = 3233 and φ(n) = 60 × 52 = 3120. Choose the public exponent e = 17; Euclid confirms gcd(17, 3120) = 1, so an inverse exists, and extended Euclid gives the private exponent d = 2753. To encrypt the message 65, fast exponentiation computes 65^17 mod 3233 = 2790. To decrypt, 2790^2753 mod 3233 returns 65. The round trip works because e·d ≡ 1 (mod φ(n)), and Euler's theorem then collapses m^(e·d) back to m.

This is a teaching example only. Real RSA uses moduli of 2,048 bits or more, randomised padding such as OAEP, and constant-time implementations that do not leak the exponent through timing. Never deploy textbook RSA or your own modular code for cryptography; use a vetted library.

Primality and factoring at scale

Trial division to √n is fine below about 10^12. Beyond that, use the Miller-Rabin test, which needs only fast exponentiation and, with a known set of bases, is deterministic for all 64-bit integers. To factor numbers too large for a sieve, Pollard's rho finds a factor in roughly n^(1/4) steps by looking for a cycle in a pseudo-random sequence modulo the hidden factor, using the same tortoise-and-hare idea as Floyd's cycle detection.

Failure modes

• Silent overflow. Multiplying two residues of a large modulus in 64-bit arithmetic wraps around and gives plausible but wrong answers. Test with moduli near your limit.
• Negative remainders. % in C-family languages returns negative values for negative inputs.
• Inverting a non-unit. Asking for the inverse of 4 modulo 6 must fail loudly, not return 0.
• Fermat with a composite modulus. a^(m-2) is not an inverse unless m is prime.
• Coprime-only CRT. Code that ignores the gcd gives wrong results when moduli share factors.
• Sieve memory. A full sieve to 10^9 exhausts memory on modest machines; segment it.
• Weak hashing moduli. A composite or small modulus in a polynomial hash increases collisions; see hash table design.

What to do next

1. Implement iterative extended Euclid and a modular inverse that raises on non-coprime input; check them against brute force for random small values.
2. Write square-and-multiply in your main language and test it with a modulus near the overflow limit of that language.
3. Build a linear SPF sieve to 10^6 and use it to count divisors of every number in the range.
4. Implement general CRT and test both a solvable and a contradictory non-coprime pair.
5. Reproduce the toy RSA example by hand, then read the Miller-Rabin page and add a deterministic 64-bit primality test.
6. Audit existing code for % on negative values and for products that can exceed 63 bits.