Lesson 10 · Core algorithms

Number Theory Basics for Algorithms

Number-theory tools reduce large arithmetic problems to invariants: Euclid preserves common divisors, modular arithmetic preserves congruence, and exponentiation by squaring follows the exponent's binary digits.

Euclid's Algorithm

The greatest common divisor of two integers is the largest number dividing both. Finding it by testing every candidate from the smaller number downward is O(min(a, b)) — impractical for large values and unnecessary.

Euclid's algorithm rests on one identity: gcd(a, b) = gcd(b, a mod b), with gcd(a, 0) = a as the base case.

The identity holds because any common divisor of a and b also divides a mod b. Writing a = q·b + r, a divisor d of both a and b divides a − q·b, which is r. The argument runs in reverse too, so the pair (a, b) and the pair (b, r) have exactly the same set of common divisors — and therefore the same greatest one. The remainder throws away no information relevant to the answer.

The algorithm is three lines: while b is non-zero, replace (a, b) with (b, a mod b); return a. It can be written recursively or iteratively with equal ease.

The complexity is O(log(min(a, b))). The bound comes from noticing that a mod b is less than half of a whenever b ≤ a/2, and when b > a/2 the remainder is a − b, which is also less than a/2 — so the larger value at least halves every two steps. The worst case is consecutive Fibonacci numbers, which is a useful fact to know since it is the standard exam question about the bound.

Two boundary cases need deciding. gcd(0, 0) is conventionally 0, though it is arguably undefined. And negative inputs should normally be taken as absolute values first, since divisibility is unaffected by sign but languages differ on the sign of % — C, C++ and Java follow the dividend, Python follows the divisor.

gcd(a, b) = gcd(b, a mod b), with gcd(a, 0) = a
The remainder preserves the entire set of common divisors
O(log min(a, b)); the worst case is consecutive Fibonacci numbers
Take absolute values first — % sign conventions differ by language

LCM and the Overflow Trap

The least common multiple is the smallest positive integer divisible by both operands, and it relates to the gcd by the identity gcd(a, b) × lcm(a, b) = a × b.

That gives lcm(a, b) = a × b / gcd(a, b) — and writing it that way is a bug.

The problem is that a × b is computed first, and for two values near the 64-bit limit that product overflows long before the division shrinks it back. The multiplication happens at full magnitude even though the final answer may be small.

The fix is to divide before multiplying: lcm = a / gcd(a, b) * b. This is exact rather than approximate, because the gcd divides a evenly by definition, so no precision is lost. The intermediate value is now at most the size of the answer itself.

The order of operations matters and is easy to get wrong when transcribing — a / (gcd * b) is a different and incorrect expression. Parenthesise deliberately.

For more than two numbers, both operations fold associatively: gcd(a, b, c) = gcd(gcd(a, b), c), and the same for lcm. Note that the gcd of a list can only shrink as elements are added, so an early exit is available once it reaches 1 — no further element can reduce it.

The lcm, by contrast, grows quickly, and the lcm of the numbers 1 through 20 already exceeds 200 million. Problems computing an lcm over many values usually require it modulo something, which introduces its own difficulty since the identity above uses division and division does not survive a modulus.

Two edge cases: lcm(a, 0) is 0, and both operands being zero makes the gcd 0, so the division must be guarded.

gcd × lcm = a × b, but the naive formula overflows
Write a / gcd(a,b) * b — exact, since the gcd divides a evenly
gcd and lcm both fold associatively over a list
gcd shrinks toward 1; lcm grows very fast

Modular Arithmetic

Modular arithmetic keeps values inside a bounded range by taking remainders, and it is what makes large-number computations tractable. Problems asking for an answer 'modulo 10⁹ + 7' are asking you to work entirely within it.

Three operations distribute over the modulus, and these are the rules to know: (a + b) mod m = ((a mod m) + (b mod m)) mod m, and the same for subtraction and multiplication. So a long computation can take the modulus at every step, keeping every intermediate small and avoiding overflow entirely.

Division does not distribute. (a / b) mod m is not ((a mod m) / (b mod m)) mod m, and this is the single most important thing to know here. Division must be replaced by multiplication by the modular inverse: a / b mod m becomes a × b⁻¹ mod m, where b⁻¹ is the value satisfying b × b⁻¹ ≡ 1 (mod m).

The inverse exists only when b and m are coprime. When m is prime — as 10⁹ + 7 is, which is why it is chosen — Fermat's little theorem gives it directly: b⁻¹ ≡ b^(m−2) (mod m), computable by fast exponentiation. For a composite modulus, the extended Euclidean algorithm finds the inverse when one exists.

Two implementation traps recur. Subtraction can go negative in languages where % follows the sign of the dividend, so normalise with ((x % m) + m) % m. Omitting this produces sporadic wrong answers that are hard to trace, and it is the same fix needed in rolling hashes.

And multiplication can overflow before the modulus is applied. With m near 10⁹, the product of two reduced values approaches 10¹⁸ — which fits a signed 64-bit integer but not a 32-bit one. Using int where long is required silently corrupts everything downstream.

