Modular Exponentiation & Primality

The previous lesson built the arithmetic of $Z_{m}$ : addition, multiplication, and the modular inverse via the extended Euclidean algorithm. One operation remains: exponentiation — given $a$ , $n$ , and a modulus $m$ , compute $a^{n} mod m$ . The obvious loop multiplies $a$ into an accumulator $n$ times, which is $Θ (n)$ , far too slow when $n$ is a 1024-bit number, as it routinely is in cryptography. Repeated squaring does it in $O (log n)$ multiplications, and underlies primality testing, the modular inverse, and the public-key primitives that secure the internet.

Binary exponentiation: repeated squaring

The idea is to read $n$ in binary. Write $n = \sum_{i} b_{i} 2^{i}$ with $b_{i} \in {0, 1}$ . Then

a^{n} = a^{\sum_{i} b_{i} 2^{i}} = i : b_{i} = 1 \prod a^{2^{i}} .

The numbers $a^{2^{i}}$ are just the repeated squares of $a$ : each one is the square of the previous, since $a^{2^{i + 1}} = (a^{2^{i}})^{2}$ . So we sweep the bits of $n$ from least to most significant, keeping a running square $a^{2^{i}}$ , and whenever the current bit is $1$ we multiply that square into the result.

Square at each step; multiply the running square into the accumulator where the bit of

n

1

Each square is one multiplication and there are $⌊ log_{2} n ⌋ + 1$ of them; each bit costs at most one more multiplication into the accumulator. So the total is at most $2 ⌊ log_{2} n ⌋ + O (1)$ multiplications, which is $O (log n)$ .¹ Reducing modulo $m$ after every multiplication keeps every intermediate value below $m^{2}$ .

Algorithm:

\textsc{ModPow}(a, n, m)

— iterative bit-scan, returns

a^n \bmod m

1
$result \gets 1$
2
$a \gets a \bmod m$
3
while $n > 0$ do
4
if $n \bmod 2 = 1$ then
low bit set
5
$result \gets (result \cdot a) \bmod m$
6
$a \gets (a \cdot a) \bmod m$
square for next bit
7
$n \gets \lfloor n / 2 \rfloor$
8
return $result$

mod_pow.pypython

def mod_pow(base: int, exponent: int, modulus: int) -> int:
  """
    Compute `base ** exponent % modulus` in O(log exponent) multiplications\n
    by square-and-multiply. The exponent must be non-negative; the modulus\n
    must be positive. A modulus of 1 collapses everything to 0.\n
  """
  if modulus <= 0:
    raise ValueError("modulus must be positive")
  if exponent < 0:
    raise ValueError("exponent must be non-negative")

  # sweep the exponent's bits low-to-high, squaring as we go.
  result: int = 1 % modulus
  square: int = base % modulus
  remaining: int = exponent
  while remaining > 0:

    # a set low bit contributes the current square to the product.
    if remaining & 1:
      result = (result * square) % modulus

    # square for the next, more significant bit.
    square = (square * square) % modulus
    remaining >>= 1
  return result

def mod_pow_recursive(base: int, exponent: int, modulus: int) -> int:
  """
    Recursive halving variant of `mod_pow`: square the child's value, and\n
    fold in one extra factor of `base` when the exponent is odd.\n
  """
  if modulus <= 0:
    raise ValueError("modulus must be positive")
  if exponent < 0:
    raise ValueError("exponent must be non-negative")
  return _mod_pow_recursive(base % modulus, exponent, modulus)

