Number-Theoretic Algorithms: A Comprehensive Guide with Examples

I want to introduce number-theoretic algorithms as the collection of computational techniques I use to solve problems rooted in the properties of integers — things like divisibility, primality, modular arithmetic, and greatest common divisors. I consider these algorithms foundational because they underlie so much of modern computing, especially cryptography, hashing, and error-correcting codes. In this guide I focus mainly on the Euclidean algorithm for GCD, the extended Euclidean algorithm, and modular exponentiation, since I see these as the core building blocks I return to again and again.

History and Background

I trace the Euclidean algorithm back more than two thousand years to Euclid’s Elements (Book VII), making it one of the oldest algorithms I know of that is still in everyday use. The extended version, which also produces the coefficients of Bézout’s identity, was developed much later as number theory matured, with contributions I associate with mathematicians such as Étienne Bézout in the 18th century. Modular exponentiation, particularly the “fast” or “binary” method, became especially important in the 20th century once public-key cryptography (RSA, Diffie-Hellman) was invented by Rivest, Shamir, Adleman, and Diffie-Hellman in the 1970s, because these cryptosystems require efficiently raising large numbers to large powers modulo another large number.

Problem Statement

I group several closely related problems under this topic:

1. GCD problem: given two integers a and b, find their greatest common divisor.
2. Bézout’s identity problem: find integers x and y such that ax + by = gcd(a, b).
3. Modular exponentiation problem: given a, b, and n, compute a^b mod n efficiently, without computing the (potentially astronomically large) value a^b directly.
4. Primality testing problem: determine whether a given integer is prime.

Core Concepts

I rely on these definitions throughout:

• Divisibility: I say a divides b (written a | b) if there exists an integer k such that b = a·k.
• GCD: the greatest common divisor of a and b is the largest integer that divides both.
• Modular arithmetic: arithmetic performed with a “wraparound” at a fixed modulus n, where I only care about remainders after division by n.
• Congruence: I write a ≡ b (mod n) when n divides a - b.
• Coprime integers: two integers are coprime if their GCD is 1.
• Multiplicative inverse: an integer x such that a·x ≡ 1 (mod n), which exists precisely when a and n are coprime.

How It Works

Euclidean Algorithm (GCD):

1. I take two integers a and b, with a ≥ b ≥ 0.
2. If b = 0, I return a as the GCD.
3. Otherwise, I replace (a, b) with (b, a mod b) and repeat.

Extended Euclidean Algorithm:

1. I run the same recursive process as above.
2. At each step, I also track coefficients so that, once I reach the base case, I can work back up the recursion to express the GCD as ax + by.

Modular Exponentiation (fast/binary method):

1. I write the exponent b in binary.
2. I repeatedly square the base modulo n, and I multiply the running result by the current squared value whenever the corresponding bit of b is 1.
3. I take the result modulo n at every multiplication to keep the numbers small.

Working Principle

I find the Euclidean algorithm’s correctness rests on a simple invariant: gcd(a, b) = gcd(b, a mod b). Because I always reduce the pair (a, b) to a strictly smaller pair, the algorithm is guaranteed to terminate, and it does so in a number of steps proportional to the number of digits of the smaller number (related to the Fibonacci sequence in the worst case).

For modular exponentiation, my internal logic exploits the identity a^b = (a^(b/2))^2 (with an extra factor of a if b is odd), letting me compute the power in logarithmic time relative to the exponent instead of linear time, and I never let intermediate values grow beyond roughly n^2 in size because I reduce modulo n after every multiplication.

Mathematical Foundation

The Euclidean algorithm rests on the identity:

$$ \gcd(a, b) = \gcd(b, ; a \bmod b) $$

Bézout’s identity guarantees that for any integers a and b, there exist integers x and y such that:

$$ ax + by = \gcd(a, b) $$

For modular exponentiation, I use the recursive definition:

$$ a^b \bmod n = \begin{cases} 1 & \text{if } b = 0 \ \left( a^{b/2} \bmod n \right)^2 \bmod n & \text{if } b \text{ is even} \ a \cdot \left( a^{b-1} \bmod n \right) \bmod n & \text{if } b \text{ is odd} \end{cases} $$

This recursive halving gives me a running time of:

$$ O(\log b) $$

modular multiplications.

Diagrams

flowchart TD
    A["Start: gcd(a, b)"] --> B{"b == 0?"}
    B -- Yes --> C["Return a"]
    B -- No --> D["r = a mod b"]
    D --> E["a = b, b = r"]
    E --> B

Pseudocode

GCD(a, b)
    while b != 0
        r = a mod b
        a = b
        b = r
    return a

EXTENDED-GCD(a, b)
    if b == 0
        return (a, 1, 0)
    (g, x1, y1) = EXTENDED-GCD(b, a mod b)
    x = y1
    y = x1 - (a div b) * y1
    return (g, x, y)

