Binary exponentiation

Binary exponentiation computes a power such as $3^{13}$ 3^{13}, or its remainder after dividing by $1 0^{9} + 7$ 10^9 + 7, with at most two multiplications per binary digit of the exponent instead of one multiplication per factor. It is the standard tool when a problem asks for a power with an exponent up to $1 0^{9}$ 10^9 or $1 0^{18}$ 10^{18}.

Explanation

To compute $3^{8}$ 3^8, we multiply 3 by itself until there are 8 factors of 3. One factor at a time, that takes 7 multiplications. Squaring is faster: 3 · 3 = 9, then 9 · 9 = 81, then 81 · 81 = 6561. Each squaring doubles the number of factors, from 1 to 2 to 4 to 8, so 3 multiplications are enough. The same way, $3^{1024}$ 3^{1024} takes 10 squarings instead of 1023 multiplications.

Four boxes in a row holding 3, 9, 81 and 6561, labelled 3 to the power 1, 2, 4 and 8; each number is the square of the one before it.

Fig. 1. The squares of 3 up to $3^{8}$ 3^8.

Most exponents are not powers of 2. Take $3^{13}$ 3^{13}. The squares go 3, 9, 81, 6561, each the square of the one before (3 itself counts as the first), that is $3^{1}$ 3^1, $3^{2}$ 3^2, $3^{4}$ 3^4 and $3^{8}$ 3^8. The next square is $3^{16}$ 3^{16}, already past 13.

The first thing to try is to square while it fits and then finish one factor at a time. After 6561 we need 5 more factors of 3: 19683, 59049, 177147, 531441, 1594323. That is 3 + 5 = 8 multiplications instead of 12, and the answer is right.

A row of nine cells: 3, 9, 81 and 6561 in blue, labelled 3 to the powers 1, 2, 4 and 8, made by squaring; then 19683, 59049, 177147, 531441 and 1594323 in red, labelled 3 to the powers 9 to 13, made by multiplying by 3 one at a time.

Fig. 2. Squaring up to $3^{8}$ 3^8, then multiplying by 3 five times (red).

Look at the red tail in Fig. 2: it is the slow method again, only on a smaller exponent. For an exponent of $1 0^{18}$ 10^{18}, the squares get as far as an exponent of about $5.8 \cdot 1 0^{17}$ 5.8 \cdot 10^{17}, and about $4.2 \cdot 1 0^{17}$ 4.2 \cdot 10^{17} factors are left to multiply one at a time. At a rough rule of thumb of $1 0^{8}$ 10^8 simple steps per second, that is more than a hundred years.

So the question is how to get the missing $3^{5}$ 3^5 fast. It is a power too, so the same trick works on it. The largest square that fits into 5 factors is $3^{4}$ 3^4, and it leaves $3^{1}$ 3^1, and both are already in the row of squares. So $3^{13} = 3^{8} \cdot 3^{4} \cdot 3^{1}$ 3^{13} = 3^8 \cdot 3^4 \cdot 3^1: multiply 6561, 81 and 3, and get 531441 and then 1594323. That is 3 squarings and 2 more multiplications, 5 in total.

The squares we used are the ones whose exponents add up to 13: 8 + 4 + 1. That is 13 written in binary, 1101, where the digits 1 stand at the places worth 8, 4 and 1. Every whole number can be written in binary, so this works for any exponent.

We don't need to build the whole row first. Starting from the largest square needs the whole row in advance, but the last binary digit of a number is easy to read: it is the remainder after dividing by 2. So we go from the smallest square, and keep two numbers: the answer so far, starting at 1, and the current square, starting at 3. Each round takes the digit: if the last binary digit of the exponent left is 1, the current square is taken, that is, multiplied into the answer so far; if it is 0, the square is skipped. Then the round halves the exponent left, which drops that digit, and squares the current square. The last digit of 13 (1101) is the digit worth 1, and it is 1, so 3 is taken and the answer so far becomes 3. What is still missing is $3^{12}$ 3^{12}, which is $9^{6}$ 9^6: counted in factors of the new current square, 9, the exponent left is 6. That is 1101 with the last digit dropped, 110, and its last digit is the digit of 13 worth 2.

This is binary exponentiation, also called fast power or exponentiation by squaring.

Fig. 3 runs the whole computation of $3^{13}$ 3^{13}, one round per binary digit of 13.

