Fast powers

How would you work out 3¹³? The obvious way is to start with 3 and multiply by 3 twelve more times. But there is a shortcut. 3¹³ is 3 × 3¹², and 3¹² is the square of 3⁶, and 3⁶ is the square of 3³. Each step halves the exponent. Watch how few multiplications that leaves.

1#include <stdio.h>2 3static int mults = 0;4 5long long power(long long x, int n)6{7    if (n == 0)8        return 1;9    long long half = power(x, n / 2);10    long long result = half * half;11    mults++;12    if (n % 2 == 1) {13        result = result * x;14        mults++;15    }16    return result;17}18 19long long power_loop(long long x, int n)20{21    long long result = x;22    for (int i = 1; i < n; i++) {23        result = result * x;24        mults++;25    }26    return result;27}28 29int main(void)30{31    long long a = power(3, 13);32    printf("power: %lld in %d multiplications\n", a, mults);33    mults = 0;34    long long b = power_loop(3, 13);35    printf("loop: %lld in %d multiplications\n", b, mults);36    return 0;37}

Halving the exponent

Here is the rule the program uses. When n is even, xⁿ is x^n/2 times itself. When n is odd, take one x out first:

The brackets ⌊ ⌋ mean “rounded down”, which is what n / 2 does in C when n is an int: 13 / 2 is 6. So it works for every n, not only powers of 2. In C:

long long power(long long x, int n)
{
    if (n == 0)
        return 1;
    long long half = power(x, n / 2);
    long long result = half * half;
    mults++;
    if (n % 2 == 1) {
        result = result * x;
        mults++;
    }
    return result;
}

In the clip, 13 halves to 6, then 3, then 1, then 0, and the base case returns 1. Then the cards return, the deepest first:

Seven multiplications, against the loop’s twelve.

The code and its cost

The bits of n

Look at the row of bits in the clip. Each call’s n is the one above it halved, so each call’s n % 2, odd or even, is one bit of 13 written in binary. Read from the deepest call to the first, the bits are 1, 1, 0, 1, and 1101 in binary is 8 + 4 + 0 + 1 = 13.

That gives the count. Every call with n above 0 makes one squaring, and there is one call for each bit of n. Every 1 bit adds one multiplication by x. So:

For 13, that is 4 + 3 = 7. The number of bits of n is ⌊log₂ n⌋ + 1, so the count is at most about 2 log₂ n. The loop always needs n − 1. In the notation of Time and space complexity (3.3), the halving version takes O(log n) time and the loop O(n). The stack holds one frame per bit, O(log n) as well.

The work here is the multiplications: everything else a call does is a test, a halving and a return. So fewer multiplications really is less work.

The loop that reads the bits

The same idea works without recursion. Go through the bits of n from the right. Keep x squared as you go, so at each bit x holds 3, 9, 81, 6561: the powers 3¹, 3², 3⁴, 3⁸. Whenever the bit is 1, multiply the result by the current x:

long long power_bits(long long x, int n)
{
    long long result = 1;
    while (n > 0) {
        if (n % 2 == 1)
            result = result * x;
        x = x * x;
        n = n / 2;
    }
    return result;
}

This uses one frame. It reads the same bits as the recursion, in the opposite order.

The plain loop

power_loop in the program is the obvious method: start with x and multiply by x, n − 1 times. It only works for n of 1 or more; for n = 0 it would return x instead of 1.

Running it yourself

gcc -std=c11 -Wall -Wextra -pedantic power.c -o power
./power

It prints power: 1594323 in 7 multiplications, then loop: 1594323 in 12 multiplications.

Worked example

The bit loop for 3¹³, by hand. Each row is the state at the top of the while loop:

At n = 13 the bit is 1, so result becomes 1 × 3 = 3. At n = 6 it is 0, and only x changes. At n = 3, result becomes 3 × 81 = 243. At n = 1, result becomes 243 × 6,561 = 1,594,323.

A million. In binary, 1,000,000 is 11110100001001000000: 20 bits, 7 of them 1. So power makes 20 + 7 = 27 multiplications, where the loop would make 999,999. A program that counts them, power-mod.c, agrees: 27.

Overflow. 3^1,000,000 has far too many digits for any C integer. The largest long long here is 9,223,372,036,854,775,807. 3³⁹ = 4,052,555,153,018,976,267 fits; 3⁴⁰ = 12,157,665,459,056,928,801 does not, and multiplying past it is undefined in C.

Remainders. To get past that, many programs keep only the remainder after dividing by some number m, at every step. The remainder of a product depends only on the remainders of the factors, so the answer comes out right. With m = 7, power-mod.c gets 3¹³ mod 7 = 3. Check: 1,594,323 = 7 × 227,760 + 3. It also gets 3^1,000,000 mod 1,000,000,007 = 64,935,414, in those 27 multiplications.

Where it’s used

Public-key cryptography. Encryption methods such as RSA raise a number to a power with hundreds of digits, keeping remainders at every step. The loop would never finish. Halving needs a few thousand multiplications at most.

Powers of a matrix. The same rule squares matrices instead of numbers. In linear algebra, Diagonalisation gives another way to raise a matrix to a power. A 2 × 2 matrix raised to the power n by halving gives fib(n) in O(log n) multiplications, against the exponential calls of Tree recursion (4.3).

Maths libraries. Functions that raise a number to a whole-number power use this method, because it needs so few multiplications.

Reference card

Platform: C11, gcc on x86-64 Linux: int 4 bytes, long long 8. The cards in the picture show names and values only.