4.28.9 — GCD, Primes & Modular Arithmetic

Three sentences turn up in problem statements again and again: "return the greatest common divisor", "count the primes less than n", and, most often of all, "the answer may be very large, so return it modulo 10⁹ + 7". Each has one standard technique behind it, and each has a trap that costs people the problem even when they know the technique. This page builds all three from the ground up, in Python and TypeScript, with the traps measured.

1. The greatest common divisor, and Euclid's shortcut

The greatest common divisor of two whole numbers, written \gcd(a, b), is the largest number that divides both exactly. \gcd(12, 18) = 6. It is what you need to reduce a fraction, to find when two repeating schedules line up again, or to solve "Greatest Common Divisor of Strings" (LeetCode 1071), where the answer is a prefix whose length is the gcd of the two lengths.

The slow way tries every candidate from \min(a, b) downward until one divides both. For two numbers around 10^{12} that can be a trillion divisions.

Euclid's rule (about 300 BC, and still the method computers use): \gcd(a, b) = \gcd(b, a \bmod b), where a \bmod b ("a mod b") is the remainder after dividing a by b. Keep applying it until the second number is 0; the first number is then the answer.

Why the rule is true. Write a = q \cdot b + r, where q is how many whole times b fits into a and r is the remainder. Any number that divides both a and b also divides a - q \cdot b, which is r. And any number that divides both b and r also divides q \cdot b + r, which is a. So the pairs (a, b) and (b, r) have exactly the same common divisors, and therefore the same greatest one. The second pair is smaller, which is the whole point.

Traced on \gcd(1071, 462):

a	b	a = q \cdot b + r	next pair
1071	462	1071 = 2 \cdot 462 + 147	(462, 147)
462	147	462 = 3 \cdot 147 + 21	(147, 21)
147	21	147 = 7 \cdot 21 + 0	(21, 0)

The second number reached 0, so the answer is 21. Check it: 1071 = 21 \times 51 and 462 = 21 \times 22, and 51 and 22 share no factor.

How fast, derived. Look at any two consecutive steps. If b \le a/2, the remainder a \bmod b is smaller than b, so it is below a/2. If b > a/2, then b fits into a exactly once and the remainder is a - b, which is again below a/2. Either way, every two steps at least halve the larger number, so the number of steps is at most about 2 \log_2 \min(a, b). For numbers near 10^{12} that is at most about 80 steps, instead of a trillion.

The least common multiple, the smallest number both divide, follows directly: \text{lcm}(a, b) = a / \gcd(a, b) \times b. Divide before multiplying, so the intermediate value never exceeds the answer.

2. Counting primes: the sieve

A prime is a whole number greater than 1 whose only divisors are 1 and itself. "Count the primes less than n" (LeetCode 204) allows n up to 5 \times 10^6.

Testing one number only needs divisors up to its square root. If x = d \times e with both factors above \sqrt{x}, their product would be above x, which is impossible; so any composite x has a factor no larger than \sqrt{x}. Testing each number up to 5 \times 10^6 this way costs up to about 2,200 divisions each, around 10^{10} operations in total. Far too slow, and wasteful: the same small primes are tried against every number.

The sieve of Eratosthenes turns the work around. Instead of asking of each number "who divides you?", it takes each prime and crosses off everything it divides. Write down 2, 3, \dots, n - 1. The first number not crossed off, 2, is prime; cross off its multiples. The next survivor, 3, is prime; cross off its multiples. Continue. Whatever survives is prime.

Two details make it fast:

• Start crossing at p^2, not 2p. A smaller multiple k \cdot p with k < p has a prime factor smaller than p, so an earlier prime already crossed it off.
• Stop once p^2 \ge n. By the square-root argument above, every composite below n has a prime factor below \sqrt{n}, so it has already been crossed off.

For n = 30: prime 2 crosses 4, 6, 8, …, 28; prime 3 crosses 9, 15, 21, 27 (12, 18 and 24 were already gone); 4 is crossed, so it is skipped; prime 5 crosses 25; and 6^2 = 36 is past 30, so the sieve stops. The survivors are 2, 3, 5, 7, 11, 13, 17, 19, 23, 29: 10 primes.