def _mod_pow_recursive(base: int, exponent: int, modulus: int) -> int:
  """
    Core recursion assuming `base` is already reduced mod `modulus`.\n
  """
  if exponent == 0:
    return 1 % modulus

  # square the half-exponent result, then fold in one base for an odd power.
  half: int = _mod_pow_recursive(base, exponent // 2, modulus)
  half = (half * half) % modulus
  if exponent & 1:
    half = (half * base) % modulus
  return half

For example, run this on $3^{13} mod 7$ : the running square is $3, 2, 4, 2$ across the bits of $13 = 110 1_{2}$ , and the accumulator picks up a factor at each set bit, finishing at $3$ .

ModPow (3, 13, 7)

traced — bits of

13 = 110 1_{2}

drive square-and-multiply

The same computation has a clean recursive shape, splitting the exponent in half:

a^{n} = ⎩ ⎨ ⎧ 1 (a^{n /2})^{2} (a^{⌊ n /2 ⌋})^{2} \cdot a n = 0, n even, n odd .

Algorithm:

\textsc{ModPow-Rec}(a, n, m)

— recursive halving

1
if $n = 0$ then return $1$
2
$h \gets \textsc{ModPow-Rec}(a, \lfloor n/2 \rfloor, m)$
3
$h \gets (h \cdot h) \bmod m$
4
if $n \bmod 2 = 1$ then
5
$h \gets (h \cdot a) \bmod m$
6
return $h$

The recursion descends by halving the exponent and rebuilds the answer on the way back up: each return squares the child's value, and an odd exponent multiplies one extra copy of $a$ . For $a^{13}$ the four downward halvings $13 \to 6 \to 3 \to 1 \to 0$ unwind into four squarings, three of them carrying the extra $\cdot a$ from an odd level.

Recursive halving for

a^{13}

— each level squares, odd exponents multiply one extra

a

The doubling structure is not special to integers. Replace multiply with matrix multiply and the identity $[1110]^{n} = [F_{n + 1} F_{n} F_{n} F_{n - 1}]$ gives the $n$ -th Fibonacci number in $O (log n)$ matrix multiplications by the very same repeated-squaring loop.

fibonacci_matrix.pypython

from typing import Optional

# A 2x2 integer matrix as a flat (a, b, c, d) tuple = [[a, b], [c, d]].
Matrix = tuple[int, int, int, int]

_IDENTITY: Matrix = (1, 0, 0, 1)
_FIBONACCI_BASE: Matrix = (1, 1, 1, 0)

def _multiply(left: Matrix, right: Matrix, modulus: Optional[int]) -> Matrix:
  """
    Product of two 2x2 matrices, optionally reduced entrywise mod `modulus`.\n
  """
  # the four dot products of the row-by-column rule.
  left_a, left_b, left_c, left_d = left
  right_a, right_b, right_c, right_d = right
  product: Matrix = (
    left_a * right_a + left_b * right_c,
    left_a * right_b + left_b * right_d,
    left_c * right_a + left_d * right_c,
    left_c * right_b + left_d * right_d,
  )

  # reduce entrywise when a modulus keeps the entries bounded.
  if modulus is not None:
    return (
      product[0] % modulus,
      product[1] % modulus,
      product[2] % modulus,
      product[3] % modulus,
    )
  return product

def _matrix_power(matrix: Matrix, exponent: int, modulus: Optional[int]) -> Matrix:
  """
    `matrix` raised to `exponent` by repeated squaring — the same bit-scan\n
    as scalar `mod_pow`, with matrix multiply as the operation.\n
  """
  # accumulate into the identity, squaring the base each bit.
  result: Matrix = _IDENTITY
  square: Matrix = matrix
  remaining: int = exponent
  while remaining > 0:

    # a set low bit folds the current square into the running product.
    if remaining & 1:
      result = _multiply(result, square, modulus)

    # square for the next, more significant bit.
    square = _multiply(square, square, modulus)
    remaining >>= 1

  return result

def fibonacci(index: int, modulus: Optional[int] = None) -> int:
  """
    The Fibonacci number F(index), with F(0) = 0 and F(1) = 1, computed in\n
    O(log index) matrix multiplications. Pass `modulus` to reduce the result\n
    (and every intermediate) mod that value. The index must be non-negative.\n
  """
  if index < 0:
    raise ValueError("index must be non-negative")
  if index == 0:
    return 0 if modulus is None else 0 % modulus

  # [[1,1],[1,0]]^index has F(index) in its top-right (and bottom-left) entry.
  powered: Matrix = _matrix_power(_FIBONACCI_BASE, index, modulus)
  return powered[1]

Fermat's little theorem

Repeated squaring lets us compute large powers; number theory tells us what those powers are modulo a prime.

A corollary recovers the modular inverse from the previous lesson without the extended Euclidean algorithm: multiplying $a^{p - 1} \equiv 1$ by $a^{- 1}$ gives

a^{- 1} \equiv a^{p - 2} (mod p),

a single $ModPow$ call. This only works for a prime modulus, but that is exactly the common case in competitive programming, where arithmetic is done modulo a fixed prime such as $1 0^{9} + 7$ .² For a general modulus, Euler's theorem generalizes Fermat: if $gcd (a, m) = 1$ then $a^{φ (m)} \equiv 1 (mod m)$ , where $φ$ is Euler's totient, giving $a^{- 1} \equiv a^{φ (m) - 1}$ .

modular_inverse.pypython

from math import gcd

from mod_pow import mod_pow

def inverse_modulo_prime(value: int, prime: int) -> int:
  """
    The inverse of `value` modulo a prime `prime`, computed as\n
    `value ** (prime - 2) % prime` by Fermat's little theorem. Raises if\n
    `value` is a multiple of `prime` (it has no inverse there). The caller\n
    is responsible for `prime` actually being prime.\n
  """
  if prime < 2:
    raise ValueError("modulus must be a prime at least 2")
  if value % prime == 0:
    raise ValueError("value has no inverse: it is a multiple of the modulus")
  return mod_pow(value, prime - 2, prime)

def euler_totient(number: int) -> int:
  """
    Euler's totient phi(number): the count of integers in 1..number that are\n
    coprime to `number`. Computed from its prime factorization by trial\n
    division. Defined for positive integers.\n
  """
  if number < 1:
    raise ValueError("number must be positive")
  # peel prime factors out of `remaining`, scaling the count as we go.
  result: int = number
  remaining: int = number
  factor: int = 2
  while factor * factor <= remaining:
    if remaining % factor == 0:

      # each distinct prime factor p scales the count by (1 - 1/p).
      while remaining % factor == 0:
        remaining //= factor
      result -= result // factor
    factor += 1

  # a leftover above 1 is the final large prime factor.
  if remaining > 1:
    result -= result // remaining
  return result

def inverse_modulo(value: int, modulus: int) -> int:
  """
    The inverse of `value` modulo any `modulus` coprime to it, via Euler's\n
    theorem: a^-1 = a^(phi(modulus) - 1) (mod modulus). Raises when `value`\n
    and `modulus` share a common factor, since no inverse exists then.\n
  """
  if modulus < 1:
    raise ValueError("modulus must be positive")
  if gcd(value % modulus, modulus) != 1:
    raise ValueError("value and modulus must be coprime for an inverse")
  totient: int = euler_totient(modulus)
  return mod_pow(value, totient - 1, modulus)

Primality testing

How do we decide whether a number $n$ is prime? Three approaches, in increasing power.

Trial division — $O (n)$

If $n$ has a nontrivial factor it has one no larger than $n$ (factors come in pairs $d \cdot (n / d)$ , and the smaller is $\leq n$ ). So testing every candidate divisor up to $⌊ n ⌋$ settles the question.

Algorithm:

\textsc{IsPrime-Trial}(n)

—

O(\sqrt n)

deterministic test

1
if $n < 2$ then return false
2
$d \gets 2$
3
while $d \cdot d \le n$ do
4
if $n \bmod d = 0$ then return false
5
$d \gets d + 1$
6
return true

primality_trial_division.pypython

from typing import Optional

def is_prime_trial(number: int) -> bool:
  """
    Whether `number` is prime, by testing every divisor up to its square\n
    root. Numbers below 2 are not prime. Runs in O(sqrt number) time.\n
  """
  if number < 2:
    return False
  if number < 4:
    return True
  if number % 2 == 0:
    return False

  # only odd divisors remain once 2 is ruled out; step by 2.
  divisor: int = 3
  while divisor * divisor <= number:
    if number % divisor == 0:
      return False
    divisor += 2
  return True

def smallest_factor(number: int) -> Optional[int]:
  """
    The smallest prime factor of `number`, or None when `number` is prime\n
    (or below 2). The witness a trial-division proof of compositeness hands\n
    back, found in O(sqrt number) time.\n
  """
  if number < 2:
    return None
  if number % 2 == 0:
    return 2

  # the first odd divisor up to sqrt(number) is its smallest prime factor.
  divisor: int = 3
  while divisor * divisor <= number:
    if number % divisor == 0:
      return divisor
    divisor += 2
  return None

This is perfectly adequate for one moderate number (say $n \leq 1 0^{14}$ ), and is the right tool when a problem hands you a single value. It is hopeless for a 200-digit cryptographic number, where $n$ is astronomically large.

The Fermat test — probabilistic

Fermat's little theorem runs in reverse as a compositeness detector. If $n$ is prime, then $a^{n - 1} \equiv 1 (mod n)$ for every $a$ coprime to $n$ . So if we find a single witness $a$ with $a^{n - 1} \neq \equiv 1 (mod n)$ , then $n$ is certainly composite, and one $ModPow$ call refutes primality. If instead $a^{n - 1} \equiv 1$ , then $n$ is only probably prime; repeat with several random $a$ to raise confidence.

This asymmetry mirrors the soundness versus completeness distinction from the foundations. Read as a primality test (declare prime when $a^{n - 1} \equiv 1$ ), it is complete: Fermat's little theorem guarantees that every prime passes, so no prime is ever wrongly rejected. But it is not sound as a primality certifier: some composites pass too, so a prime verdict can be a false positive. Read as a compositeness test (declare composite when $a^{n - 1} \neq \equiv 1$ ), the verdicts flip roles: now it is sound (a failing base proves compositeness, so a composite verdict is never wrong) but incomplete (it can miss composites that happen to pass). Which property holds depends on which answer is trusted.

The test is unsound because of the Carmichael numbers: composites such as $561 = 3 \cdot 11 \cdot 17$ for which $a^{n - 1} \equiv 1 (mod n)$ holds for every $a$ coprime to $n$ . No choice of coprime witness exposes them, so the Fermat test declares them prime no matter how many rounds are run, and there are infinitely many of them. A stronger test is needed.

fermat_test.pypython

import random
from math import gcd

from mod_pow import mod_pow

def is_fermat_witness(number: int, base: int) -> bool:
  """
    Whether `base` witnesses that `number` is composite, i.e. whether\n
    base^(number-1) != 1 (mod number). A True verdict is airtight: a Fermat\n
    witness convicts `number` with certainty.\n
  """
  if number < 2:
    raise ValueError("number must be at least 2")
  return mod_pow(base, number - 1, number) != 1

def is_probably_prime_fermat(
  number: int,
  rounds: int = 20,
  rng: random.Random | None = None,
) -> bool:
  """
    Whether `number` survives `rounds` random Fermat bases. False means\n
    `number` is certainly composite; True means it is only probably prime,\n
    and is unconditionally wrong for Carmichael numbers, which pass for\n
    every coprime base.\n
  """
  # settle the small and even cases before any random work.
  if number < 2:
    return False
  if number < 4:
    return True
  if number % 2 == 0:
    return False

  # try fresh random bases, each a potential proof of compositeness.
  generator: random.Random = rng if rng is not None else random.Random()
  for _ in range(rounds):
    base: int = generator.randrange(2, number - 1)

    # a base sharing a factor with n is itself a proof of compositeness.
    if gcd(base, number) != 1:
      return False
    if is_fermat_witness(number, base):
      return False
  return True

The full square-root chain for $561$ shows what the Fermat test misses. With witness $a = 2$ and $560 = 2^{4} \cdot 35$ , the chain ends at $1$ , so the Fermat test, which checks only that last entry, passes $561$ . But the chain reaches $1$ from $67$ , a value that is neither $+ 1$ nor $- 1$ . A prime can never square a non- $\pm 1$ residue to $1$ , so this one extra observation proves $561$ composite.

The Carmichael number

561

passes Fermat (

a^{560} \equiv 1

) but fails Miller–Rabin: the chain hits

1

from

67

, a nontrivial square root of

1

Miller–Rabin

Miller–Rabin strengthens the Fermat test by exploiting a second fact about primes: in $Z_{p}$ , the only square roots of $1$ are $\pm 1$ . Write the even number $n - 1$ as

n - 1 = 2^{s} d, d odd .

For a witness $a$ , consider the chain obtained by computing $a^{d} mod n$ and then squaring it $s$ times:

a^{d}, a^{2 d}, a^{4 d}, \dots, a^{2^{s - 1} d}, a^{2^{s} d} = a^{n - 1} (mod n) .

If $n$ is prime, this chain must end at $1$ (Fermat), and the first time it reaches $1$ it must arrive from $- 1$ , because $1$ has no square root other than $\pm 1$ . So a prime forces one of two patterns: either $a^{d} \equiv 1$ , or some $a^{2^{r} d} \equiv - 1$ for $0 \leq r < s$ . If neither holds, we have found a $2^{r} d$ -step value that is a nontrivial square root of $1$ (it squares to $1$ but is not $\pm 1$ ), which a prime can never have, so $a$ is a witness that $n$ is composite.³

A nontrivial square root of

1

in the chain betrays a composite

Algorithm:

\textsc{Miller-Rabin}(n, a)

— true if

a

fails to witness compositeness

1
write $n - 1 = 2^{s} d$ with $d$ odd
2
$x \gets \textsc{ModPow}(a, d, n)$
3
if $x = 1$ or $x = n - 1$ then return true
probably prime
4
repeat $s - 1$ times
5
$x \gets (x \cdot x) \bmod n$
6
if $x = n - 1$ then return true
hit $-1$
7
return false
nontrivial $\sqrt 1$ ⇒ composite

miller_rabin.pypython

import random
from typing import Iterable

from mod_pow import mod_pow

# A proven-correct deterministic base set for every n < 3.317e24.
_DETERMINISTIC_BASES: tuple[int, ...] = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)

