Fast exponentiation

Compute \(a^n\) with about \(\log_2 n\) multiplications instead of \(n - 1\), by squaring: \(a, a^2, a^4, a^8, \dots\) and multiplying together the ones that the binary digits of \(n\) ask for. With remainders taken at each step, it computes \(a^n \bmod m\) for numbers hundreds of digits long.

Time	Extra space	Output
\(O(\log n)\) multiplications	\(O(1)\)	\(a^n\), \(a^n \bmod m\), matrix powers

The idea

\(3^{13}\) the slow way is \(3 \cdot 3 \cdot \ldots \cdot 3\), twelve multiplications. But \(13 = 8 + 4 + 1\) (binary 1101), so

\[3^{13} = 3^{8} \cdot 3^{4} \cdot 3^{1},\]

and \(3^2, 3^4, 3^8\) come from squaring: \(3^2 = 3 \cdot 3\), \(3^4 = 3^2 \cdot 3^2\), \(3^8 = 3^4 \cdot 3^4\). That is 3 squarings and 2 more multiplications: 5 instead of 12. For \(n\) around a billion it is about 60 instead of a billion.

The same idea works for anything you can multiply: numbers modulo \(m\) (reduce after every multiplication so nothing grows), matrices (which gives Fibonacci numbers in \(O(\log n)\)), even polynomials.

How it works

1. result = 1.
2. While \(n > 0\): if the lowest binary digit of \(n\) is 1, result *= a. Then a *= a and \(n \leftarrow \lfloor n / 2 \rfloor\).
3. For \(a^n \bmod m\): take % m after each multiplication.

Code

PythonJavaC++

def power(a, n):
    """a to the power n with O(log n) multiplications (square and multiply)."""
    result = 1
    while n > 0:
        if n & 1:                            # this binary digit of n is 1: include the current square
            result *= a
        a *= a                               # a, a^2, a^4, a^8, ...
        n >>= 1
    return result


def power_mod(a, n, m):
    """a^n mod m, keeping every number below m."""
    result, a = 1 % m, a % m
    while n > 0:
        if n & 1:
            result = result * a % m
        a = a * a % m
        n >>= 1
    return result


def mat_mult(X, Y):
    return [[X[0][0] * Y[0][0] + X[0][1] * Y[1][0], X[0][0] * Y[0][1] + X[0][1] * Y[1][1]],
            [X[1][0] * Y[0][0] + X[1][1] * Y[1][0], X[1][0] * Y[0][1] + X[1][1] * Y[1][1]]]


def fib(n):
    """The nth Fibonacci number from the matrix [[1, 1], [1, 0]] to the power n."""
    result, M = [[1, 0], [0, 1]], [[1, 1], [1, 0]]
    while n > 0:
        if n & 1:
            result = mat_mult(result, M)
        M = mat_mult(M, M)
        n >>= 1
    return result[0][1]


if __name__ == "__main__":
    print(power(2, 10), power(3, 13), power(7, 0))
    MOD = 1_000_000_007
    print(power_mod(2, 1_000_000, MOD), power_mod(3, MOD - 2, MOD) * 3 % MOD)
    print(power_mod(5, 117, 19), fib(10), fib(90))

public class FastPower {
    /** a to the power n with O(log n) multiplications (square and multiply). */
    static long power(long a, long n) {
        long result = 1;
        while (n > 0) {
            if ((n & 1) == 1) result *= a;            // this binary digit of n is 1: include the current square
            a *= a;                                   // a, a^2, a^4, a^8, ...
            n >>= 1;
        }
        return result;
    }

    /** a^n mod m, keeping every number below m. */
    static long powerMod(long a, long n, long m) {
        long result = 1 % m;
        a %= m;
        while (n > 0) {
            if ((n & 1) == 1) result = result * a % m;
            a = a * a % m;
            n >>= 1;
        }
        return result;
    }

    static long[][] matMult(long[][] X, long[][] Y) {
        return new long[][]{
            {X[0][0] * Y[0][0] + X[0][1] * Y[1][0], X[0][0] * Y[0][1] + X[0][1] * Y[1][1]},
            {X[1][0] * Y[0][0] + X[1][1] * Y[1][0], X[1][0] * Y[0][1] + X[1][1] * Y[1][1]}};
    }

