Modular Arithmetic & Fast Exponentiation

Intuition, Vocabulary and Foundations

Modular arithmetic keeps only the remainder after division by a positive modulus. A familiar analogy is a clock: starting at hour $10$ and moving forward $5$ hours gives hour $3$, because $15$ and $3$ occupy the same position on a $12$-hour clock.

For an integer modulus $m\ge 1$, integers $a$ and $b$ are congruent modulo $m$, written

$$ a\equiv b\pmod m, $$

when $m$ divides their difference:

$$ m\mid(a-b). $$

Thus $17\equiv5\pmod{12}$ because $12\mid(17-5)$. Congruence is an equivalence relation: it is reflexive, symmetric, and transitive. It partitions the integers into residue classes. Modulo $5$, every integer belongs to one of the five classes represented by $0,1,2,3,4$.

The canonical residue of $a$ modulo $m$ is the unique integer $r$ satisfying

$$ a\equiv r\pmod m \qquad\text{and}\qquad 0\le r\lt m. $$

Consequently, the canonical residue of $-17$ modulo $5$ is $3$, since $-17=5(-4)+3$. Programming languages do not all assign the same sign to a remainder, so portable code should normalize explicitly:

$$ \operatorname{norm}_m(a)=((a\bmod m)+m)\bmod m. $$

Congruence is preserved by addition and multiplication. These replacement rules permit intermediate values to be reduced without changing the final residue.

A modular inverse of $a$ modulo $m$ is a value $x$ satisfying

$$ ax\equiv1\pmod m. $$

It exists exactly when $\gcd(a,m)=1$. For example, $3^{-1}\equiv4\pmod{11}$ because $3\cdot4\equiv1\pmod{11}$. In contrast, $6$ has no inverse modulo $9$ because $\gcd(6,9)=3$.

Multiplication of residues forms a monoid: it is associative and has identity $1\bmod m$, but not every residue is invertible. The invertible residues form the group $(\mathbb Z/m\mathbb Z)^\times$. If $\gcd(a,m)=1$, Euler's theorem states $a^{\varphi(m)}\equiv1\pmod m$. For prime $p$ and $p\nmid a$, Fermat's little theorem gives $a^{p-1}\equiv1\pmod p$.

Step-by-Step Derivations and Algorithms

A direct computation of $a^e$ performs $e-1$ multiplications and creates an enormous intermediate integer. Binary exponentiation instead uses the binary expansion

$$ e=\sum_{i=0}^{k}b_i2^i, \qquad b_i\in\{0,1\}, $$

so that

$$ a^e=\prod_{i:b_i=1}a^{2^i}. $$

Successive squaring generates $a,a^2,a^4,a^8,\ldots$. At each iteration, the lowest exponent bit decides whether the current power contributes to the result.

For variables result, base, and remaining, the loop invariant is

$$ \text{result}\cdot\text{base}^{\text{remaining}} \equiv a^e\pmod m. $$

Initially result is the multiplicative identity and remaining is $e$. When remaining is odd, multiplying result by base accounts for its lowest binary bit. Squaring base and halving remaining preserves the invariant. At termination, the invariant yields the answer. The algorithm performs $O(\log e)$ modular multiplications with $O(1)$ auxiliary storage.

In a fixed-width integer type, evaluating $(ab)\bmod m$ may overflow before reduction. Double-and-add multiplication uses the binary expansion of $b$ to conditionally add successive doublings of $a$. This costs $O(\log b)$ modular additions. Wider intermediate types, Barrett reduction, and Montgomery multiplication can provide faster workload-specific alternatives.

The extended Euclidean algorithm constructs integers $x$ and $y$ satisfying

$$ ax+my=\gcd(a,m). $$

When the gcd is $1$, reduction modulo $m$ gives $ax\equiv1\pmod m$, so the normalized coefficient $x$ is the inverse of $a$. This works for prime and composite moduli alike.

Exponent reduction is optional. If $\gcd(a,m)=1$, Euler's theorem permits $a^e\equiv a^{e\bmod\varphi(m)}\pmod m$. For prime $p$, Fermat permits reduction modulo $p-1$ when $p\nmid a$. Without the relevant coprimality condition, blindly reducing the exponent can change the answer.