def _decompose(even_number: int) -> tuple[int, int]:
  """
    Factor `even_number` as 2^power * odd_part, returning (power, odd_part).\n
    Used to split n-1 into the exponent ladder Miller-Rabin walks.\n
  """
  power: int = 0
  odd_part: int = even_number
  while odd_part % 2 == 0:
    odd_part //= 2
    power += 1
  return power, odd_part

def passes_miller_rabin(number: int, base: int) -> bool:
  """
    Whether `base` fails to witness that `number` is composite — True when\n
    the square-root chain for this base looks prime. A False verdict is an\n
    airtight proof of compositeness. Assumes `number` is an odd integer\n
    greater than 2 and `base` lies in 2..number-2.\n
  """
  power, odd_part = _decompose(number - 1)
  value: int = mod_pow(base, odd_part, number)

  # a^d already 1, or already -1: the chain starts in a prime-compatible state.
  if value == 1 or value == number - 1:
    return True

  # square up the ladder; reaching -1 anywhere means this base looks prime.
  for _ in range(power - 1):
    value = (value * value) % number
    if value == number - 1:
      return True

  # never hit -1 before landing on 1: a nontrivial square root of 1.
  return False

def is_prime_deterministic(number: int) -> bool:
  """
    Deterministic primality for `number` using the fixed twelve-prime base\n
    set. Provably correct for every number that fits in 64 (indeed 80) bits,\n
    and exact for all inputs here since Python integers are unbounded only\n
    in principle — for n below 3.3 * 10^24 the answer is certified.\n
  """
  if number < 2:
    return False

  # a base equal to n proves prime; a base dividing n proves composite.
  for small_prime in _DETERMINISTIC_BASES:
    if number == small_prime:
      return True
    if number % small_prime == 0:
      return False

  # number is odd and shares no factor with any base; run the witness test.
  for base in _DETERMINISTIC_BASES:
    if not passes_miller_rabin(number, base):
      return False
  return True

