Fast exponentiation

Compute 3 to the power 10¹⁸, or rather the remainder of this number divided by a billion and seven. Multiplying by three over and over would take a billion billion multiplications: at a billion multiplications a second, more than thirty years. The algorithm below does it in 84 multiplications. It looks at the exponent in binary, and each binary digit costs it one squaring and, if the digit is a one, one extra multiplication.

Change the exponent. There are as many steps as digits in the binary notation, and as many highlighted rungs as ones in it. A billion billion has 60 binary digits, 24 of them ones, hence 60 + 24 = 84.

The idea

a¹³ = a⁸ · a⁴ · a¹, because 13 = 1101₂. Powers of the form a, a², a⁴, a⁸, … are obtained from one another by a single squaring, and the ones in the binary notation of the exponent choose which of them are needed.

The same can be said recursively: aⁿ = (a^n/2)² for even n and aⁿ = a · aⁿ⁻¹ for odd n. Every two steps the exponent is at least halved: this is divide and conquer in one line. The iterative version is more convenient: go over the bits from right to left, keeping the current square and the current answer.

The main thing in problems “modulo”: take the remainder after every multiplication. The remainder of a product equals the remainder of the product of the remainders, so the answer does not change, but the numbers never grow beyond m². Without this, 3^10¹⁸ is a number of about 4.8·10¹⁷ decimal digits: no memory could hold it.

Code

In Python this is built in: pow(a, n, m) with three arguments does exactly this. In Go there is big.Int.Exp for big numbers. For ordinary uint64 the function is written by hand, and there is a trap: the product of two remainders modulo about 10⁹ fits in 64 bits, but modulo about 10¹⁸ it no longer does. Large moduli need a multiplication with a 128-bit result (bits.Mul64) or math/big.

The last line of the output is the inverse element. By Fermat’s little theorem, for a prime p the number a^p−1 leaves remainder 1, so a^p−2 is the inverse of a. This is a second way to divide modulo a number besides the extended Euclidean algorithm, but only for a prime modulus.

Why so few multiplications

Let the binary notation of n have k = ⌊log₂ n⌋ + 1 digits, s of them ones. The loop makes k iterations, each with one squaring, and the extra multiplication happens only on the ones: k + s ≤ 2k multiplications. For n about 10¹⁸ this is at most 120.

In general one can hardly do with fewer: one multiplication can at most double the exponent of a power, so getting aⁿ takes at least log₂ n multiplications. There are also clever addition chains that save some of the extra multiplications, and they are sought for fixed large exponents in cryptography, but the order log n remains.

A key agreed in plain sight

The most important application of fast exponentiation works every time you open a secure connection. Two parties need to agree on a shared secret key, but any exchange between them may be overheard. In 1976 Whitfield Diffie and Martin Hellman proposed a way:

1. Everyone knows a large prime p and a number g.
2. Alice picks a secret a and openly sends A = gᵃ mod p. Bob picks a secret b and sends B = gᵇ mod p.
3. Alice computes Bᵃ, Bob computes Aᵇ. Both get g^ab mod p, the shared key.

Below is a real prime from the standard RFC 3526: 2048 bits, more than six hundred decimal digits. SSH still uses it. Secrets of the same length are chosen by your browser’s random number generator.

Look at the time: four powers of a 2048-bit number take a fraction of a second right in the browser, about three thousand multiplications each. And Eve sees p, g, A and B, and to find the key Eve has to recover a from A. This is the discrete logarithm problem, and for numbers like these nobody knows how to solve it in reasonable time. The whole scheme rests on an asymmetry: raising to a power modulo a number is easy, taking the logarithm is hard.

RSA uses the same operation: encryption is mᵉ mod N, decryption is c^d mod N. The second problem is about this.

What it costs

O(log n) multiplications. If the numbers fit in a machine word, each multiplication is O(1), and the whole power is a few dozen clock cycles. For long numbers of L bits schoolbook multiplication costs O(L²), so raising to an L-bit exponent modulo an L-bit number is O(L³). For 2048 bits this is billions of elementary operations, and still a fraction of a second.

Not only numbers

The algorithm does not care what it multiplies, as long as the multiplication is associative. Matrices, by all means: raising the matrix [[1, 1], [1, 0]] to the power n gives Fibonacci numbers, and the n-th Fibonacci number modulo a number is computed in O(log n); that is the first problem, and in detail matrices and recurrences. Permutations: in log n steps you can find out where a card ends up after a billion identical shuffles. Polynomials, transformations, paths of length n in a graph: wherever there is “do the same thing n times”, the same ladder of squares fits.

History

The technique is ancient. Donald Knuth in “The Art of Computer Programming” points out that it appears as early as Pingala’s treatise on prosody, written before 200 BC: there 2ⁿ, the number of possible rhythms of n syllables, is computed this way. It came into computing with the first computers and became an unnoticed foundation of public-key cryptography.

Problems

In the first problem it is a matrix that has to be raised to a power. In the second, a whole toy RSA has to be assembled: Euclid, fast exponentiation and one trap with overflow.