Number-Theoretic Algorithms & RSA

By the end of this lesson you will be able to trace euclid, extendedEuclid (including its back-substitution of x and y), modularLinearEquationSolver, the Chinese remainder theorem, modularExponentiation (bit by bit), RSA key generation / encryption / decryption, pseudoprime, witness, millerRabin and pollardRho line-by-line with every number visible; read a modular clock, a CRT table and the group Z_n^* with the order of each element; implement each procedure in idiomatic Dart (and know exactly where 64-bit int overflows and BigInt is needed); and explain why fast gcd, fast powers and fast primality testing are what make public-key cryptography practical while factoring stays hard. Nothing is assumed beyond school arithmetic and the earlier lessons: What an algorithm is (and its cost), Insertion and merge sort (loop invariants) and Growth of functions (O-notation).

Divisibility, primes and the modular clock (Z_n)

A number-theoretic algorithm is a procedure whose inputs are whole numbers (integers) and whose cleverness comes from properties of divisibility and remainders. Before any procedure, four words: divides, prime, gcd and modulus.

Divisibility. We write d | a ("d divides a") when a = k·d for some integer k; for example 3 | 12 because 12 = 4·3. The greatest common divisor gcd(a, b) is the largest integer dividing both (gcd(12, 8) = 4). A number p > 1 is prime if its only positive divisors are 1 and p; otherwise it is composite (12 = 3·4). Two numbers are relatively prime (coprime) when gcd(a, b) = 1.

The remainder operator. a mod n is the remainder r in 0 ≤ r < n with a = ⌊a/n⌋·n + r (⌊x⌋ means "round down"). Two numbers are congruent modulo n, written a ≡ b (mod n), when they have the same remainder, that is when n | (a − b). The number n is called the modulus and the set of clock positions {0, 1, …, n−1} is written Z_n.

Small example. On the 12-clock: 10 + 5 = 15 ≡ 3 and 7 · 8 = 56 ≡ 8 (56 = 4·12 + 8). Also 15 ≡ 27 ≡ 3 (mod 12) because 15, 27 and 3 all leave remainder 3. Try the player below: it walks around a clock in jumps of a, which is exactly the picture behind modularLinearEquationSolver later.

The same walk in Dart (the loop ends when the hand returns to the start, which happens after n / gcd(step, n) jumps):

What to notice. Jumping a hours at a time around an n-hour clock visits exactly n / gcd(a, n) different hours, and then repeats. Jump 5 on the 12-clock: gcd(5, 12) = 1, so all 12 hours are visited. Jump 4: gcd(4, 12) = 4, so only 3 hours (0, 4, 8). The animation's last frame states that count; it is also the "order" of a in the additive group Z_n.

In Dart. % always returns a value in 0 ≤ r < n for n > 0 (so -7 % 5 == 3), whereas .remainder and ~/ truncate toward zero ((-7).remainder(5) == -2, -7 ~/ 5 == -1) — the maths "mod" of this lesson is Dart's %.

The group Z_n^* and the order of an element