def is_probably_prime(
  number: int,
  rounds: int = 20,
  rng: random.Random | None = None,
) -> bool:
  """
    Probabilistic Miller-Rabin over `rounds` random bases. False is certain\n
    compositeness; True leaves a false-"prime" probability below 4^-rounds.\n
    Unlike the Fermat test, Carmichael numbers are not immune.\n
  """
  # settle the small and even cases before any random work.
  if number < 2:
    return False
  if number < 4:
    return True
  if number % 2 == 0:
    return False

  # any base that fails the witness test convicts n as composite.
  generator: random.Random = rng if rng is not None else random.Random()
  for _ in range(rounds):
    base: int = generator.randrange(2, number - 1)
    if not passes_miller_rabin(number, base):
      return False
  return True

def is_prime_with_bases(number: int, bases: Iterable[int]) -> bool:
  """
    Run Miller-Rabin with a caller-supplied set of `bases` — the building\n
    block both the deterministic and probabilistic entry points share.\n
    Bases that are multiples of `number` are reduced and skipped.\n
  """
  # settle the small and even cases the witness test cannot handle.
  if number < 2:
    return False
  if number == 2 or number == 3:
    return True
  if number % 2 == 0:
    return False

  # trivial residues carry no information; test the rest as witnesses.
  for base in bases:
    reduced: int = base % number
    if reduced in (0, 1, number - 1):
      continue
    if not passes_miller_rabin(number, reduced):
      return False
  return True