MOD-EXP(a, b, n)
    result = 1
    a = a mod n
    while b > 0
        if b is odd
            result = (result * a) mod n
        b = b >> 1
        a = (a * a) mod n
    return result

Step-by-Step Example

GCD example: I compute gcd(48, 18).

• 48 mod 18 = 12 → pair becomes (18, 12)
• 18 mod 12 = 6 → pair becomes (12, 6)
• 12 mod 6 = 0 → pair becomes (6, 0)
• Since the second value is 0, I return 6.

So gcd(48, 18) = 6.

Modular exponentiation example: I compute 3^13 mod 7.

• I write 13 in binary as 1101.
• I initialize result = 1, a = 3 mod 7 = 3.
• Bit 1 (LSB, value 1): result = (1 * 3) mod 7 = 3; then a = 3*3 mod 7 = 2.
• Bit 0: I don’t multiply; a = 2*2 mod 7 = 4.
• Bit 1: result = (3 * 4) mod 7 = 12 mod 7 = 5; then a = 4*4 mod 7 = 2.
• Bit 1: result = (5 * 2) mod 7 = 10 mod 7 = 3; then a = 2*2 mod 7 = 4.

Final result: 3^13 mod 7 = 3. I can verify: 3^13 = 1594323, and 1594323 mod 7 = 3 — it checks out.

Time Complexity

• Euclidean algorithm: $O(\log(\min(a, b)))$ in all cases (best, average, worst), since each step reduces the values roughly by a factor related to the golden ratio; the true worst case occurs on consecutive Fibonacci numbers.
• Extended Euclidean algorithm: same asymptotic bound, $O(\log(\min(a, b)))$, since it performs the same recursive reduction with constant extra work per step.
• Modular exponentiation: $O(\log b)$ modular multiplications; if I also account for the cost of multiplying numbers of size up to n, the total bit-complexity is $O(\log b \cdot M(\log n))$, where M denotes the multiplication cost.

Space Complexity

• I need only $O(1)$ auxiliary space for the iterative Euclidean algorithm and iterative modular exponentiation.
• For the recursive extended Euclidean algorithm, I use $O(\log(\min(a,b)))$ space on the call stack, unless I convert it to an iterative version, in which case it also becomes $O(1)$.

Correctness Analysis

I justify the Euclidean algorithm’s correctness through the identity gcd(a, b) = gcd(b, a mod b), which I can prove because any common divisor of a and b must also divide a mod b, and conversely any common divisor of b and a mod b must divide a. Since the algorithm terminates when b = 0 (and gcd(a, 0) = a by definition), the returned value is correct by induction on the sequence of reductions.

For modular exponentiation, I justify correctness by induction on the exponent: the recursive definition exactly mirrors the mathematical identity for powers, and reducing modulo n at each multiplication does not change the final result modulo n, since modular arithmetic is compatible with multiplication ((x mod n)(y mod n) mod n = xy mod n).

Advantages

• I get extremely fast, efficient algorithms — logarithmic time for both GCD and modular exponentiation.
• These algorithms use very little memory.
• They are foundational for cryptographic systems I might build or study, such as RSA.
• They generalize cleanly to related problems, like modular inverses and Chinese Remainder Theorem computations.

Disadvantages

• Naive (non-fast) exponentiation, if I mistakenly use it instead of the binary method, becomes prohibitively slow for cryptographic-sized exponents.
• The extended Euclidean algorithm’s recursive form can be harder for me to reason about than the plain GCD algorithm.
• Working with very large integers (hundreds of digits, as in real cryptography) requires big-integer arithmetic libraries, since standard C integer types overflow quickly.

Applications

• Cryptography: RSA key generation and encryption/decryption rely directly on modular exponentiation and modular inverses.
• Hashing: many hash functions use modular arithmetic to keep values within a fixed range.
• Computer algebra systems: simplifying fractions requires the GCD.
• Error-correcting codes: number-theoretic structures underlie codes like Reed-Solomon.
• Competitive programming: GCD, modular inverse, and modular exponentiation appear constantly in algorithmic problem sets.

Implementation in C

#include <stdio.h>

/* I compute the GCD iteratively using the Euclidean algorithm */
long long gcd(long long a, long long b) {
    while (b != 0) {
        long long r = a % b;
        a = b;
        b = r;
    }
    return a;
}

/* I compute the extended GCD, filling in x and y such that a*x + b*y = gcd(a,b) */
long long extendedGcd(long long a, long long b, long long *x, long long *y) {
    if (b == 0) {
        *x = 1;
        *y = 0;
        return a;
    }
    long long x1, y1;
    long long g = extendedGcd(b, a % b, &x1, &y1);
    *x = y1;
    *y = x1 - (a / b) * y1;
    return g;
}