    /** The nth Fibonacci number from the matrix [[1, 1], [1, 0]] to the power n. */
    static long fib(int n) {
        long[][] result = {{1, 0}, {0, 1}}, M = {{1, 1}, {1, 0}};
        while (n > 0) {
            if ((n & 1) == 1) result = matMult(result, M);
            n >>= 1;
            if (n > 0) M = matMult(M, M);             // skip the last squaring, which would overflow
        }
        return result[0][1];
    }

    public static void main(String[] args) {
        System.out.println(power(2, 10) + " " + power(3, 13) + " " + power(7, 0));
        final long MOD = 1_000_000_007L;
        System.out.println(powerMod(2, 1_000_000, MOD) + " " + powerMod(3, MOD - 2, MOD) * 3 % MOD);
        System.out.println(powerMod(5, 117, 19) + " " + fib(10) + " " + fib(90));
    }
}

#include <array>
#include <iostream>
using namespace std;

// a to the power n with O(log n) multiplications (square and multiply).
long long power(long long a, long long n) {
    long long result = 1;
    while (n > 0) {
        if (n & 1) result *= a;                      // this binary digit of n is 1: include the current square
        a *= a;                                      // a, a^2, a^4, a^8, ...
        n >>= 1;
    }
    return result;
}

// a^n mod m, keeping every number below m.
long long powerMod(long long a, long long n, long long m) {
    long long result = 1 % m;
    a %= m;
    while (n > 0) {
        if (n & 1) result = result * a % m;
        a = a * a % m;
        n >>= 1;
    }
    return result;
}

using Mat = array<array<long long, 2>, 2>;

Mat matMult(const Mat& X, const Mat& Y) {
    return {{{X[0][0] * Y[0][0] + X[0][1] * Y[1][0], X[0][0] * Y[0][1] + X[0][1] * Y[1][1]},
             {X[1][0] * Y[0][0] + X[1][1] * Y[1][0], X[1][0] * Y[0][1] + X[1][1] * Y[1][1]}}};
}

// The nth Fibonacci number from the matrix [[1, 1], [1, 0]] to the power n.
long long fib(int n) {
    Mat result = {{{1, 0}, {0, 1}}}, M = {{{1, 1}, {1, 0}}};
    while (n > 0) {
        if (n & 1) result = matMult(result, M);
        n >>= 1;
        if (n > 0) M = matMult(M, M);                // skip the last squaring, which would overflow
    }
    return result[0][1];
}

int main() {
    cout << power(2, 10) << " " << power(3, 13) << " " << power(7, 0) << "\n";
    const long long MOD = 1000000007LL;
    cout << powerMod(2, 1000000, MOD) << " " << powerMod(3, MOD - 2, MOD) * 3 % MOD << "\n";
    cout << powerMod(5, 117, 19) << " " << fib(10) << " " << fib(90) << "\n";
}

Output:

• 1024 1594323 1: \(2^{10}\), \(3^{13}\) and \(7^0\).
• 235042059 1: \(2^{1{,}000{,}000} \bmod (10^9 + 7)\), a number with 301,030 digits reduced to 9 digits, in about 26 multiplications; then \(3 \cdot 3^{p-2} \bmod p = 1\), showing that \(3^{p-2}\) is the inverse of 3 modulo the prime \(p\).
• 1 55 2880067194370816120: \(5^{117} \bmod 19 = 1\), and the Fibonacci numbers \(F(10)\) and \(F(90)\) from matrix powers.

When to use it

• Cryptography: RSA, Diffie–Hellman and digital signatures all compute \(a^n \bmod m\) with numbers of 600 or more digits; see cryptography.akadian.ca.
• Modular inverses modulo a prime: \(a^{-1} \equiv a^{p-2} \pmod p\) (Fermat).
• Primality tests such as Fermat and Miller–Rabin, which check \(a^{n-1} \bmod n\).
• Linear recurrences (Fibonacci and others) and counting paths of length \(n\) in a graph, through matrix powers.
• Hashing: the base powers in Rabin-Karp.

Common mistakes