Miller–Rabin remains a sound compositeness test (a witness still proves $n$ composite), but its error is now one-sided and bounded. With random witnesses, each composite $n$ is exposed by at least three quarters of the possible $a$ , so $k$ independent rounds leave a false-:qprime probability below $4^{- k}$ . Carmichael numbers are not immune, because the square-root check detects structure the Fermat test cannot. A prime verdict is still not a certificate — it reads probably prime — but the false-positive rate can be driven arbitrarily low by adding rounds, unlike the unbounded error the Fermat test suffers on Carmichael numbers. The test can also be made deterministic for bounded inputs: there is a fixed small set of bases that never errs below a threshold. Testing against the first twelve primes ${2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}$ is a proven-correct deterministic primality test for all $n < 3.3 \times 1 0^{24}$ , covering every 64-bit (indeed every 80-bit) integer with a dozen $ModPow$ calls.⁴

Miller–Rabin decides whether $n$ is prime but never produces a factor. Factoring a large composite is a separate, much harder problem; Pollard's rho finds a nontrivial factor in expected $O (n^{1/4})$ time using a cycle-detection trick on $x \mapsto x^{2} + c mod n$ , and is the standard tool for splitting numbers too big for trial division. (The next lesson handles factoring small numbers wholesale with a sieve.)