How fast. Prime p crosses off about n/p numbers, so the total is n \times (\tfrac12 + \tfrac13 + \tfrac15 + \tfrac17 + \dots), the sum running over primes up to \sqrt{n}. A result of Franz Mertens (1874) says that sum of 1/p grows like \ln \ln n, the logarithm of the logarithm, which is almost constant. So the sieve does about n \ln \ln n steps. Measured for n = 5 \times 10^6: it crosses off 11.2 million entries, about 2.2 per number, and the Python version below counts the primes below 5 \times 10^6 (there are 348,513) in about 0.03 seconds. Memory is one byte per number, 5 MB.

3. Modular arithmetic, and why 10⁹ + 7

"How many ways are there to climb 10,000 stairs?" has an answer with more than 2,000 digits. Problems that would produce such numbers ask for the answer modulo a fixed number m, meaning the remainder after dividing by m, so every intermediate value stays small. The rules that make this possible:

(a + b) \bmod m = \big((a \bmod m) + (b \bmod m)\big) \bmod m

(a \times b) \bmod m = \big((a \bmod m) \times (b \bmod m)\big) \bmod m

Read aloud: you may take the remainder after every addition and every multiplication, as often as you like, and the final remainder comes out the same. So a dynamic program that adds counts can apply % MOD after each addition and never hold a big number.

Why 10⁹ + 7 in particular? Three properties. It is prime, which makes division possible (below). It is just under 2^{30}, so the sum of two remainders, below 2 \times 10^9 + 14, still fits in a 32-bit signed integer (largest value about 2.147 \times 10^9). And the product of two remainders is below (10^9 + 7)^2 \approx 10^{18}, which fits in a 64-bit signed integer (largest about 9.2 \times 10^{18}). It is large enough that answers rarely collide, and small enough that the arithmetic never overflows in C++ or Java.

Subtraction needs care, because a remainder can go negative: (a - b) \bmod m should be computed as ((a - b) % m + m) % m. And here the two languages genuinely disagree. In Python, -7 % 3 is 2: the result takes the sign of the divisor, so it is always non-negative for a positive modulus. In JavaScript and TypeScript, -7 % 3 is -1: the result takes the sign of the dividend. A JavaScript solution that subtracts without adding m back returns negative answers and fails.

Multiplication is where JavaScript silently breaks. A JavaScript number is a 64-bit floating-point value, exact only for whole numbers up to 2^{53} \approx 9.0 \times 10^{15}. The product of two remainders can reach 10^{18}, far beyond that, so the low digits are lost before % even runs. Measured: 1_000_000_006 * 1_000_000_006 % 1_000_000_007 prints 999999972 in Node.js. The true answer is 1 (since 10^9 + 6 \equiv -1, and (-1)^2 = 1). The fix is BigInt, JavaScript's arbitrary-size integers, written with a trailing n; with BigInt the same expression gives 1n. Python's integers never overflow, so Python never has this problem.

Division is the one operation that does not carry over directly. You cannot divide remainders: 10 / 2 = 5, but working modulo 7, 10 \bmod 7 = 3, and 3 / 2 is not a whole number. Instead, multiply by a modular inverse: a number x with a \cdot x \bmod m = 1, which plays the part of 1/a. When m is prime and a is not a multiple of m, Fermat's little theorem says a^{m-1} \bmod m = 1, so a \cdot a^{m-2} \bmod m = 1, and the inverse is a^{m-2} \bmod m. That power is computed in O(\log m) multiplications with the fast-power loop of 4.28.6.

Check it on a small prime by hand: the inverse of 3 modulo 7 is 3^5 \bmod 7 = 243 \bmod 7 = 5 (because 7 \times 34 = 238), and indeed 3 \times 5 = 15 = 2 \times 7 + 1, remainder 1.

The classic use is the binomial coefficient \binom{n}{k} = \frac{n!}{k!\,(n-k)!}, "n choose k", the number of ways to pick k items from n. Modulo a prime it becomes n! \times \text{inverse}(k!) \times \text{inverse}((n-k)!), all modulo m.

The code

Python

MOD = 10**9 + 7


def gcd(a: int, b: int) -> int:
    while b:                          # (1)
        a, b = b, a % b               # (2)
    return a


def lcm(a: int, b: int) -> int:
    return a // gcd(a, b) * b         # (3) divide first, then multiply


def count_primes_below(n: int) -> int:
    if n < 3:
        return 0
    is_prime = bytearray([1]) * n     # (4) one byte per number
    is_prime[0] = is_prime[1] = 0
    p = 2
    while p * p < n:                  # (5)
        if is_prime[p]:
            is_prime[p * p::p] = bytes(len(range(p * p, n, p)))   # (6) cross off p², p²+p, ...
        p += 1
    return sum(is_prime)