• Overflow in a * a % m: with \(m\) near \(10^9\) the product fits in 64 bits; with \(m\) near \(10^{18}\) it does not; use 128-bit multiplication (__int128 in C++, Math.multiplyHigh or BigInteger in Java).
• 1 % m: for \(m = 1\) the answer is 0, not 1; starting with result = 1 % m handles it.
• Negative exponents: they need an inverse, which exists only when \(\gcd(a, m) = 1\).
• Squaring once too many: without a modulus, the last a *= a can overflow even though it is never used; the Java and C++ fib skip it.

Going deeper

Each part below starts closed; open the ones you want.

Sample walkthrough: 3^13

\(13 = 1101_2\). Reading bits from the lowest:

Bit	Value	Current square	Result
0	1	\(3^1 = 3\)	\(1 \cdot 3 = 3\)
1	0	\(3^2 = 9\)	3 (skip)
2	1	\(3^4 = 81\)	\(3 \cdot 81 = 243\)
3	1	\(3^8 = 6561\)	\(243 \cdot 6561 = 1{,}594{,}323\)

Modular example, \(5^{117} \bmod 19\): by Fermat's little theorem \(5^{18} \equiv 1\), and \(117 = 6 \cdot 18 + 9\), so it equals \(5^9 \bmod 19\). Squares: \(5^2 = 25 \equiv 6\), \(5^4 \equiv 36 \equiv 17\), \(5^8 \equiv 289 \equiv 4\); then \(5^9 \equiv 4 \cdot 5 = 20 \equiv 1\).

The math: the number of multiplications, and Fermat

Count. Let \(n\) have \(k = \lfloor \log_2 n \rfloor + 1\) binary digits, \(s\) of them equal to 1. The loop does \(k - 1\) useful squarings and \(s\) multiplications into the result: at most \(2 \lfloor \log_2 n \rfloor + 1\). For \(n = 1{,}000{,}000\) (20 binary digits, 7 of them 1) that is 19 squarings and 7 multiplications.

Invariant. At the start of each loop, \(\text{result} \cdot a^n\) equals the original \(a_0^{n_0}\). If \(n\) is even, \((a^2)^{n/2} = a^n\); if odd, \(\text{result} \cdot a \cdot (a^2)^{(n-1)/2} = \text{result} \cdot a^n\). When \(n = 0\), the result is the answer.

Not always optimal. The fewest multiplications for \(a^n\) is the length of the shortest addition chain for \(n\); for \(n = 15\), square-and-multiply uses 6 but 5 suffice (\(a^2, a^3, a^6, a^{12}, a^{15}\)). Finding the shortest chain is hard in general, and square-and-multiply is within a factor 2 of it.

Fermat's little theorem. If \(p\) is prime and \(p \nmid a\), then \(a^{p-1} \equiv 1 \pmod p\). So \(a \cdot a^{p-2} \equiv 1\): the inverse of \(a\) is \(a^{p-2} \bmod p\). Euler's generalisation: \(a^{\varphi(m)} \equiv 1 \pmod m\) when \(\gcd(a, m) = 1\).

Watch it step by step

Orange: the bit being read; red: the square that gets multiplied in; faded: a square that is skipped.

Variations and related ideas

• Left-to-right binary method: read bits from the top: square, then multiply by \(a\) if the bit is 1; the same count, and the multiplier is always the small \(a\).
• Sliding windows (precompute \(a^1, a^3, \dots, a^{2^w - 1}\)) use fewer multiplications for long exponents, such as those in RSA.
• Montgomery ladder: the same steps whatever the bits, so the time does not leak the secret exponent (a defence against timing attacks).
• Montgomery multiplication and Barrett reduction make each % m fast for big numbers.

See also

• Euclid's GCD: inverses for any modulus
• Fibonacci: \(F(n)\) with matrix powers
• Rabin-Karp

References

• D. E. Knuth, The Art of Computer Programming, Vol. 2, 3rd ed., Section 4.6.3 (evaluation of powers).
• A. J. Menezes, P. C. van Oorschot and S. A. Vanstone, Handbook of Applied Cryptography, CRC Press, 1996, Chapter 14 (efficient implementation).