Number Theory for CS

The foundation: what does it mean for one integer to divide another?

We say $a \mid b$ (read: a divides b) if there exists an integer $k$ such that $b = ak$. For example, $3 \mid 12$ because $12 = 3 \cdot 4$, but $5 \nmid 12$.

The Division Algorithm (actually a theorem, not an algorithm in the modern sense) states: for any integers $a$ and $d > 0$, there exist unique integers $q$ (quotient) and $r$ (remainder) such that: $$a = qd + r \quad \text{with } 0 \le r < d$$

This uniqueness is what makes modular arithmetic well-defined. The remainder $r$ is what Python computes as a % d. In C/C++, % can return a negative value when a < 0, which is a common source of bugs — always normalize if you need the mathematical remainder.

Key divisibility properties (provable directly from the definition):

• If $a \mid b$ and $a \mid c$, then $a \mid (b + c)$ and $a \mid (b - c)$
• If $a \mid b$ and $b \mid c$, then $a \mid c$ (transitivity)
• If $a \mid b$ and $a \mid c$, then $a \mid (mb + nc)$ for all integers $m, n$ — this linear combination closure is crucial for the GCD proof

These properties might seem trivial but they carry enormous weight: the entire theory of GCDs, primes, and modular arithmetic is built on them.

# Division algorithm in Python
def div_algorithm(a: int, d: int) -> tuple[int, int]:
    """Returns (q, r) such that a = q*d + r, 0 <= r < d."""
    assert d > 0
    q = a // d  # floor division — always rounds toward -inf
    r = a % d   # always non-negative in Python
    assert a == q * d + r and 0 <= r < d
    return q, r

print(div_algorithm(17, 5))   # (3, 2): 17 = 3*5 + 2
print(div_algorithm(-7, 3))   # (-3, 2): -7 = -3*3 + 2

Thinking on a clock: the integers mod n

Fix a positive integer $n$ called the modulus. We write $a \equiv b \pmod{n}$ (read: a is congruent to b mod n) when $n \mid (a - b)$, i.e., $a$ and $b$ have the same remainder when divided by $n$.

The key insight is that congruence mod $n$ is an equivalence relation: reflexive ($a \equiv a$), symmetric, and transitive. This partitions the integers into exactly $n$ residue classes: ${0, 1, 2, \ldots, n-1}$. Every integer belongs to exactly one class.

Arithmetic works inside residue classes — this is what makes the whole edifice useful:

• If $a \equiv a' \pmod{n}$ and $b \equiv b' \pmod{n}$, then:
  • $a + b \equiv a' + b' \pmod{n}$
  • $a \cdot b \equiv a' \cdot b' \pmod{n}$

Proof sketch: $a = a' + kn$ and $b = b' + ln$, so $ab = a'b' + n(a'l + b'k + kln)$, which is divisible by $n$.

This compatibility means we can reduce at any point in a computation. Want to compute $2^{100} \bmod 13$? You do not need to compute the 30-digit number $2^{100}$ first — reduce after every multiplication.

The ring $\mathbb{Z}/n\mathbb{Z}$: the set ${0, 1, \ldots, n-1}$ with addition and multiplication mod $n$ forms a ring. When $n$ is prime it is actually a field — every nonzero element has a multiplicative inverse. This algebraic structure is why prime moduli appear everywhere in cryptography.

def mod_exp(base: int, exp: int, mod: int) -> int:
    """Compute base^exp % mod using repeated squaring.
    Time: O(log exp) multiplications — essential for crypto."""
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:          # odd exponent: multiply in current base
            result = result * base % mod
        base = base * base % mod  # square
        exp >>= 1            # halve exponent
    return result

# 2^100 mod 13 without ever computing 2^100
print(mod_exp(2, 100, 13))   # 3

# Verify: 2^12 ≡ 1 (mod 13) by Fermat, 100 = 8*12 + 4, so 2^100 ≡ 2^4 = 16 ≡ 3
print(16 % 13)               # 3 ✓

The oldest nontrivial algorithm in mathematics

The greatest common divisor $\gcd(a, b)$ is the largest positive integer dividing both $a$ and $b$. A crucial equivalent characterization: $\gcd(a, b)$ is the smallest positive integer expressible as a linear combination $ma + nb$ for integers $m, n$. This is Bezout's theorem, and it is constructive — the extended Euclidean algorithm finds those $m, n$.

