Offensive Security Researcher & Quantum Cryptography Analyst

Number-Theoretic Algorithms: A Comprehensive Guide with Examples

Abdul Wahab Junaid
Number-Theoretic Algorithms: A Comprehensive Guide

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:

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).

So gcd(48, 18) = 6.

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

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

Time Complexity

Space Complexity

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

Disadvantages

Applications

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;
}

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

Common Mistakes

Further Reading

Exit mobile version