Number Theory Algorithms in Python: GCD, Sieve, RSA

By Engr. Henry Akpan

Old Mathematics, Modern Computing

Euclid described his method for the greatest common divisor over 2,300 years ago, and Eratosthenes sieved for primes in ancient Alexandria. Today the same ideas protect every online bank transfer, run inside competitive-programming solutions, and generate the keys that secure websites. This tutorial covers number theory algorithms in Python: the Euclidean algorithm for GCD and LCM, the extended Euclidean algorithm and modular inverses, the Sieve of Eratosthenes, prime factorisation, fast modular exponentiation, primality testing with Fermat's test and Miller–Rabin, and a small working example of RSA encryption. It builds on the tutorials on modular arithmetic and Big O notation.

The Greatest Common Divisor: Euclid's Algorithm

The greatest common divisor (GCD), also called the highest common factor, of two numbers is the largest number that divides both. Euclid's algorithm is based on one fact: gcd(a, b) = gcd(b, a mod b), because any number dividing both a and b also divides the remainder when a is divided by b. Repeat until the remainder is 0. The last non-zero value is the GCD.

def gcd(a, b):
    while b:
        print(f"gcd({a}, {b})")
        a, b = b, a % b
    return a

print(gcd(252, 105))

Output:

gcd(252, 105)
gcd(105, 42)
gcd(42, 21)
21

The numbers shrink very quickly: the algorithm takes O(log(min(a, b))) steps, so even numbers with hundreds of digits are handled in a fraction of a second. Python provides it ready-made as math.gcd.

The lowest common multiple

Since gcd(a, b) × lcm(a, b) = a × b, the LCM follows directly from the GCD:

import math

def lcm(a, b):
    return a // math.gcd(a, b) * b         # divide first to keep numbers small

print(lcm(12, 18), math.lcm(4, 6, 10))
print(math.gcd(1071, 462))

Output:

36 60
21

A classic use: two buses leave a park together, one every 12 minutes and one every 18 minutes. They next leave together after lcm(12, 18) = 36 minutes.

The Extended Euclidean Algorithm

The extended algorithm finds not only g = gcd(a, b), but also whole numbers x and y such that:

a × x + b × y = g

This is called Bézout's identity. It works backwards through the steps of Euclid's algorithm.

