Modular exponentiation

The following is an example in pseudocode based on Applied Cryptography by Bruce Schneier.[2] The inputs base, exponent, and modulus correspond to b, e, and m in the equations given above.

function modular_pow(base, exponent, modulus) is
    if modulus = 1 then
        return 0
    Assert :: (modulus - 1) * (modulus - 1) does not overflow base
    result := 1
    base := base mod modulus
    while exponent > 0 do
        if (exponent mod 2 == 1) then
            result := (result * base) mod modulus
        exponent := exponent >> 1
        base := (base * base) mod modulus
    return result

Note that upon entering the loop for the first time, the code variable base is equivalent to b. However, the repeated squaring in the third line of code ensures that at the completion of every loop, the variable base is equivalent to b²ⁱ mod m, where i is the number of times the loop has been iterated. (This makes i the next working bit of the binary exponent exponent, where the least-significant bit is exponent₀).

The first line of code simply carries out the multiplication in ${\displaystyle \prod _{i=0}^{n-1}b^{a_{i}2^{i}}{\pmod {m}}}$. If a is zero, no code executes since this effectively multiplies the running total by one. If a instead is one, the variable base (containing the value b²ⁱ mod m of the original base) is simply multiplied in.

In this example, the base b is raised to the exponent e = 13. The exponent is 1101 in binary. There are four binary digits, so the loop executes four times, with values a₀ = 1, a₁ = 0, a₂ = 1, and a₃ = 1.

First, initialize the result ${\displaystyle R}$ to 1 and preserve the value of b in the variable x:

${\displaystyle R\leftarrow 1\,(=b^{0}){\text{ and }}x\leftarrow b}$.

Step 1) bit 1 is 1, so set ${\displaystyle R\leftarrow R\cdot x{\text{ }}(=b^{1})}$;

set ${\displaystyle x\leftarrow x^{2}{\text{ }}(=b^{2})}$.

Step 2) bit 2 is 0, so do not reset R;

set ${\displaystyle x\leftarrow x^{2}{\text{ }}(=b^{4})}$.

Step 3) bit 3 is 1, so set ${\displaystyle R\leftarrow R\cdot x{\text{ }}(=b^{5})}$;

set ${\displaystyle x\leftarrow x^{2}{\text{ }}(=b^{8})}$.

Step 4) bit 4 is 1, so set ${\displaystyle R\leftarrow R\cdot x{\text{ }}(=b^{13})}$;

This is the last step so we don't need to square x.

We are done: R is now ${\displaystyle b^{13}}$.

Here is the above calculation, where we compute b = 4 to the power e = 13, performed modulo 497.

Initialize:

${\displaystyle R\leftarrow 1\,(=b^{0})}$ and ${\displaystyle x\leftarrow b=4}$.

Step 1) bit 1 is 1, so set ${\displaystyle R\leftarrow R\cdot 4\equiv 4{\pmod {497}}}$;

set ${\displaystyle x\leftarrow x^{2}{\text{ }}(=b^{2})\equiv 4^{2}\equiv 16{\pmod {497}}}$.

Step 2) bit 2 is 0, so do not reset R;

set ${\displaystyle x\leftarrow x^{2}{\text{ }}(=b^{4})\equiv 16^{2}\equiv 256{\pmod {497}}}$.

Step 3) bit 3 is 1, so set ${\displaystyle R\leftarrow R\cdot x{\text{ }}(=b^{5})\equiv 4\cdot 256\equiv 30{\pmod {497}}}$;

set ${\displaystyle x\leftarrow x^{2}{\text{ }}(=b^{8})\equiv 256^{2}\equiv 429{\pmod {497}}}$.

Step 4) bit 4 is 1, so set ${\displaystyle R\leftarrow R\cdot x{\text{ }}(=b^{13})\equiv 30\cdot 429\equiv 445{\pmod {497}}}$;

We are done: R is now ${\displaystyle 4^{13}\equiv 445{\pmod {497}}}$, the same result obtained in the previous algorithms.

The running time of this algorithm is O(log exponent). When working with large values of exponent, this offers a substantial speed benefit over the previous two algorithms, whose time is O(exponent). For example, if the exponent was 2²⁰ = 1048576, this algorithm would have 20 steps instead of 1048576 steps.

function modPow(b, e, m)
  if m == 1 then
    return 0
  else
    local r = 1
    b = b % m
    while e > 0 do
      if e % 2 == 1 then
        r = (r*b) % m
      end
      b = (b*b) % m
      e = e >> 1     --use 'e = math.floor(e / 2)' on Lua 5.2 or older
    end
    return r
  end
end

In The Art of Computer Programming, Vol. 2, Seminumerical Algorithms, page 463, Donald Knuth notes that contrary to some assertions, this method does not always give the minimum possible number of multiplications. The smallest counterexample is for a power of 15, when the binary method needs six multiplications. Instead, form x³ in two multiplications, then x⁶ by squaring x³, then x¹² by squaring x⁶, and finally x¹⁵ by multiplying x¹² and x³, thereby achieving the desired result with only five multiplications. However, many pages follow describing how such sequences might be contrived in general.

The m-th term of any constant-recursive sequence (such as Fibonacci numbers or Perrin numbers) where each term is a linear function of k previous terms can be computed efficiently modulo n by computing A^m mod n, where A is the corresponding k×k companion matrix. The above methods adapt easily to this application. This can be used for primality testing of large numbers n, for example.

Pseudocode

A recursive algorithm for ModExp(A, b, c) = A^b mod c, where A is a square matrix.

function Matrix_ModExp(Matrix A, int b, int c) is
    if b == 0 then
        return I  // The identity matrix
    if (b mod 2 == 1) then
        return (A * Matrix_ModExp(A, b - 1, c)) mod c
    Matrix D := Matrix_ModExp(A, b / 2, c)
    return (D * D) mod c

Diffie–Hellman key exchange uses exponentiation in finite cyclic groups. The above methods for modular matrix exponentiation clearly extend to this context. The modular matrix multiplication C ≡ AB (mod n) is simply replaced everywhere by the group multiplication c = ab.

In quantum computing, modular exponentiation appears as the bottleneck of Shor's algorithm, where it must be computed by a circuit consisting of reversible gates, which can be further broken down into quantum gates appropriate for a specific physical device. Furthermore, in Shor's algorithm it is possible to know the base and the modulus of exponentiation at every call, which enables various circuit optimizations.[3]