Formally Z_n^* = {a in 0..n−1 : gcd(a, n) = 1}. It is closed under multiplication mod n, has identity 1, and every element has a multiplicative inverse (a number x with a·x ≡ 1 mod n), so it is a group of size φ(n) (Euler's totient); Z₁₂^* = {1, 5, 7, 11}, φ(12) = 4, and for a prime p, Z_p^* = {1, …, p−1} has φ(p) = p − 1 members. The order of a, ord(a), is the smallest t ≥ 1 with a^t ≡ 1 (mod n): the number of multiplications until you are back at 1. Lagrange: ord(a) always divides φ(n) — that is why Euler's theorem holds: a^φ(n) ≡ 1 (mod n) (prime special case, Fermat: a^p−1 ≡ 1 (mod p)). If one element has order exactly φ(n) it is a generator (primitive root) and the group is cyclic.

Small example. Z₇^*: powers of 3 are 3, 2, 6, 4, 5, 1 — six steps to reach 1, so ord(3) = 6 = φ(7), and 3 is a generator. Powers of 2 are 2, 4, 1: ord(2) = 3, which divides 6 ✓. In Z₁₂^* = {1, 5, 7, 11} every non-identity element has order 2 and none has order 4, so Z₁₂^* has no generator.

What to notice. Each frame lists the powers of one element until they return to 1, then colours the element by whether its order equals φ(n) (green = generator). Every order in the last frame divides φ(n).

Dart (helper used by the group player; the loop stops when the power returns to 1, which is guaranteed for gcd(a, n) = 1):

Z_n^* is a group of order φ(n), and every element order divides it (group laws and Lagrange's theorem), checked for every n ≤ 60:

The cyclic-group theorem (when is Z_n^* cyclic?) as two functions that must agree for every n ≤ 300:

euclid(a, b)

The whole algorithm rests on the gcd recursion theorem: gcd(a, b) = gcd(b, a mod b). Every common divisor of a and b also divides the leftover a mod b, and vice versa, so replacing the pair by a smaller pair never changes the gcd — and when b reaches 0, gcd(a, 0) = a. Dart index note: nothing here is an array, so there is no index shift; the three lines of the panel map to one if and one recursive return.

Small example. gcd(12, 8): 12 mod 8 = 4, so gcd(12, 8) = gcd(8, 4); 8 mod 4 = 0, so = gcd(4, 0) = 4.

What to notice. The second argument shrinks every call (it becomes a mod b, which is always smaller than b) so the recursion must stop, and when b = 0 the first argument is the answer. When a < b (Example 2) the first call just swaps the numbers, since a mod b = a.

Correctness idea (loop-invariant style): the invariant "gcd of the current pair = gcd of the original pair" holds initially and is preserved by the gcd recursion theorem; the recursion terminates because b strictly decreases (a mod b < b).

Dart implementation:

Input size → what is feasible. a, b ≤ 1018 → at most ≈ 87 divisions (Lamé), instant; 106 such gcds ≈ 9·10⁷ steps ≈ 1 s. Only the naive divisor scan (10¹⁸ steps) is hopeless.

Lamé’s theorem and the O(lg b) bound in code (Fibonacci list indexed by the subscript):

The GCD recursion theorem, Bézout and the common-divisor corollary checked against the definition of gcd:

Euler’s theorem, Fermat’s theorem and φ(n) = n·Π(1 − 1/p) as executable checks:

extendedEuclid(a, b)

Bézout's identity says gcd(a, b) is always an integer combination d = a·x + b·y. The procedure recurses exactly like euclid, then on the way back up converts the recipe for (b, a mod b) into a recipe for (a, b): since a mod b = a − ⌊a/b⌋·b, we get d = b·x′ + (a − ⌊a/b⌋b)·y′ = a·y′ + b·(x′ − ⌊a/b⌋·y′). Those are lines 4–5, and it is the back-substitution you may know from school: rewrite each remainder in terms of the two numbers before it. This is also how modular inverses are found: if gcd(a, n) = 1 then x (reduced mod n) is a⁻¹ mod n.

Small example. gcd(12, 8) = 4. Calls: (12, 8) → (8, 4) → (4, 0). Bottom: (4, 1, 0). Level (8, 4): ⌊8/4⌋ = 2, so (d, x, y) = (4, 0, 1 − 2·0) = (4, 0, 1), i.e. 4 = 8·0 + 4·1. Level (12, 8): ⌊12/8⌋ = 1 so (4, 1, 0 − 1·1) = (4, 1, −1): 4 = 12·1 + 8·(−1) ✓.

First the recursion itself, with the back-substitution column filled in as each call returns (the last column checks d = a·x + b·y at every level):

Now the same computation read the way it is done on paper: write the division equations rk−1 = qk·rk + rk+1 going down, then start from the last non-zero remainder (the gcd) and substitute upward, keeping the expression in terms of the two most recent remainders until only a and b are left. The coefficient update in each substitution is exactly lines 4–5 of the procedure.

What to notice. The two views give identical (x, y). In the paper view, the pair of coefficients (s, t) moves from (1, −q) to (t, s − q′·t) at each substitution — the same "swap and subtract q times" as lines 4–5.

Dart implementation (returns a record; ~/ is integer division, which equals ⌊a/b⌋ for non-negative a, b) and an iterative form that never recurses (it keeps the two most recent coefficient pairs — the same back-substitution done on the way down):

Input size → what is feasible. a, b ≤ 1018 → ≤ 87 recursion levels, coefficients stay below b, so no overflow; recursion depth 87 is safe.

modularLinearEquationSolver(a, b, n)

The equation a·x ≡ b (mod n) has solutions iff d = gcd(a, n) divides b, and then it has exactly d solutions modulo n: x₀, x₀ + n/d, x₀ + 2n/d, … The procedure gets d and the Bézout number x′ (with d = a·x′ + n·y′) from extendedEuclid, then x₀ = x′·(b/d) mod n is the first solution.

Small example. 3x ≡ 6 (mod 9): d = gcd(3, 9) = 3 divides 6, so there are 3 solutions. extendedEuclid(3, 9) = (3, 1, 0), so x₀ = 1·(6/3) mod 9 = 2, and the solutions are 2, 2 + 9/3 = 5, 8 (check: 3·8 = 24 = 2·9 + 6 ✓). When d = 1 there is a single solution, x = a⁻¹·b mod n.

What to notice. In the "a·x mod n" column every solution gives b; solutions are spaced n/d apart; in Example 2 the test "d | b" fails so nothing is printed.

Dart implementation (returns the list instead of printing; empty list means "no solutions"; the multiplication x1 * (b ~/ d) can overflow for n near 2³¹ or more, in which case use BigInt):

Input size → what is feasible. n ≤ 2·109 → one extendedEuclid plus d ≤ n solutions to list; if d is huge (e.g. n = 10⁹, a = 0) the output itself has d numbers.

Count of solutions: exactly d = gcd(a, n) solutions or none, compared with brute force for every n ≤ 40:

The Chinese remainder theorem (the CRT table)

The Chinese remainder theorem. If n₁, …, n_k are pairwise relatively prime and n = n₁n₂⋯n_k, then the map a ↦ (a mod n₁, …, a mod n_k) is a one-to-one correspondence (a bijection: every input has its own output and every output is hit) between Z_n and Z_n₁ × ⋯ × Z_nk. It can be proved constructively with m_i = n/n_i and c_i = m_i·(m_i⁻¹ mod n_i): c_i is ≡ 1 mod n_i and ≡ 0 mod every other n_j, so a = (Σ a_ic_i) mod n has the right remainder on every clock. The panel below writes the construction as steps, and the inverse in line 5 comes from extendedEuclid.

Small example. x ≡ 1 (mod 2) and x ≡ 2 (mod 3): n = 6, m₁ = 3, 3⁻¹ mod 2 = 1, c₁ = 3; m₂ = 2, 2⁻¹ mod 3 = 2, c₂ = 4; x = (1·3 + 2·4) mod 6 = 11 mod 6 = 5. Check: 5 mod 2 = 1, 5 mod 3 = 2 ✓.

The CRT table

Picture Z_n ↔ Z_n₁ × Z_n₂ as a grid: row = a mod n₁, column = a mod n₂, and the cell holds the unique a in 0..n−1. Every cell is filled exactly once — that is the bijection. The next players draw the grid (for n ≤ 130), pick an a, show the cell (a mod n₁, a mod n₂) it lands in, then go backwards from the cell to a with the c_i formula.

And the Z₁₅ ↔ Z₃ × Z₅ table in full (computed by the page, not typed): each a appears once, and adding or multiplying a's mod 15 corresponds to adding or multiplying the residue pairs componentwise.

The construction itself, step by step, in the animation below:

What to notice. Each c_i is 1 on its own clock and 0 on the others; the answer is a sum of "one clock at a time" contributions. Example 2 (moduli 4 and 6) shows the theorem refusing to apply: they share the factor 2.

Dart implementation (the product big *= ni is where 64-bit overflow would occur if the moduli multiply past 2⁶³; the page's own inputs are far smaller):

Input size → what is feasible. moduli ≤ 10⁹ with n = Πn_i ≤ 9·10¹⁸ → fits an int; beyond that use BigInt. Each modulus costs one extendedEuclid.

The CRT isomorphism Z₁₀₅ ≅ Z₃×Z₅×Z₇ (bijection, respects + and ·) and its failure for the non-coprime moduli 4 and 6:

modularExponentiation(a, b, n)

Reading the bits of b from the most significant, each iteration squares the running result d (doubling the exponent seen so far, tracked by c) and, if the bit is 1, multiplies once by a (adding 1 to the exponent). After all bits, d = a^b mod n. The variable c is not needed to get the answer — it is kept only to state the loop invariant "d = a^c mod n"; the Dart panel drops it. Binary numbers: b = 13 is 1101 because 13 = 8 + 4 + 0 + 1; the leftmost bit is the "most significant".

Small example. 3¹³ mod 7, bits 1101: start d = 1. Bit 1: square (1), multiply by 3 → 3. Bit 1: square 9 mod 7 = 2, multiply by 3 → 6. Bit 0: square 36 mod 7 = 1, no multiply → 1. Bit 1: square 1, multiply by 3 → 3. Answer 3 (four squarings and three multiplies instead of twelve multiplications). The pictures below show the bits of b as boxes; the box being read is highlighted and the S / M strip records "square" and "multiply".

What to notice. Every bit causes exactly one square (S); only the 1-bits also cause a multiply (M) — so the S/M strip is literally the binary expansion of b read left to right. c follows the bits already read, and d always equals a^c mod n.

Dart implementation (bitLength gives k + 1, the number of bits; for b = 0 the loop simply does not run and d stays 1), and two overflow-safe variants — one for moduli up to 2⁶² using the doubling mulMod from the question bank, and one with BigInt, which has no size limit (BigInt.modPow is the library's own modularExponentiation):

Input size → what is feasible. b ≤ 1018 → ≤ 120 multiplications; but every product must fit: n ≤ 3·109 with plain int, larger n needs mulMod or BigInt (question 20).

The multiplication count β + popcount(b) ≤ 2β, in code:

The RSA public-key cryptosystem

Key generation: pick two distinct primes p and q; n = pq; φ(n) = (p−1)(q−1); pick a small odd e with gcd(e, φ(n)) = 1; compute d = e⁻¹ mod φ(n) with extendedEuclid. Public key P = (e, n), secret key S = (d, n). Encrypt message M < n: C = P(M) = M^e mod n. Decrypt: S(C) = C^d mod n = M^ed mod n = M, because ed ≡ 1 (mod φ(n)) and Euler's theorem (the full correctness argument also covers a message M that shares a factor with n). All the heavy lifting is modularExponentiation. The first player below uses p = 17, q = 23, e = 3, d = 235.

Small example. p = 3, q = 11: n = 33, φ = 2·10 = 20, e = 3 (gcd(3, 20) = 1), d = 7 because 3·7 = 21 ≡ 1 (mod 20). Message M = 4: C = 4³ mod 33 = 64 mod 33 = 31. Decrypt: 31⁷ mod 33 = 4 ✓. The players below repeat this with p = 17 and q = 23, an edge case, your own numbers, and finally a whole word one letter at a time.

Round trip on text: each letter A..Z becomes 1..26, is encrypted separately and decrypted back (the modulus must exceed 26). This is only to show the mechanism; separate deterministic letters would be trivially breakable by frequency counting in real life.

What to notice. Only two heavy operations exist: extendedEuclid once (to get d) and modularExponentiation twice per message. Ciphertext letters are unrelated-looking numbers; only d turns them back.

Dart implementation:

Input size → what is feasible. toy primes below 5·10⁴ fit int; real keys (2048-bit n) need BigInt (question 37).

RSA correctness for every message M in Z_n, and why n must be squarefree:

Primality testing: pseudoprime, witness and millerRabin

Fermat's test says: if n is prime then aⁿ⁻¹ ≡ 1 (mod n) for every a not divisible by n. Testing this with a = 2 gives the first procedure, pseudoprime. Its errors are called base-2 pseudoprimes: composite n with 2ⁿ⁻¹ ≡ 1 (mod n). The converse fails badly for Carmichael numbers such as 561 = 3·11·17 that satisfy a⁵⁶⁰ ≡ 1 for every a coprime to 561 yet are composite (for example 2⁵⁶⁰ mod 561 = 1). witness repairs this with a second check: writing n − 1 = 2^t·u with u odd, it computes a^u, then squares t times. If a squaring produces 1 from a number other than 1 or n − 1, that number is a nontrivial square root of 1, which cannot exist modulo a prime — so n is definitely composite. If the final value is not 1, Fermat's test failed — also definitely composite.

pseudoprime(n)

The procedure is three lines: if 2ⁿ⁻¹ mod n ≠ 1 return COMPOSITE (definitely); else return PRIME ("we hope!"). Only one kind of error is possible: answering PRIME for a composite. Small example. n = 15: 2¹⁴ mod 15 = 4 ≠ 1 → COMPOSITE, certain. n = 341 = 11·31: 2³⁴⁰ mod 341 = 1 → answers PRIME although 341 is composite (the smallest base-2 pseudoprime). The player shows the modular exponentiation of 2ⁿ⁻¹ bit by bit.

What to notice. "COMPOSITE" always comes with a proof (the value 2ⁿ⁻¹ mod n ≠ 1); "PRIME" is only a guess. The table below (computed by the page) counts the guesses that are wrong.

Dart implementation:

Input size → what is feasible. n < 3·109 with int arithmetic; one modular exponentiation ≈ 60 multiplications.

Prime number theorem π(n) ~ n/ln n and the density of primes among random odd numbers:

How often does pseudoprime lie? Counting base-2 pseudoprimes and Carmichael numbers (Korselt) below 10⁶:

witness(a, n)

Small example. n = 13 gives n − 1 = 12 = 2²·3, so t = 2, u = 3. For a = 5: x₀ = 5³ mod 13 = 8, x₁ = 8² mod 13 = 12 (= n − 1), x₂ = 12² mod 13 = 1 — Fermat holds and the only "1 after a square" came from n − 1 = −1, a legal square root, so a is not a witness. Here is the sequence a^u, a^2u, …, a^2tu drawn as boxes:

What to notice. The boxes x₀ … x_t are successive squarings. A witness answer true comes from either (1) a box equal to 1 right after a box that is neither 1 nor n − 1 (red pair: nontrivial square root of 1) or (2) a last box not equal to 1 (Fermat fails).

Dart implementation:

Input size → what is feasible. n < 3·109 in plain int; use BigInt for 64-bit n. Cost ≈ log n multiplications.

Square roots of 1, the fact witness exploits:

millerRabin(n, s)

millerRabin tries s random bases. Witness-count theorem: if n is an odd composite, at least (n − 1)/2 of the possible values of a are witnesses, so each round misses with probability ≤ 1/2 and the answer "PRIME" is wrong with probability ≤ 2^−s. "COMPOSITE" is always right. The page uses a small seeded pseudo-random generator, so a given (n, s, seed) always gives the same trace. Small example. n = 15 with a = 2: 14 = 2·7, x₀ = 2⁷ mod 15 = 8, x₁ = 64 mod 15 = 4 ≠ 1, so true at once; millerRabin(15, s) returns COMPOSITE in the first round.

What to notice. A composite is usually exposed in round 1 (see the witness counts below); primes are never exposed, so all s rounds run and the answer is PRIME.

How many witnesses does a composite really have? Counted by brute force by the page (bases 1..n−1):

Dart implementation (an exact all-witness check is used in the verify file; here rng supplies the randomness, and the x * x inside witness stays exact only while n < 3·10⁹; the deterministic 64-bit variant in the question bank uses BigInt):

Input size → what is feasible. s = 30 rounds on a 64-bit n ≈ 30·64 ≈ 2000 multiplications; error ≤ 2⁻³⁰ ≈ 10⁻⁹.

Error bound: error ≤ 2^−s, from the exact witness counts of every odd composite below 2000, plus a seeded simulation:

pollardRho factoring (the ρ shape)

The sequence is x_i = (x_i−1² − 1) mod n. Rather than compare every pair, the procedure remembers only y = x_k for k = 1, 2, 4, 8, … (lines 2, 12–13) and takes gcd(y − x_i, n) each step (line 8). The loop is while true — it may print several factors and never halts on its own; the animations stop after the first printed factor, or after a step cap.

Small example. n = 91 = 7·13, x₁ = 2: the sequence is 2, 3, 8, 63, 55, 21, 76, 42, 34, 63, 55, 21, … — it comes back to 63 at x₁₀ = x₄: a tail 2, 3, 8 and a cycle 63, 55, 21, 76, 42, 34. Modulo 7 the same numbers are 2, 3, 1, 0, 6, 0, 6, … which loops much earlier (tail 2, 3, 1, then the cycle 0, 6), and at that moment the difference is a multiple of 7 while it is not a multiple of 13, so gcd = 7.

First the ρ shape itself: the sequence reduced mod the smallest prime factor p of n, laid out with the tail on a line and the cycle on a circle (nodes appear in order of i). The frame where a value repeats is the moment when pollardRho's gcd can reveal p:

Then the full procedure, with the checkpoint y doubling at i = 2, 4, 8, …:

What to notice. In the ρ picture the loop mod p closes after about √p steps; the checkpoint y at k = 2, 4, 8, … eventually sits inside the cycle so some x_i ≡ y (mod p) is found while x_i ≢ y (mod n/p) — then d = gcd(y − x_i, n) is a proper divisor. For a prime n there is nothing to find.

Dart implementation (maxSteps bounds the otherwise infinite loop and returns null on failure):

Input size → what is feasible. n < 262 → ≈ n^1/4 ≈ 4.6·10⁴ iterations (heuristic); a prime n never prints anything, so test primality first.

The √p birthday bound for x ↦ x² − 1 mod p, and why other constants are worse:

Quiz

Interview questions

Cheat sheet