flowchart TD
    A["Normalize operands and validate modulus"] --> B{"Requested operation"}
    B -->|Power| C{"Exponent nonnegative?"}
    B -->|Inverse or division| D["Compute extended gcd"]
    B -->|Product| E{"Native product overflow-safe?"}
    C -->|Yes| F["Binary exponentiation"]
    C -->|No| D
    D --> G{"gcd equals one?"}
    G -->|Yes| H["Normalize Bezout coefficient"]
    G -->|No| I["Report no inverse"]
    H --> F
    E -->|Yes| J["Multiply and reduce"]
    E -->|No| K["Double-and-add or wide reduction"]

Implementation, Complexity and Engineering

The public contract should require a positive modulus and return canonical residues. Modulus $1$ deserves explicit treatment: every nonnegative power and product is congruent to $0$, but no inverse exists. A negative exponent means inversion followed by ordinary exponentiation, so it must fail when the base is not a unit.

Python integers do not overflow, but the double-and-add implementation below models the technique required when a fixed-width product cannot hold $ab$. For high-throughput native code, a wider intermediate type or Montgomery reduction is usually faster.

def normalize(value: int, modulus: int) -> int:
    if modulus <= 0:
        raise ValueError("modulus must be positive")
    return value % modulus


def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    old_r, r = abs(a), abs(b)
    old_s, s = 1, 0
    old_t, t = 0, 1

    while r:
        quotient = old_r // r
        old_r, r = r, old_r - quotient * r
        old_s, s = s, old_s - quotient * s
        old_t, t = t, old_t - quotient * t

    x = -old_s if a < 0 else old_s
    y = -old_t if b < 0 else old_t
    return old_r, x, y


def modular_inverse(value: int, modulus: int) -> int:
    if modulus <= 1:
        raise ValueError("an inverse requires modulus greater than one")
    value = normalize(value, modulus)
    gcd, coefficient, _ = extended_gcd(value, modulus)
    if gcd != 1:
        raise ValueError("inverse does not exist")
    return coefficient % modulus


def multiply_mod(a: int, b: int, modulus: int) -> int:
    """Overflow-safe double-and-add multiplication."""
    if modulus <= 0:
        raise ValueError("modulus must be positive")
    a = normalize(a, modulus)
    b = normalize(b, modulus)
    result = 0

    while b:
        if b & 1:
            result = (result + a) % modulus
        a = (a + a) % modulus
        b >>= 1
    return result


def power_mod(base: int, exponent: int, modulus: int) -> int:
    """Compute base raised to exponent modulo modulus."""
    if modulus <= 0:
        raise ValueError("modulus must be positive")
    if exponent < 0:
        base = modular_inverse(base, modulus)
        exponent = -exponent

    base = normalize(base, modulus)
    result = 1 % modulus
    while exponent:
        if exponent & 1:
            result = multiply_mod(result, base, modulus)
        base = multiply_mod(base, base, modulus)
        exponent >>= 1
    return result


def divide_mod(numerator: int, denominator: int, modulus: int) -> int:
    inverse = modular_inverse(denominator, modulus)
    return multiply_mod(numerator, inverse, modulus)


assert normalize(-17, 5) == 3
assert extended_gcd(30, 18)[0] == 6
assert modular_inverse(3, 11) == 4
assert multiply_mod(-7, 13, 5) == 4
assert power_mod(3, 13, 17) == 12
assert power_mod(3, -1, 11) == 4
assert power_mod(5, 0, 1) == 0
assert divide_mod(8, 3, 11) == 10

large_a = 9_223_372_036_854_775_123
large_b = 8_765_432_109_876_543_211
large_m = 9_223_372_036_854_775_583
assert multiply_mod(large_a, large_b, large_m) == (large_a * large_b) % large_m

for operation in (
    lambda: normalize(5, 0),
    lambda: power_mod(2, 10, -7),
    lambda: modular_inverse(6, 9),
    lambda: modular_inverse(0, 11),
    lambda: modular_inverse(1, 1),
    lambda: power_mod(6, -1, 9),
):
    try:
        operation()
        raise AssertionError("expected ValueError")
    except ValueError:
        pass

print("all modular-arithmetic tests passed")

Applications, Case Studies and Exercises

Modular arithmetic appears in cyclic schedules, checksums, rolling hashes, combinatorics, cryptography, pseudorandom generators, and algorithms whose exact counts are too large to store. The implementation strategy depends on the number of operations, the available product width, whether inverses are needed, and whether the modulus is known to be prime.

A useful testing oracle compares optimized operations with direct arbitrary-precision arithmetic on small random inputs. Important properties include normalization bounds, agreement with a trusted power function for nonnegative exponents, the identity $a\cdot a^{-1}\equiv1\pmod m$, and preservation of congruence under addition and multiplication.