Binary Exponentiation: Compute a^b mod m in O(log b)

Many problems ask for an answer “modulo 10⁹ + 7”, and many of those need a large power such as 2ⁿ or a^b with b up to 10¹⁸. Multiplying a by itself b times is hopeless. Binary exponentiation does it in about 60 multiplications.

The idea: square instead of multiply

Notice that a⁸ = ((a²)²)². That is only three multiplications instead of seven. Repeated squaring gives us a, a², a⁴, a⁸, a¹⁶, … very cheaply.

Any exponent can be written as a sum of powers of two, which is exactly its binary representation. For example 13 = 8 + 4 + 1, which is 1101 in binary. So

a¹³ = a⁸ × a⁴ × a¹

We walk through the bits of b from the lowest. At each step, the current a holds a^2k. If the bit is 1, multiply it into the result. Then square a and shift b right.

Tracing a = 3, b = 13

b (binary)	lowest bit	current a	result
1101	1	3	1 × 3 = 3
110	0	9	3
11	1	81	3 × 81 = 243
1	1	6561	243 × 6561 = 1 594 323

And 3¹³ is indeed 1 594 323, found with a handful of multiplications.

Code (C++)

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 1'000'000'007;

// Computes (a^b) mod m in O(log b) multiplications.
long long power(long long a, long long b, long long m) {
    long long result = 1 % m;
    a %= m;
    while (b > 0) {
        if (b & 1) result = result * a % m;   // this bit of b is set
        a = a * a % m;                        // a, a^2, a^4, a^8, ...
        b >>= 1;
    }
    return result;
}

int main() {
    int q;
    cin >> q;
    while (q--) {
        long long a, b;
        cin >> a >> b;
        cout << power(a, b, MOD) << '\n';
    }
    return 0;
}

The input format matches CSES 1095 – Exponentiation: a number of queries, then a and b on each line.

Complexity

The loop runs once per bit of b, so O(log b) multiplications. For b up to 10¹⁸ that is at most 60 iterations.

Details that cause Wrong Answer

• Overflow. Both factors are below 10⁹ + 7, so their product is below about 10¹⁸, which fits in long long. With int it overflows. If the modulus is close to 10¹⁸, even long long is not enough and you need __int128.
• Reduce a first. If a itself is larger than the modulus, take a %= m before the loop.
• b = 0. Anything to the power 0 is 1, including 0⁰ by the usual convention in these problems. The code starts result at 1 and handles it without a special case.

Where it shows up

• Modular inverse. When the modulus p is prime, Fermat's little theorem gives a⁻¹ ≡ a^p−2 (mod p). That is how you “divide” under a modulus, for example when computing binomial coefficients.
• Huge exponents. To compute a^bc mod p, reduce the exponent b^c modulo p − 1 first (again by Fermat), then use binary exponentiation twice.
• Matrix exponentiation. The same loop works with matrices instead of numbers, which computes the n-th Fibonacci number in O(log n).

Practice problems

• CSES 1095 – Exponentiation: the problem above.
• CSES 1712 – Exponentiation II: the a^bc version.