/* I compute (base^exp) mod m using fast binary exponentiation */
long long modExp(long long base, long long exp, long long mod) {
    long long result = 1;
    base = base % mod;

    while (exp > 0) {
        /* If the current bit is set, I fold base into the result */
        if (exp & 1) {
            result = (result * base) % mod;
        }
        exp = exp >> 1;
        base = (base * base) % mod;
    }

    return result;
}

int main(void) {
    long long a = 48, b = 18;
    printf("gcd(%lld, %lld) = %lld\n", a, b, gcd(a, b));

    long long x, y;
    long long g = extendedGcd(a, b, &x, &y);
    printf("extendedGcd: g=%lld, x=%lld, y=%lld (check: %lld*x + %lld*y = %lld)\n",
           g, x, y, a, b, a * x + b * y);

    long long base = 3, exp = 13, mod = 7;
    printf("%lld^%lld mod %lld = %lld\n", base, exp, mod, modExp(base, exp, mod));

    return 0;
}

#include <stdio.h>

/* I compute the GCD iteratively using the Euclidean algorithm */
long long gcd(long long a, long long b) {
    while (b != 0) {
        long long r = a % b;
        a = b;
        b = r;
    }
    return a;
}

/* I compute the extended GCD, filling in x and y such that a*x + b*y = gcd(a,b) */
long long extendedGcd(long long a, long long b, long long *x, long long *y) {
    if (b == 0) {
        *x = 1;
        *y = 0;
        return a;
    }
    long long x1, y1;
    long long g = extendedGcd(b, a % b, &x1, &y1);
    *x = y1;
    *y = x1 - (a / b) * y1;
    return g;
}

/* I compute (base^exp) mod m using fast binary exponentiation */
long long modExp(long long base, long long exp, long long mod) {
    long long result = 1;
    base = base % mod;

    while (exp > 0) {
        /* If the current bit is set, I fold base into the result */
        if (exp & 1) {
            result = (result * base) % mod;
        }
        exp = exp >> 1;
        base = (base * base) % mod;
    }

    return result;
}

int main(void) {
    long long a = 48, b = 18;
    printf("gcd(%lld, %lld) = %lld\n", a, b, gcd(a, b));

    long long x, y;
    long long g = extendedGcd(a, b, &x, &y);
    printf("extendedGcd: g=%lld, x=%lld, y=%lld (check: %lld*x + %lld*y = %lld)\n",
           g, x, y, a, b, a * x + b * y);

    long long base = 3, exp = 13, mod = 7;
    printf("%lld^%lld mod %lld = %lld\n", base, exp, mod, modExp(base, exp, mod));

    return 0;
}

I designed each function to be self-contained: gcd for the plain GCD, extendedGcd for Bézout coefficients, and modExp for fast modular exponentiation. I used long long throughout so the program can handle reasonably large numbers without overflowing 32-bit integers.

Sample Input and Output

Input: a = 48, b = 18, and for modular exponentiation base = 3, exp = 13, mod = 7.

Output:

gcd(48, 18) = 6
extendedGcd: g=6, x=-1, y=3 (check: 48*x + 18*y = 6)
3^13 mod 7 = 3

Optimization Techniques

• Binary GCD (Stein’s algorithm): I can replace division-based GCD with a version that uses only subtraction and bit-shifts, which is faster on hardware where division is expensive.
• Precomputed modular inverses: when I need many modular inverses under the same modulus, I precompute them in a batch using a linear-time technique rather than calling extended GCD repeatedly.
• Montgomery multiplication: for very large moduli, I can use Montgomery form to speed up repeated modular multiplications in modular exponentiation.
• Sliding-window exponentiation: I can process multiple bits of the exponent at once instead of one bit at a time, reducing the number of multiplications for large exponents.

Common Mistakes

• Forgetting to reduce base = base % mod before starting modular exponentiation, which can cause overflow.
• Using int instead of long long (or a big-integer type), causing silent overflow on large inputs.
• Swapping the order of x and y when reading off Bézout coefficients from the extended Euclidean recursion.
• Assuming gcd(a, 0) = 0 instead of the correct value, a.
• Implementing exponentiation with a simple loop that multiplies b times instead of using the fast binary method, leading to unnecessarily slow code for large exponents.

Further Reading

• Cormen, Leiserson, Rivest, Stein — Introduction to Algorithms, Chapter 31 (Number-Theoretic Algorithms): https://mitpress.mit.edu/9780262046305/introduction-to-algorithms/
• Rivest, Shamir, Adleman — “A Method for Obtaining Digital Signatures and Public-Key Cryptosystems,” Communications of the ACM, 1978: https://dl.acm.org/doi/10.1145/359340.359342
• Knuth — The Art of Computer Programming, Volume 2: Seminumerical Algorithms: https://www-cs-faculty.stanford.edu/~knuth/taocp.html
• NIST — Digital Signature Standard (FIPS 186-5), background on modular exponentiation in cryptography: https://csrc.nist.gov/pubs/fips/186-5/final