def mod_pow(base: int, exp: int, m: int = MOD) -> int:
    result, base = 1, base % m
    while exp:                        # (7) the loop from 4.28.6, with % m after each multiply
        if exp & 1:
            result = result * base % m
        base = base * base % m
        exp >>= 1
    return result


def mod_inverse(a: int, m: int = MOD) -> int:
    return mod_pow(a, m - 2, m)       # (8) Fermat: valid only when m is prime and a % m != 0


def n_choose_k(n: int, k: int, m: int = MOD) -> int:
    if k < 0 or k > n:
        return 0
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i % m                         # (9)
    return fact[n] * mod_inverse(fact[k], m) % m * mod_inverse(fact[n - k], m) % m   # (10)


if __name__ == "__main__":
    print(gcd(1071, 462), lcm(4, 6))                  # → 21 12
    print(count_primes_below(30))                     # → 10
    print(count_primes_below(5_000_000))              # → 348513
    print(mod_inverse(3, 7), 3 * mod_inverse(3, 7) % 7)   # → 5 1
    print(n_choose_k(10, 3), n_choose_k(100_000, 50_000))
    print(-7 % 3)                                     # → 2 in Python

TypeScript

const MOD = 1_000_000_007n;

function gcd(a: number, b: number): number {
  while (b !== 0) [a, b] = [b, a % b];              // (1) (2)
  return a;
}

const lcm = (a: number, b: number): number => (a / gcd(a, b)) * b;   // (3)

function countPrimesBelow(n: number): number {
  if (n < 3) return 0;
  const isPrime = new Uint8Array(n).fill(1);         // (4)
  isPrime[0] = isPrime[1] = 0;
  for (let p = 2; p * p < n; p++) {                  // (5)
    if (!isPrime[p]) continue;
    for (let q = p * p; q < n; q += p) isPrime[q] = 0;   // (6)
  }
  let count = 0;
  for (const flag of isPrime) count += flag;
  return count;
}

function modPow(base: bigint, exp: bigint, m: bigint = MOD): bigint {
  let result = 1n;
  base %= m;
  while (exp > 0n) {                                 // (7)
    if (exp & 1n) result = (result * base) % m;
    base = (base * base) % m;
    exp >>= 1n;
  }
  return result;
}

const modInverse = (a: bigint, m: bigint = MOD): bigint => modPow(a, m - 2n, m);   // (8)

function nChooseK(n: number, k: number, m: bigint = MOD): bigint {
  if (k < 0 || k > n) return 0n;
  const fact: bigint[] = [1n];
  for (let i = 1; i <= n; i++) fact.push((fact[i - 1] * BigInt(i)) % m);     // (9)
  return (fact[n] * modInverse(fact[k], m) % m) * modInverse(fact[n - k], m) % m;   // (10)
}

console.log(gcd(1071, 462), lcm(4, 6));            // → 21 12
console.log(countPrimesBelow(30), countPrimesBelow(5_000_000));   // → 10 348513
console.log(modInverse(3n, 7n), nChooseK(10, 3), nChooseK(100_000, 50_000));   // → 5n 120n 149033233n
console.log(-7 % 3);                                // → -1 in JavaScript
console.log(1_000_000_006 * 1_000_000_006 % 1_000_000_007);   // → wrong: Number loses precision
console.log((1_000_000_006n * 1_000_000_006n) % 1_000_000_007n);  // → 1n, the right answer

1. Euclid's loop runs until the second number is 0.
2. One step of the rule: the pair (a, b) becomes (b, a \bmod b). Both languages swap with a single tuple or array assignment so neither value is overwritten before it is read.
3. Dividing first keeps the intermediate value no larger than the final answer.
4. The sieve's table: one byte per number, all assumed prime to start.
5. Stop once p^2 reaches n, by the square-root argument.
6. Cross off multiples starting at p^2. The Python version assigns a whole slice of zero bytes at once, which runs in C and is much faster than a Python loop; the TypeScript version loops, which is already fast in JavaScript.
7. The fast-power loop from 4.28.6, with % m after every multiplication so no value exceeds m^2. In TypeScript every value is a BigInt (1n, exp & 1n), because the products would otherwise lose precision.
8. Fermat's inverse. The comment is part of the code: if m is not prime, this silently returns a wrong answer.
9. Factorials modulo m, each from the previous one.
10. n! times the inverses of k! and (n-k)!, taking the remainder after each multiplication.