Ten steps of computing 3 to the power 13: a row of squares 3, 9, 81, 6561 and 43046721 grows on top, with the taken squares 3, 81 and 6561 green and the skipped ones grey, and below it n goes 13, 6, 3, 1, 0 with its binary digits, the answer so far goes 1, 3, 243, 1594323, and the current square goes 3, 9, 81, 6561, 43046721.

Start: the exponent left is 13, binary 1101, the answer so far is 1 and the current square is 3.

1/10

Fig. 3. Step through binary exponentiation with the buttons or the arrow keys. Green squares are taken, grey ones are skipped.

The run made 4 rounds, one for each binary digit of 13, and each round had at most two multiplications: one to take the digit, one to square. An exponent of $1 0^{18}$ 10^{18} has 60 binary digits, so it needs at most 120 multiplications.

Remainders of huge powers

$3^{1 0^{18}}$ 3^{10^{18}} has about $4.8 \cdot 1 0^{17}$ 4.8 \cdot 10^{17} digits, and no computer can store it. So problems ask for the power modulo some number, that is, its remainder after dividing by that number: modulo $1 0^{9} + 7$ 10^9 + 7, or modulo 100 to get the last two digits. The remainder of a product depends only on the remainders of the factors (modular arithmetic). For example, 6561 · 81 = 531441 and 61 · 81 = 4941 end in the same two digits, 41. So we can take the remainder after every multiplication, and the numbers never grow past the square of the modulus.

Fig. 4 repeats the run on $3^{13}$ 3^{13}, keeping only the last two digits after every multiplication. It ends with 23, and 1594323 does end in 23.

The squares of 3 modulo 100, labelled 3 to the power 1, 2, 4 and 8: 3 taken, 9 skipped, 81 taken, and 6561 crossed out and replaced by its last two digits 61, taken; under each taken square the answer so far, 3, then 43, then the final answer 23.

Fig. 4. The run on $3^{13}$ 3^{13} modulo 100.

Implementation

In code, the exponent left is n, the current square is base, and the answer so far is result. The input is a, n and m, and the program prints a to the power n modulo m. All of them are long long, because n can be $1 0^{18}$ 10^{18} and a product of two remainders below $1 0^{9} + 7$ 10^9 + 7 reaches about $1 0^{18}$ 10^{18}.

Start: base = a % m, result = 1.
Take the digit: if n % 2 == 1, the last binary digit of n is 1, so result = result * base % m.
Square: base = base * base % m.
Halve: n /= 2 drops the last binary digit of n.
Stop when n == 0; result is the answer. Otherwise repeat from step 2.

n % 2 and n / 2 are the remainder and the integer division by 2. The loop needs n ≥ 0 and m ≥ 1. For a negative n it doesn't run at all and prints 1; negative powers modulo a prime need the modular inverse. For $3^{13}$ 3^{13} modulo $1 0^{9} + 7$ 10^9 + 7 the input is:

3 13 1000000007

#include <iostream>
using namespace std;

int main() {
    long long a, n, m;
    cin >> a >> n >> m;

    // start
    long long base = (a % m + m) % m;   // a can be negative, see Common mistakes
    long long result = 1 % m;           // 1 % m is 0 when m is 1

    while (n > 0) {                     // stop when n is 0
        // take the digit
        if (n % 2 == 1) result = result * base % m;
        // square
        base = base * base % m;
        // halve
        n /= 2;
    }

    cout << result << '\n';
}

All the numbers in the run on $3^{13}$ 3^{13} stay below $1 0^{9} + 7$ 10^9 + 7, so the remainders change nothing and the program prints 1594323. The loop runs once per binary digit of n, so it takes $O (lo g n)$ O(\log n) time.

Recursive version

The same idea can be written as a recursive function. If n is even, a to the power n is the square of a to the power n / 2. If n is odd, the integer division n / 2 loses one factor, so multiply by a once more. For $3^{13}$ 3^{13} the calls go down to the exponents 6, 3, 1 and 0, and on the way back the results are 1, 3, 27, 729 and 1594323.

long long power(long long a, long long n, long long m) {   // a already reduced: 0 <= a < m
    if (n == 0) return 1 % m;
    long long half = power(a, n / 2, m);
    long long result = half * half % m;
    if (n % 2 == 1) result = result * a % m;   // the factor lost by n / 2
    return result;
}

C++

It reads the binary digits of n from the left, so the values it builds on the way differ from the loop's values of result, but it also makes at most two multiplications per digit. The recursion is at most 61 calls deep for n up to $1 0^{18}$ 10^{18}.