pollard_rho_factor.pypython

import random
from math import gcd

from miller_rabin import is_prime_deterministic

def pollard_rho(number: int, rng: random.Random | None = None) -> int:
  """
    A nontrivial factor of composite `number` (not necessarily prime),\n
    found by Floyd cycle detection on x -> x^2 + c (mod number). Returns\n
    `number` itself for primes and small inputs that cannot be split this\n
    way (2 and below). Restarts with a fresh constant on a degenerate run.\n
  """
  if number % 2 == 0:
    return 2
  if number < 2 or is_prime_deterministic(number):
    return number

  # retry with a fresh polynomial constant until a split sticks.
  generator: random.Random = rng if rng is not None else random.Random()
  while True:

    # seed both walkers at the same random point on the new curve.
    constant: int = generator.randrange(1, number)
    start: int = generator.randrange(2, number)
    slow: int = start
    fast: int = start
    divisor: int = 1

    while divisor == 1:

      # the tortoise advances one step, the hare two.
      slow = (slow * slow + constant) % number
      fast = (fast * fast + constant) % number
      fast = (fast * fast + constant) % number
      divisor = gcd(abs(slow - fast), number)

    # divisor == number is a failed run; try a new constant.
    if divisor != number:
      return divisor

def factorize(number: int, rng: random.Random | None = None) -> list[int]:
  """
    The full prime factorization of `number` as a sorted list with\n
    multiplicity, e.g. 360 -> [2, 2, 2, 3, 3, 5]. Splits composites with\n
    Pollard's rho and certifies the pieces with deterministic Miller-Rabin.\n
    Defined for integers at least 1; 1 factorizes to the empty list.\n
  """
  if number < 1:
    raise ValueError("number must be positive")
  if number == 1:
    return []

  # work a stack of pieces, peeling primes off and splitting composites.
  factors: list[int] = []
  pending: list[int] = [number]
  while pending:

    # a prime piece is a final factor; 1 contributes nothing.
    current: int = pending.pop()
    if current == 1:
      continue
    if is_prime_deterministic(current):
      factors.append(current)
      continue

    # split off a nontrivial factor and recurse on both halves.
    factor: int = pollard_rho(current, rng)
    pending.append(factor)
    pending.append(current // factor)

  factors.sort()
  return factors

Why this matters: cryptography

Fast modular exponentiation and fast primality testing together underlie public-key cryptography. RSA picks two large random primes (found by Miller–Rabin), and both encryption and decryption are a single modular exponentiation $m \mapsto m^{e} mod N$ . Diffie–Hellman key exchange computes $g^{x} mod p$ the same way. In every case the security rests on a power being easy to compute but its inverse (factoring $N$ , or the discrete logarithm) being infeasible, and $ModPow$ is what makes the easy direction $O (log n)$ .

For contrast with the Carmichael trace, run the same square-root chain on a genuine prime, $n = 97$ . Write $n - 1 = 96 = 2^{5} \cdot 3$ , so $s = 5$ , $d = 3$ ; take witness $a = 5$ . The chain is $a^{d} = 5^{3} = 125 \equiv 28 (mod 97)$ , then repeated squares $2 8^{2} \equiv 63$ , $6 3^{2} \equiv 96 \equiv - 1$ . It hits $- 1$ at the third link, so Miller–Rabin returns probably prime immediately — and because $- 1$ appeared, every later square is $+ 1$ , exactly the pattern a prime is forced into. A composite would either miss $\pm 1$ entirely or reach $1$ from a non- $\pm 1$ value, as $561$ did.

Miller-Rabin on the prime

97

with

a = 5

: the chain reaches

- 1

(here

96

), the pattern a prime must show; contrast

561

, which never does.

Primality, factoring, and the quantum frontier

Primality and factoring have sharply different complexity, and the modern results are worth knowing.

Primality is in P (AKS). For decades, is $n$ prime? had fast randomized answers (Miller–Rabin) and a fast answer assuming the Riemann hypothesis (deterministic Miller under GRH), but no unconditional deterministic polynomial algorithm was known. In 2002 Agrawal, Kayal, and Saxena settled it: the AKS primality test decides primality in deterministic polynomial time, $\tilde{O} (log^{6} n)$ in the original paper, placing PRIMES firmly in P.⁵ It is a theoretical landmark rather than a practical tool — deterministic Miller–Rabin with a fixed base set is far faster for the sizes anyone actually tests — but it closed a question open since antiquity.

Factoring is (believed) hard. No polynomial algorithm is known for the reverse problem of splitting a composite. Pollard's rho finds a factor in $O (n^{1/4})$ ; the quadratic sieve and the general number field sieve (GNFS) do far better for large inputs, GNFS running in sub-exponential $exp (O ((log n)^{1/3} (log log n)^{2/3}))$ time — still super-polynomial, the reason RSA moduli of 2048+ bits remain secure.⁶ The current public factoring record (RSA-250, an 829-bit number, 2020) took thousands of CPU-core-years.

Shor's algorithm. The asymmetry that protects RSA does not survive quantum computation. In 1994 Peter Shor gave a quantum algorithm that factors an $n$ -bit integer in $O (n^{2} log n log log n)$ time by reducing factoring to period-finding and using the quantum Fourier transform.⁷ A large fault-tolerant quantum computer would break RSA and Diffie–Hellman outright, which is the entire motivation for post-quantum cryptography (lattice- and code-based schemes) now being standardized. The Fast Fourier Transform that appears later in this module is the classical analogue of the quantum Fourier transform in Shor's algorithm.

Takeaways

Binary exponentiation computes $a^{n} mod m$ in $O (log n)$ multiplications by repeated squaring, reading the bits of $n$ and multiplying in each square whose bit is $1$ ; reduce mod $m$ every step, and guard against overflow with 128-bit or $MulMod$ arithmetic. The same doubling gives $O (log n)$ Fibonacci via matrix powers.
Fermat's little theorem ( $a^{p - 1} \equiv 1 (mod p)$ ) gives the modular inverse $a^{- 1} \equiv a^{p - 2}$ for a prime modulus; Euler's theorem generalizes it to $a^{φ (m)} \equiv 1$ for any coprime $a$ .
Trial division tests divisors up to $n$ in $O (n)$ time, deterministic, fine for one moderate number.
The Fermat test detects composites probabilistically but is fooled by Carmichael numbers for every coprime witness.
Miller–Rabin writes $n - 1 = 2^{s} d$ and watches the square-root chain $a^{d}, a^{2 d}, \dots$ collapse to $1$ , catching the nontrivial square roots that betray a composite; it is probabilistic with random witnesses and deterministic below $3.3 \times 1 0^{24}$ with the first twelve primes as bases.
Modular exponentiation and primality testing power RSA and Diffie–Hellman; Pollard's rho handles factoring of large numbers when a factor is actually needed.

CLRS, Ch. 31 — Number-Theoretic Algorithms (§31.6): modular exponentiation by repeated squaring in $O (log n)$ multiplications, reducing mod $m$ at each step. ↩
Skiena, § — Number Theory: Fermat's little theorem and modular inverse via $a^{p - 2}$ for a prime modulus. ↩
CLRS, Ch. 31 — Number-Theoretic Algorithms (§31.8): the Miller–Rabin witness test built on nontrivial square roots of $1$ , with error below $2^{- 2 k}$ over $k$ rounds. ↩
Skiena, § — Number Theory: deterministic Miller–Rabin with a fixed small base set, and Pollard's rho for factoring. ↩
M. Agrawal, N. Kayal, N. Saxena, PRIMES is in P, Annals of Mathematics 160(2), 2004 (announced 2002): the first unconditional deterministic polynomial-time primality test. ↩
A. K. Lenstra, H. W. Lenstra Jr. (eds.), The Development of the Number Field Sieve, Springer LNM 1554, 1993; and Pomerance, A tale of two sieves, Notices AMS 43(12), 1996, for the quadratic sieve and GNFS running times. ↩
P. W. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM J. Computing 26(5), 1997 (conference version 1994). ↩

Binary exponentiation: repeated squaring

Fermat's little theorem

Primality testing

Trial division — O(n​)

The Fermat test — probabilistic

Miller–Rabin

Why this matters: cryptography

Primality, factoring, and the quantum frontier

Takeaways

Footnotes

Trial division — $O (n)$