Addition, subtraction and multiplication distribute over the modulus
Division does not — multiply by the modular inverse instead
With a prime modulus, b⁻¹ = b^(m−2) by Fermat's little theorem
Normalise negatives with ((x % m) + m) % m, and use 64-bit types

Implementation

Euclid, LCM, and binary modular power

import math

a, b = 84, 30
gcd = math.gcd(a, b)
lcm = a // gcd * b

print(f"gcd={gcd}, lcm={lcm}, 3^13 mod 7={pow(3, 13, 7)}")

#include <iostream>
#include <vector>
#include <queue>
#include <map>
#include <set>
#include <string>
#include <algorithm>
#include <functional>
#include <tuple>
#include <array>
#include <numeric>
using namespace std;
long long pw(long long a,int n,int m)
{
    long long r=1;
    for(;n;n/=2,a=a*a%m)if(n%2)r=r*a%m;
    return r;
}
int main()
{
    long long a=84,b=30,g=gcd(a,b);
    cout<<"gcd="<<g<<", lcm="<<a/g*b<<", 3^13 mod 7="<<pw(3,13,7);
}

class Main
{
    static long gcd(long a,long b)
    {
        while(b!=0)
        {
            long r=a%b;
            a=b;
            b=r;
        }
        return Math.abs(a);
    }
    static long pw(long a,int n,long m)
    {
        long r=1;
        for(;n>0;n/=2,a=a*a%m)if(n%2==1)r=r*a%m;
        return r;
    }
    public static void main(String[]z)
    {
        long a=84,b=30,g=gcd(a,b);
        System.out.print("gcd="+g+", lcm="+(a/g*b)+", 3^13 mod 7="+pw(3,13,7));
    }
}

Watch it run

Step through it

Running on 84, 30; 3^13 mod 7

Output

Read all 15 Steps

Start Euclid Begin with gcd(84,30).
First remainder 84 mod 30 is 24.
Replace the pair The invariant gives gcd(30,24).
Second remainder 30 mod 24 is 6.
Replace again Continue with gcd(24,6).
Third remainder 24 mod 6 is zero.
GCD reached The last nonzero divisor is 6.
Compute LCM safely Divide first: 84/6×30=420.
Read exponent bits 13 is binary 1101.
Consume bit 1 The low bit is one, so result becomes 3 mod 7.
Square base The base squares to 2 mod 7.
Consume bit 2 This bit is zero; result stays 3.
Square again The base squares to 4, and the bit contributes: result becomes 5.
Consume bit 4 Square to 2; the next bit contributes and result returns to 3.
Consume bit 8 The exponent is exhausted; report all three values.

Exponentiation by Squaring

Computing a^n by multiplying a by itself n times is O(n), which is hopeless when n is 10¹⁸ — a value that appears routinely in problems using modular arithmetic.

Exponentiation by squaring — also called binary or fast exponentiation — reduces this to O(log n) using one observation: a^n = (a^(n/2))² when n is even, and a × a^(n−1) when n is odd.

Each step either halves the exponent or decrements it, and a decrement is always followed by a halving, so the exponent reaches zero in at most 2·log₂ n steps. Computing a^(10¹⁸) takes about 120 multiplications instead of a quintillion.

The iterative form is the one to write, since it avoids recursion depth entirely. Keep a result at 1 and a running base at a. While the exponent is non-zero: if its lowest bit is set, multiply result by base; then square base and shift the exponent right by one. Taking the modulus after every multiplication keeps all values bounded.

Reading it in terms of bits makes it obvious rather than clever: the loop walks the binary representation of n, and base holds a^(2^k) at step k. Multiplying in exactly the powers corresponding to set bits reconstructs a^n, because n is the sum of those powers of two.

The same doubling structure applies beyond integers. Matrix exponentiation raises a transition matrix to the nth power in O(log n) matrix multiplications, which computes the nth Fibonacci number in O(log n) and solves linear recurrences generally. The requirement is only that the operation be associative, which is why the technique also underlies fast computation on monoids.

Where these tools go together: modular exponentiation is the core operation of RSA, which relies on it being fast to compute a^e mod n and infeasible to reverse without the factorisation. Fermat's little theorem built on it gives both the modular inverse above and the Miller–Rabin primality test. And the extended Euclidean algorithm, which returns coefficients x and y satisfying ax + by = gcd(a, b), is what solves linear Diophantine equations and the Chinese remainder theorem.

The costs worth remembering
Operation	Naive	Better
gcd(a, b)	O(min(a,b))	O(log min(a,b)) — Euclid
a^n	O(n)	O(log n) — squaring
a^n mod m	O(n)	O(log n) — modular squaring
Modular inverse, prime m	—	O(log m) — Fermat
nth Fibonacci	O(n)	O(log n) — matrix power

a^n = (a^(n/2))² when n is even — halving gives O(log n)
The loop walks n's binary digits, squaring the base each step
Take the modulus after every multiplication to stay bounded
Matrix exponentiation extends it to any linear recurrence