def extended_gcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extended_gcd(b, a % b)
    return g, y1, x1 - (a // b) * y1

g, x, y = extended_gcd(240, 46)
print(g, x, y, 240 * x + 46 * y)

Output:

2 -9 47 2

Modular inverses

The modular inverse of a modulo m is a number x with a × x ≡ 1 (mod m). It exists exactly when gcd(a, m) = 1, and the extended algorithm finds it: if a × x + m × y = 1, then a × x ≡ 1 (mod m).

def mod_inverse(a, m):
    g, x, _ = extended_gcd(a, m)
    if g != 1:
        raise ValueError(f"{a} has no inverse modulo {m}")
    return x % m

print(mod_inverse(3, 11), (3 * 4) % 11)
print(pow(3, -1, 11))                       # Python 3.8+ does this built in
try:
    mod_inverse(6, 9)
except ValueError as e:
    print(e)

Output:

4 1
4
6 has no inverse modulo 9

Modular inverses let you "divide" in modular arithmetic, and they are essential in cryptography.

The Sieve of Eratosthenes

To find all the primes up to n, testing each number separately is wasteful. The sieve instead crosses out multiples:

  1. Write down all the numbers from 2 to n, assuming each one is prime.
  2. Take the smallest number not yet crossed out, p. It is prime. Cross out all its multiples, starting from p², because smaller multiples have already been crossed out by smaller primes.
  3. Repeat until p² is greater than n.
def sieve(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    p = 2
    while p * p <= n:
        if is_prime[p]:
            for multiple in range(p * p, n + 1, p):
                is_prime[multiple] = False
        p += 1
    return [i for i, prime in enumerate(is_prime) if prime]

print(sieve(50))
print(len(sieve(1_000_000)))

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
78498

The sieve finds all 78,498 primes below one million in well under a second. Its running time is O(n log log n), which is almost linear.

Prime Factorisation

To factorise n, divide by 2, then 3, then 5, and so on, as long as each one divides. Only divisors up to √n need to be tried: whatever remains at the end, if it is greater than 1, must itself be prime.

def factorise(n):
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1 if d == 2 else 2               # after 2, try odd numbers only
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

print(factorise(360), factorise(97), factorise(1001))

def num_divisors(n):
    count = 1
    for exponent in factorise(n).values():
        count *= exponent + 1
    return count

print(num_divisors(360))

Output:

{2: 3, 3: 2, 5: 1} {97: 1} {7: 1, 11: 1, 13: 1}
24

So 360 = 2³ × 3² × 5, and it has (3 + 1)(2 + 1)(1 + 1) = 24 divisors. Trial division is fast for numbers up to about 10¹². Factorising the enormous numbers used in cryptography, with hundreds of digits, is far beyond any known method, and that difficulty is what keeps RSA secure.

Fast Modular Exponentiation

Calculating 7^560 mod 561 by first computing 7^560 would produce a number with 474 digits. Exponentiation by squaring keeps every intermediate result small by reducing modulo m at each step, and it needs only about log₂(exponent) multiplications.

The idea is that a^13 = a^8 × a^4 × a^1, because 13 is 1101 in binary. Square repeatedly to get a, a², a⁴, a⁸, …, and multiply together the ones that match the 1-bits of the exponent.

def power_mod(base, exp, mod):
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:                    # lowest bit is 1: include this power
            result = result * base % mod
        base = base * base % mod       # square for the next bit
        exp >>= 1
    return result

print(power_mod(3, 13, 7), pow(3, 13, 7), 3 ** 13 % 7)
print(power_mod(7, 560, 561))
print(pow(2, 10**18, 1_000_000_007))

Output:

3 3 3
1
719476260

Python's built-in three-argument pow(base, exp, mod) uses the same method, and it handles an exponent of 10¹⁸ instantly.

Primality Testing

Fermat's test

Fermat's little theorem states that if p is prime and a is not a multiple of p, then a^(p−1) ≡ 1 (mod p). So if a^(n−1) mod n is not 1, n is definitely composite. If it is 1, n is probably prime, but not certainly. Some composite numbers, called Carmichael numbers, fool the test for almost every base. The smallest is 561 = 3 × 11 × 17, and the previous example showed that 7^560 ≡ 1 (mod 561).

The Miller–Rabin test

Miller–Rabin strengthens Fermat's test so that no composite number can fool it for most bases. Testing the first twelve primes as bases gives a certain answer for every number below 3.3 × 10²⁴, which covers anything that fits in 64 bits.

def is_probable_prime(n):
    if n < 2:
        return False
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    for p in small_primes:
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:                  # write n - 1 as d * 2^s, with d odd
        d //= 2
        s += 1
    for a in small_primes:
        x = pow(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               # a proves that n is composite
    return True

print([n for n in range(90, 110) if is_probable_prime(n)])
print(is_probable_prime(561), is_probable_prime(2**61 - 1), is_probable_prime(2**61 + 1))

Output:

[97, 101, 103, 107, 109]
False True False

Miller–Rabin correctly rejects the Carmichael number 561, and it confirms instantly that 2⁶¹ − 1, a 19-digit Mersenne prime, is prime. Trial division would need over a billion divisions for that.

Putting It Together: Toy RSA Encryption

RSA, published in 1977, is the public-key system behind secure websites, digital signatures and much of online banking. It uses every algorithm in this tutorial:

  1. Choose two primes, p and q, and compute n = p × q and φ = (p − 1)(q − 1).
  2. Choose a public exponent e with gcd(e, φ) = 1.
  3. Compute the private exponent d, the modular inverse of e modulo φ.
  4. To encrypt a message m (a number smaller than n), compute c = m^e mod n. To decrypt, compute m = c^d mod n.

The pair (n, e) is the public key, which anyone may know. d is the private key. Finding d from n and e needs φ, which needs the factors p and q, and for real key sizes (n with over 600 digits), no one knows how to factorise n in any practical time.

p, q = 61, 53                      # tiny primes, for illustration only
n = p * q
phi = (p - 1) * (q - 1)
e = 17
assert math.gcd(e, phi) == 1
d = mod_inverse(e, phi)
print("public key:", (n, e), " private key:", d)

message = 65
cipher = pow(message, e, n)
decrypted = pow(cipher, d, n)
print(cipher, decrypted)

text = "HI"
encrypted = [pow(ord(ch), e, n) for ch in text]
print(encrypted, "".join(chr(pow(c, d, n)) for c in encrypted))

Output:

public key: (3233, 17)  private key: 2753
2790 65
[3000, 1486] HI

Real systems use enormous primes, plus carefully designed padding schemes. Never write your own cryptography for real use. Use established libraries, such as Python's cryptography package.

Common Mistakes with Number Theory Algorithms

  • Computing a huge power and then taking the remainder, as in a ** b % m, which is extremely slow for large b. Use pow(a, b, m).
  • Testing divisors beyond √n in trial division or the sieve, which wastes time.
  • Trusting Fermat's test alone. Carmichael numbers fool it. Use Miller–Rabin.
  • Expecting a modular inverse to always exist. It needs gcd(a, m) = 1.
  • Using floating-point division (/) in number-theory code. Use integer division (//) and %, because floats lose precision for big numbers.

Number Theory Algorithms: Practice Exercises

  1. Use the sieve to find the sum of all the primes below 100.
  2. Find the smallest positive x with 7x ≡ 1 (mod 26).
  3. Find the last three digits of 2^1000.
  4. Two numbers have GCD 6 and LCM 180. One of them is 36. Find the other.
  5. Find the largest prime factor of 600851475143.

Sample solutions:

print(sum(sieve(99)))
print(mod_inverse(7, 26))
print(pow(2, 1000, 1000))
print(6 * 180 // 36)
print(max(factorise(600851475143)))

Output:

1060
15
376
30
6857

Check question 2: 7 × 15 = 105 = 4 × 26 + 1.

Conclusion

Number theory algorithms are short, fast and powerful. Euclid's algorithm, gcd(a, b) = gcd(b, a mod b), finds the GCD in logarithmic time, and the LCM follows as ab/gcd. The extended version gives Bézout coefficients and modular inverses. The Sieve of Eratosthenes lists all the primes up to n almost linearly, and trial division up to √n factorises moderate numbers. Exponentiation by squaring, Python's pow(a, b, m), makes huge modular powers easy, enabling the Fermat and Miller–Rabin primality tests and RSA encryption, which rests on the difficulty of factorising large numbers.

Related Tutorials