Both versions print 21 12, then 10 and 348513, then the inverse 5, \binom{10}{3} = 120 and \binom{100000}{50000} \bmod (10^9+7) = 149033233. That last value was checked against Python's exact math.comb (computed with the full 30,101-digit number, then reduced). The gcd was checked against math.gcd on 2,000 random pairs, and the sieve against trial division for every n below 400.

Complexity

• gcd: O(\log \min(a, b)) steps, from the halving argument in section 1; O(1) space.
• Sieve: O(n \log \log n) time and O(n) space, from the sum over primes in section 2.
• Modular power and inverse: O(\log e) multiplications for exponent e, so O(\log m) for an inverse.
• \binom{n}{k} \bmod m as written: O(n) to build the factorials plus two inverses, O(n + \log m) in total.

Traps, collected

• JavaScript % keeps the sign of the dividend, so add m back after any subtraction.
• JavaScript numbers lose exactness above 2^{53}, so multiply remainders as BigInt. Measured above: 999999972 instead of 1.
• Fermat's inverse needs a prime modulus. For a non-prime modulus, use the extended Euclidean algorithm, and the inverse exists only when \gcd(a, m) = 1.
• Apply the modulus at every step, not only at the end. In Python the result is still correct, but it is slow, because the intermediate numbers grow to thousands of digits; in other languages the intermediate overflows.
• "Count primes less than n" means strictly less. The table has n entries, indices 0 to n-1.

Variations you will meet

• Many binomial queries: precompute inverse factorials once, from the top down: compute the inverse of n! with one power, then \text{inverse}((i-1)!) = \text{inverse}(i!) \times i. Every query is then O(1).
• Prime factorisation of many numbers: a sieve that stores each number's smallest prime factor lets any number up to n be factorised in O(\log n) steps.
• Fast power with a modulus is the core operation of RSA and Diffie-Hellman key exchange (8.2).
• Rolling hashes for string matching use the same modular arithmetic with a large prime (4.31).

What the interviewer will push on

"Why does Euclid's algorithm work, and how fast is it?" They want the common-divisor argument (r = a - qb), not just the code, and the halving argument for the step count. "It's logarithmic" with no reason is the memorised answer.

"Why start the sieve at p²?" Because every smaller multiple of p has a smaller prime factor and is already crossed off. Saying this, and why the outer loop stops at \sqrt{n}, shows you understand the sieve rather than having copied it.

"Why 10⁹ + 7?" Prime (so division by inverses works), sums fit in 32 bits, products fit in 64 bits. Three reasons in one breath is the strong answer.

"Your JavaScript solution fails on large inputs but the logic is right. Why?" Almost always one of the two traps: a product above 2^{53}, or a negative remainder after subtraction. Naming both, with the fix for each, is exactly what is being checked.

"How do you divide under a modulus?" Multiply by the modular inverse, via Fermat when the modulus is prime, and say the condition out loud: it fails for a non-prime modulus or when a is a multiple of it.

One thing to volunteer: the modulus is not only a convenience for judges. Every hash table bucket index, every rolling hash and every public-key signature is arithmetic modulo something, so these three rules are the same ones running inside a database index and a TLS handshake.

Recall

• Euclid: \gcd(a, b) = \gcd(b, a \bmod b), because (a, b) and (b, r) share the same divisors; at most about 2 \log_2 \min(a, b) steps since every two steps halve the larger number. \text{lcm} = a / \gcd \times b.
• Sieve of Eratosthenes: cross off multiples of each prime from p^2, stop when p^2 \ge n; O(n \log \log n) time, O(n) bytes.
• Work modulo m by reducing after every + and \times. 10^9 + 7: prime, sums fit 32 bits, products fit 64 bits.
• JavaScript traps: % can return negatives (add m back), and products above 2^{53} lose digits (use BigInt).
• Division means multiplying by a modular inverse: a^{m-2} \bmod m for prime m (Fermat), computed with fast power.

Self-test: Trace \gcd(252, 105) · Why does the sieve start at p^2? · Compute the inverse of 4 modulo 11 by hand and check it · Why does a * b % MOD go wrong in JavaScript but not in Python? · Compute \binom{6}{2} \bmod 7 using inverses.

Next: 4.29 drops to the level of individual bits — the idioms worth memorising, the probabilistic structures Parts 10 and 11 lean on, and what to do when a problem turns out to be NP-hard.