This function forgets the odd case:

long long pw(long long a, long long n) {
    if (n == 0) return 1;
    long long h = pw(a, n / 2);
    return h * h;
}

What does pw(3, 5) return?

243

Without a remainder

When the answer itself fits in long long (problems that promise it is at most $1 0^{18}$ 10^{18}), drop % m. If it doesn't fit, result overflows without any error, so check the limits first. The loop's last squaring computes a square that is never used, and that square can overflow even when the answer fits: for a = $1 0^{9}$ 10^9 and n = 2 the answer $1 0^{18}$ 10^{18} fits, but the loop then squares $1 0^{18}$ 10^{18}. So square only while another digit is left, before halving:

if (n % 2 == 1) result *= base;
if (n > 1) base *= base;   // skip the last square, it is never used
n /= 2;

Correctness proof

Call the input $a$ and $n_{0}$ n_0, and let $r$ , $b$ and $n$ be the values of result, base and n at the start of a round. As in the walk through $3^{13}$ 3^{13}, the answer so far, times the current square raised to the exponent left, always equals the power we want: $r \cdot b^{n} = a^{n_{0}}$ r \cdot b^n = a^{n_0}. After the first round of $3^{13}$ 3^{13}, for example, 3 · $9^{6}$ 9^6 = $3^{13}$ 3^{13}. At the start $r = 1$ r = 1, $b = a$ b = a and $n = n_{0}$ n = n_0, so it holds.

Take a round. If $n$ is even, $b^{n} = (b^{2})^{n /2}$ b^n = (b^2)^{n/2}: the round squares $b$ and halves $n$ and leaves $r$ alone, so the product doesn't change. If $n$ is odd, $b^{n} = b \cdot (b^{2})^{(n - 1) /2}$ b^n = b \cdot (b^2)^{(n-1)/2}: the round moves one factor $b$ into $r$ , then squares $b$ , and the integer division gives $(n - 1) /2$ (n - 1) / 2. The product doesn't change either.

The loop stops at $n = 0$ n = 0, where $b^{0} = 1$ b^0 = 1, so $r = a^{n_{0}}$ r = a^{n_0}. With remainders, every step only replaces a number by its remainder modulo $m$ , and that keeps the remainder of the product. The loop stops because $n$ halves every round.

Five columns, one for the start of each round of computing 3 to the power 13 and one for the end: 1 · 3^13, 3 · 9^6, 3 · 81^3, 243 · 6561^1 and 1594323 · 43046721^0, each with "= 3^13" in green below it.

Fig. 5. $r \cdot b^{n}$ r \cdot b^n at the start of each round and at the end.

Common mistakes

Forgetting the odd exponent in the recursive version. With the base case n == 1 returning a, pw(3, 5) returns 81 and pw(3, 13) returns 6561 instead of 1594323. Multiply by a once more when n is odd.
Taking % m only at the end. On 3 40 1000000007 the long long product $3^{40}$ 3^{40}, about $1.2 \cdot 1 0^{19}$ 1.2 \cdot 10^{19}, overflows, and the program printed −629072825 instead of 953271190. Signed overflow is undefined, so another compiler may print something else. Take % m after every multiplication.
Not reducing a first. On 1000000000000000000 2 1000000007 the first base * base overflows, and the program printed −908814443 instead of 2401. Start with base = a % m.
Using pow from <cmath>. It computes in double, so (long long)pow(3, 39) gives 4052555153018976256 while the true value is 4052555153018976267, and pow(a, n) % m doesn't compile. Use the integer loop.
Forgetting that % in C++ keeps the sign. With base = a % m, the input -2 3 5 prints −3 instead of 2. Start with base = (a % m + m) % m.

Similar algorithms

Matrix exponentiation runs the same loop on square matrices instead of numbers. It finds the $n$ -th Fibonacci number, or any term of a linear recurrence, for $n$ up to $1 0^{18}$ 10^{18}.

Division modulo a prime $p$ is multiplication by the modular inverse, $b^{p - 2}$ b^{p-2} modulo $p$ for $b$ not divisible by $p$ , which is computed with binary exponentiation. When the modulus is not prime, the inverse comes from the extended Euclidean algorithm.

When a problem needs every digit of a big power, such as $9^{7000}$ 9^{7000}, no remainder helps: the number has to be stored digit by digit with long arithmetic.