The Euclidean algorithm exploits one key observation: $$\gcd(a, b) = \gcd(b, a \bmod b)$$

Proof: if $a = qb + r$, then any common divisor of $a$ and $b$ also divides $r = a - qb$, and vice versa. So the set of common divisors is unchanged, and hence so is the GCD.

Apply this repeatedly until the remainder is 0: $$\gcd(48, 18) = \gcd(18, 12) = \gcd(12, 6) = \gcd(6, 0) = 6$$

Complexity: the number of steps is $O(\log \min(a,b))$ — each two steps reduce the problem size by at least half (Lame's theorem, using Fibonacci numbers for the worst case). For 2048-bit RSA moduli, this is about 2048 steps — completely practical.

The Extended Euclidean Algorithm additionally finds integers $x, y$ such that: $$ax + by = \gcd(a, b)$$

This is essential for computing modular inverses: if $\gcd(a, n) = 1$, then $ax \equiv 1 \pmod{n}$, so $x = a^{-1} \bmod n$.

Key corollary: $\gcd(a, n) = 1$ is necessary and sufficient for the inverse to exist. When working mod a prime $p$, every nonzero element has an inverse because $\gcd(a, p) = 1$ for all $1 \le a < p$.

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

def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Returns (g, x, y) with a*x + b*y == g == gcd(a, b)."""
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extended_gcd(b, a % b)
    # Back-substitute: a*x + b*y = g
    # We know: b*x1 + (a%b)*y1 = g
    #          b*x1 + (a - (a//b)*b)*y1 = g
    #          a*y1 + b*(x1 - (a//b)*y1) = g
    return g, y1, x1 - (a // b) * y1

def mod_inverse(a: int, n: int) -> int:
    """Compute a^{-1} mod n. Raises if gcd(a,n) != 1."""
    g, x, _ = extended_gcd(a % n, n)
    if g != 1:
        raise ValueError(f"{a} has no inverse mod {n}")
    return x % n

print(gcd(48, 18))             # 6
print(extended_gcd(35, 15))    # (5, 1, -2): 35*1 + 15*(-2) = 5
print(mod_inverse(3, 11))      # 4: 3*4 = 12 ≡ 1 (mod 11)

The atoms of arithmetic

An integer $p > 1$ is prime if its only positive divisors are 1 and itself. The Fundamental Theorem of Arithmetic states: every integer $n > 1$ has a unique factorization into primes (up to ordering): $$n = p_1^{e_1} p_2^{e_2} \cdots p_k^{e_k}$$

This is a non-trivial theorem — existence is proved by strong induction, uniqueness requires Euclid's lemma: if $p$ is prime and $p \mid ab$, then $p \mid a$ or $p \mid b$. This lemma fails for composite numbers (e.g., $4 \mid 2 \cdot 6$ but $4 \nmid 2$ and $4 \nmid 6$), which is why primes are special.

How many primes are there? Infinitely many — Euclid's proof is a marvel of simplicity. Suppose $p_1, \ldots, p_k$ were all primes; let $N = p_1 p_2 \cdots p_k + 1$. Then $N$ is either prime (contradicting completeness) or has a prime factor not in our list (since $p_i \nmid N$ for all $i$). Either way the list was incomplete.

The Prime Number Theorem (deep result): the number of primes up to $n$ is asymptotically $n / \ln n$. For 2048-bit numbers, primes are dense enough that you find one in about 1400 random trials — making probabilistic primality testing practical for RSA key generation.

Trial division for testing primality: check all primes up to $\sqrt{n}$. If $n$ has no factor $\le \sqrt{n}$, it is prime (since if $n = ab$ with $a \le b$, then $a \le \sqrt{n}$). This is $O(\sqrt{n})$ — too slow for large crypto numbers.

Miller-Rabin primality test: a probabilistic, polynomial-time test based on Fermat's Little Theorem. For each witness $a$, if $n$ is prime, a certain condition must hold; composites fail for most witnesses. In practice you pick several random witnesses and declare primality if all pass.

def sieve_of_eratosthenes(n: int) -> list[int]:
    """All primes up to n. O(n log log n) time, O(n) space."""
    is_prime = bytearray([1]) * (n + 1)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            # Mark all multiples of i starting at i^2
            is_prime[i*i::i] = bytearray(len(is_prime[i*i::i]))
    return [i for i, v in enumerate(is_prime) if v]

primes = sieve_of_eratosthenes(50)
print(primes)
# [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

def is_prime_trial(n: int) -> bool:
    if n < 2: return False
    if n < 4: return True
    if n % 2 == 0 or n % 3 == 0: return False
    i = 5
    while i * i <= n:
        if n % i == 0 or n % (i + 2) == 0:
            return False
        i += 6  # primes > 3 are ±1 mod 6
    return True

The theorem that makes RSA work

Theorem (Fermat, 1640): If $p$ is prime and $\gcd(a, p) = 1$ (i.e., $p \nmid a$), then: $$a^{p-1} \equiv 1 \pmod{p}$$

Equivalently, for any integer $a$: $a^p \equiv a \pmod{p}$.

Proof sketch 1 (counting): Consider the set $S = {a, 2a, 3a, \ldots, (p-1)a}$ reduced mod $p$. Since $\gcd(a, p) = 1$, multiplication by $a$ is a bijection on ${1, \ldots, p-1}$ — it permutes the residues. Therefore the product of all elements of $S$ mod $p$ equals $(p-1)!$ mod $p$. But the product of $S$ is $a^{p-1} \cdot (p-1)!$. Canceling $(p-1)!$ (which is coprime to $p$) gives $a^{p-1} \equiv 1$.

Proof sketch 2 (binomial theorem / induction): Prove $a^p \equiv a$ by induction using $(a+1)^p \equiv a^p + 1^p \pmod{p}$ (since the middle binomial coefficients $\binom{p}{k}$ are all divisible by $p$ for prime $p$).

Why does this matter for crypto?

• Modular inverse: $a^{-1} \equiv a^{p-2} \pmod{p}$ — no extended GCD needed when $p$ is prime!
• Fast exponentiation structure: $a^k \equiv a^{k \bmod (p-1)} \pmod{p}$, so exponents can be reduced mod $p-1$
• RSA: if $ed \equiv 1 \pmod{p-1}$, then $(m^e)^d = m^{ed} \equiv m^1 = m \pmod{p}$ — encryption and decryption are inverses

Fermat primality test: if $a^{n-1} \not\equiv 1 \pmod{n}$, then $n$ is definitely composite. Caveat: Carmichael numbers (e.g., 561) pass for all $a$ coprime to $n$ despite being composite. Miller-Rabin fixes this.

# Fermat's little theorem: a^(p-1) ≡ 1 (mod p) for prime p, gcd(a,p)=1
p = 17
for a in range(1, p):
    assert mod_exp(a, p - 1, p) == 1, f"Failed for a={a}"
print("Fermat's little theorem verified for p=17")

# Modular inverse via Fermat (only works when p is prime!)
def fermat_inverse(a: int, p: int) -> int:
    """Requires p prime. Returns a^{-1} mod p."""
    return mod_exp(a, p - 2, p)

print(fermat_inverse(3, 17))   # 6: 3*6 = 18 ≡ 1 (mod 17)

# Fermat primality test (probabilistic)
import random
def fermat_probably_prime(n: int, k: int = 20) -> bool:
    """Returns False if definitely composite, True if probably prime."""
    if n < 2: return False
    for _ in range(k):
        a = random.randrange(2, n)
        if mod_exp(a, n - 1, n) != 1:
            return False
    return True

Generalizing Fermat beyond primes

Fermat's theorem only applies when the modulus is prime. But RSA uses a composite modulus $n = pq$. Euler's theorem is the generalization we need.

Definition: Euler's totient function $\phi(n)$ counts the integers in ${1, \ldots, n}$ that are coprime to $n$: $$\phi(n) = |{k : 1 \le k \le n, \gcd(k, n) = 1}|$$

Key values:

• $\phi(p) = p - 1$ for prime $p$ (all of $1, \ldots, p-1$ are coprime)
• $\phi(p^k) = p^k - p^{k-1} = p^{k-1}(p-1)$ (exclude multiples of $p$)
• $\phi$ is multiplicative: if $\gcd(m, n) = 1$, then $\phi(mn) = \phi(m)\phi(n)$
• Therefore: $\phi(pq) = (p-1)(q-1)$ for distinct primes $p, q$

Theorem (Euler, 1763): For any $n$ and any $a$ with $\gcd(a, n) = 1$: $$a^{\phi(n)} \equiv 1 \pmod{n}$$

Proof: exactly the same permutation argument as Fermat, but over the $\phi(n)$ units of $\mathbb{Z}/n\mathbb{Z}$.

RSA key generation (the full picture):

1. Choose large primes $p, q$; compute $n = pq$, $\phi(n) = (p-1)(q-1)$
2. Choose public exponent $e$ with $\gcd(e, \phi(n)) = 1$ (commonly $e = 65537$)
3. Compute private exponent $d$ with $ed \equiv 1 \pmod{\phi(n)}$ via extended GCD
4. Encrypt: $c = m^e \bmod n$
5. Decrypt: $m = c^d \bmod n$
6. Correctness: $c^d = m^{ed} = m^{1 + k\phi(n)} = m \cdot (m^{\phi(n)})^k \equiv m \cdot 1^k = m \pmod{n}$

The security rests on: knowing $n$ but not $p, q$, you cannot compute $\phi(n) = (p-1)(q-1)$, so you cannot compute $d$. Factoring $n$ is believed to be hard (sub-exponential algorithms exist but are impractical for 2048+ bit $n$).

def euler_totient(n: int) -> int:
    """Compute phi(n) by trial factoring. O(sqrt(n))."""
    result = n
    p = 2
    temp = n
    while p * p <= temp:
        if temp % p == 0:
            while temp % p == 0:
                temp //= p
            result -= result // p   # multiply by (1 - 1/p)
        p += 1
    if temp > 1:  # temp is a prime factor
        result -= result // temp
    return result

print(euler_totient(12))   # 4: {1,5,7,11}
print(euler_totient(15))   # 8: (3-1)*(5-1)=8

# Minimal RSA (educational — NOT production-safe)
def rsa_keygen(p: int, q: int):
    n = p * q
    phi_n = (p - 1) * (q - 1)
    e = 65537
    assert gcd(e, phi_n) == 1
    d = mod_inverse(e, phi_n)
    return (e, n), (d, n)   # (public_key, private_key)

pub, priv = rsa_keygen(61, 53)
m = 42
c = mod_exp(m, pub[0], pub[1])
print(mod_exp(c, priv[0], priv[1]))  # 42 — decrypted!

Splitting a computation across independent moduli

The Chinese Remainder Theorem (CRT) states: if $n_1, n_2, \ldots, n_k$ are pairwise coprime (every pair has GCD 1), and $a_1, \ldots, a_k$ are arbitrary integers, then the system of congruences: $$x \equiv a_1 \pmod{n_1}, \quad x \equiv a_2 \pmod{n_2}, \quad \ldots, \quad x \equiv a_k \pmod{n_k}$$ has a unique solution modulo $N = n_1 n_2 \cdots n_k$.

Construction: let $N_i = N / n_i$ and $M_i = N_i^{-1} \bmod n_i$. Then: $$x \equiv \sum_{i=1}^k a_i N_i M_i \pmod{N}$$

Intuition: each term $a_i N_i M_i$ is "tuned" to be $a_i$ mod $n_i$ and 0 mod all other $n_j$ — a kind of basis. The sum hits all targets simultaneously.

Applications in CS:

1. RSA speedup: decrypt faster by computing $m_p = c^{d \bmod (p-1)} \bmod p$ and $m_q = c^{d \bmod (q-1)} \bmod q$, then combining with CRT. Exponents are half the size, giving a 4x speedup (Garner's algorithm).
2. Multi-precision arithmetic: represent large numbers as tuples of small residues; add/multiply component-wise; reconstruct with CRT. Avoids big-integer overflow.
3. Hash function design: hashing with a product of small primes reduces collision probability by independence.

The isomorphism view: CRT says $\mathbb{Z}/N\mathbb{Z} \cong \mathbb{Z}/n_1\mathbb{Z} \times \cdots \times \mathbb{Z}/n_k\mathbb{Z}$ as rings. Operations in the big ring decompose into independent operations in each small ring — a divide-and-conquer principle at the algebraic level.

def crt(remainders: list[int], moduli: list[int]) -> int:
    """Solve x ≡ r_i (mod m_i) for pairwise coprime m_i.
    Returns x in [0, product of moduli)."""
    N = 1
    for m in moduli:
        N *= m
    x = 0
    for r, m in zip(remainders, moduli):
        Ni = N // m
        Mi = mod_inverse(Ni, m)   # Ni^{-1} mod m
        x += r * Ni * Mi
    return x % N

# x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7)
# Unique solution mod 105
result = crt([2, 3, 2], [3, 5, 7])
print(result)          # 23
print(23 % 3, 23 % 5, 23 % 7)   # 2 3 2 ✓

# RSA-CRT decryption speedup
def rsa_decrypt_crt(c, d, p, q):
    dp = d % (p - 1)     # d mod (p-1)
    dq = d % (q - 1)     # d mod (q-1)
    mp = mod_exp(c, dp, p)
    mq = mod_exp(c, dq, q)
    # Garner's formula to combine
    qinv = mod_inverse(q, p)
    h = (qinv * (mp - mq)) % p
    return mq + h * q

When can you divide in modular arithmetic?

Division in modular arithmetic means multiplication by the inverse. The equation $ax \equiv b \pmod{n}$ has a solution if and only if $\gcd(a, n) \mid b$. When this condition holds, there are exactly $\gcd(a,n)$ distinct solutions mod $n$.

Special case: if $\gcd(a, n) = 1$, there is exactly one solution: $x \equiv a^{-1} b \pmod{n}$.

Three ways to compute $a^{-1} \bmod n$:

1. Extended Euclidean algorithm: general, works for any $n$, runs in $O(\log n)$
2. Fermat's little theorem ($a^{p-2} \bmod p$): only when $n = p$ is prime, but elegant and easy to remember
3. Euler's theorem ($a^{\phi(n)-1} \bmod n$): works for any $n$ with known $\phi(n)$, but computing $\phi(n)$ requires factoring $n$ — practical only if you know the factorization

Why inverses matter:

• Cryptography: RSA private key $d = e^{-1} \bmod \phi(n)$; ECDSA signing uses $k^{-1}$
• Linear algebra mod p: Gaussian elimination over $\mathbb{F}_p$ requires division, i.e., modular inverses
• Competitive programming: dividing factorials in combinatorics problems mod a prime (computing $\binom{n}{k} \bmod p$ for large $n$)
• Hash functions: polynomial rolling hashes sometimes need inverse for substring queries

Batch inverse computation: computing $n$ inverses naively takes $O(n \log M)$ time ($M$ = modulus). A clever trick using prefix products computes all $n$ inverses with just one modular inverse call plus $O(n)$ multiplications.

# Batch modular inverses: compute inv[i] = i^{-1} mod p for i in 1..n
# Uses: prefix products + one gcd call + suffix scan
def batch_inverses(n: int, p: int) -> list[int]:
    """Returns list where result[i] = i^{-1} mod p for i = 1..n."""
    inv = [0] * (n + 1)
    inv[1] = 1
    for i in range(2, n + 1):
        # inv[i] = -(p // i) * inv[p % i] mod p
        # Derived from: p = (p//i)*i + (p%i), reduce mod p
        inv[i] = -(p // i) * inv[p % i] % p
    return inv

p = 1_000_000_007
invs = batch_inverses(10, p)
for i in range(1, 11):
    assert i * invs[i] % p == 1, f"inv[{i}] is wrong"
print("All batch inverses verified")
print(f"inv[7] mod {p} = {invs[7]}")
# Check: 7 * inv[7] % p == 1

Number theory in every hash table and random number generator

Polynomial (rolling) hashing: treat a string $s_0 s_1 \ldots s_{n-1}$ as a polynomial evaluated at a base $b$: $$H(s) = \left(\sum_{i=0}^{n-1} s_i \cdot b^{n-1-i}\right) \bmod p$$ for a prime $p$. Why prime? Because for a prime modulus, the polynomial has at most $\deg$ roots mod $p$ — the number of "collisions" for distinct strings of length $n$ is at most $n/p$. Choosing $p \approx 10^{18}$ makes collisions astronomically unlikely.

Why use a prime table size in hash tables? If your hash table has $m$ entries and you insert keys via $k \bmod m$, prime $m$ avoids the systematic clustering that occurs with powers-of-2 table sizes (since most real key distributions have structure at powers of 2). This is a practical consequence of divisibility theory.

Linear Congruential Generator (LCG): the workhorse PRNG (used in C's rand()): $$X_{n+1} = (aX_n + c) \bmod m$$ The Hull-Dobell theorem characterizes when the LCG has full period $m$: (1) $\gcd(c, m) = 1$; (2) $a \equiv 1 \pmod{p}$ for every prime $p \mid m$; (3) if $4 \mid m$ then $a \equiv 1 \pmod{4}$. All three conditions are number-theoretic.

Diffie-Hellman key exchange is another direct application: choose a large prime $p$ and a primitive root $g$ (a generator of the multiplicative group mod $p$). Alice sends $g^a \bmod p$, Bob sends $g^b \bmod p$, both compute the shared secret $g^{ab} \bmod p$. Security rests on the discrete logarithm problem: recovering $a$ from $g^a \bmod p$ is believed to be hard.

# Polynomial rolling hash with a large prime
MOD = (1 << 61) - 1   # Mersenne prime: fast reduction
BASE = 131

def poly_hash(s: str) -> int:
    h = 0
    for c in s:
        h = (h * BASE + ord(c)) % MOD
    return h

# Rabin-Karp substring search using rolling hash
def rabin_karp(text: str, pattern: str) -> list[int]:
    n, m = len(text), len(pattern)
    if m > n: return []
    ph = poly_hash(pattern)
    # Precompute BASE^(m-1) mod MOD for rolling
    power = pow(BASE, m - 1, MOD)
    th = poly_hash(text[:m])
    matches = []
    if th == ph and text[:m] == pattern:
        matches.append(0)
    for i in range(1, n - m + 1):
        # Roll: remove leftmost char, add new rightmost
        th = (th - ord(text[i-1]) * power) % MOD
        th = (th * BASE + ord(text[i+m-1])) % MOD
        if th == ph and text[i:i+m] == pattern:
            matches.append(i)
    return matches

print(rabin_karp("abcabcabc", "abc"))  # [0, 3, 6]

The practical primality test for cryptographic numbers

Trial division is $O(\sqrt{n})$ — completely infeasible for 1024-bit numbers. Miller-Rabin is a probabilistic test running in $O(k \log^2 n)$ time (for $k$ witnesses) that is the industry standard for primality checking.

Setup: Write $n - 1 = 2^s \cdot d$ with $d$ odd (factor out all 2s). Fermat says $a^{n-1} \equiv 1 \pmod n$ for prime $n$ and $\gcd(a,n)=1$. But more: $a^{n-1} - 1 = (a^{(n-1)/2} - 1)(a^{(n-1)/2} + 1)$. Since $\mathbb{Z}/p\mathbb{Z}$ has no nontrivial square roots of 1, either $a^{(n-1)/2} \equiv 1$ or $a^{(n-1)/2} \equiv -1 \pmod p$. This propagates down the squaring sequence.

The test: given witness $a$, compute $x = a^d \bmod n$. Then repeatedly square: $x, x^2, x^4, \ldots, x^{2^s} = a^{n-1}$. If $n$ is prime, either $x = 1$ initially, or at some point $x = n-1$ (i.e., $\equiv -1$), after which all subsequent values are 1. If neither occurs before hitting a non-1, non-(-1) value and then squaring to 1, $n$ is composite for certain.

Error probability: a composite $n$ fools at most $1/4$ of all witnesses. With $k$ independent witnesses, the false-positive probability is $\le (1/4)^k$. For $k = 40$, this is $\approx 10^{-24}$ — completely negligible.

Deterministic variant: for $n < 3{,}317{,}044{,}064{,}679{,}887{,}385{,}961{,}981$, testing witnesses ${2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}$ is sufficient (proved by computation). Python's sympy.isprime and most crypto libraries use this or the BPSW test.

import random

def miller_rabin(n: int, k: int = 40) -> bool:
    """Probabilistic primality test. Returns False if definitely
    composite, True if probably prime (error prob <= 4^{-k})."""
    if n < 2: return False
    if n in (2, 3): return True
    if n % 2 == 0: return False

    # Write n-1 = 2^s * d
    s, d = 0, n - 1
    while d % 2 == 0:
        s += 1
        d //= 2

    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = mod_exp(a, d, n)
        if x in (1, n - 1):
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False   # Composite witness found!
    return True            # Probably prime

# Generate a random probable prime of bit_length bits
def gen_prime(bits: int) -> int:
    while True:
        p = random.getrandbits(bits) | (1 << bits - 1) | 1  # odd, top bit set
        if miller_rabin(p, 40):
            return p

print(miller_rabin(104729))   # True (104729 is prime)
print(miller_rabin(104728))   # False