Fast Exponentiation

autor: iamtwz, billchenchina, CBW2007, CCXXXI, chinggg, Enter-tainer, eyedeng, FFjet, gaojude, Great-designer, H-J-Granger, Henry-ZHR, hsfzLZH1, Ir1d, kenlig, Konano, ksyx, luoguyuntianming, Marcythm, Menci, NachtgeistW, ouuan, Peanut-Tang, qwqAutomaton, sshwy, StudyingFather, Tiphereth-A, TrisolarisHD, TRSWNCA, Xeonacid, Yuuko10032, Zhangjiacheng2006, Zhoier, Hszzzx, shenshuaijie, kfy666

Introduction

Fast exponentiation (fast exponentiation), also called binary exponentiation (binary exponentiation) or exponentiation by squaring (exponentiation by squaring), is a small technique for computing $a^n$ in $\Theta(\log n)$ time, whereas brute-force computation requires $\Theta(n)$ time.

This technique can be applied whenever multiplication involving $a$ is associative, such as modular exponentiation and matrix powers. See the later Applications section for details.

Process

Computing the $n$-th power of $a$ means multiplying $n$ copies of $a$ together: $a^{n} = \underbrace{a \times a \cdots \times a}_{n\text{ 个 a}}$. However, this method is not suitable when $n$ is too large or a single multiplication is too expensive. The idea of binary exponentiation is to divide the exponentiation task into smaller tasks according to the exponent's binary representation.

This is the basic idea of fast exponentiation. There are two common implementations.

Iterative version

Let the binary representation of $n$ be $(n_tn_{t-1}\cdots n_1n_0)_2$, that is,

where $n_i\in\{0,1\}$. Then,

Note that only terms with $n_i=1$ actually appear in the product.

According to this expression, we can first compute the $\Theta(\log n)$ values of the $2^k$-th powers of $a$ in $\Theta(\log n)$ time, and then spend $\Theta(\log n)$ time multiplying into the final result the powers corresponding to binary bits equal to $1$. This is the iterative implementation of fast exponentiation.

The pseudocode is as follows:

Computing fast exponentiation this way requires $\Theta(\log n)$ multiplications.

Recursive version

The same process can also be implemented recursively. Notice that the binary expansion of the exponent $n$ can be written recursively as

Therefore, $a^n$ can be computed recursively as

This is the recursive implementation of fast exponentiation.

The pseudocode is as follows:

Computing fast exponentiation this way requires $\Theta(\log n)$ recursive calls and also $\Theta(\log n)$ multiplications. Although the complexities are the same, the iterative version is faster in practice because recursion has some overhead.

Applications

Modular exponentiation

This is a very common application; for example, it can be used to compute a modular multiplicative inverse. Since taking a modulus does not interfere with multiplication, we only need to take the modulus during the computation.

First, we can implement it directly according to the recursive method above:

long long binpow(long long a, long long b, long long p) {
  if (b == 0) return 1;
  long long res = binpow(a, b / 2, p);
  if (b % 2)
    return res * res % p * a % p;
  else
    return res * res % p;
}

def binpow(a, b, p):
    if b == 0:
        return 1
    res = binpow(a, b // 2, p)
    if (b % 2) == 1:
        return res * res * a % p
    else:
        return res * res % p

The second implementation is non-recursive. It multiplies into the answer the powers corresponding to binary bits equal to 1 during the loop. Although the theoretical complexities of the two are the same, the second is faster in practice because recursion incurs some overhead.

long long binpow(long long a, long long b, long long p) {
  long long res = 1;
  while (b > 0) {
    if (b & 1) res = res * a % p;
    a = a * a % p;
    b >>= 1;
  }
  return res;
}

def binpow(a, b, p):
    res = 1
    while b > 0:
        if b & 1:
            res = res * a % p
        a = a * a % p
        b >>= 1
    return res

Computing Fibonacci numbers

According to the recurrence $F_n = F_{n-1} + F_{n-2}$ of the Fibonacci sequence, we can construct a $2\times 2$ matrix to represent the transformation from $F_i,F_{i+1}$ to $F_{i+1},F_{i+2}$. Then, when computing the $n$-th power of this matrix, we can use the idea of fast exponentiation to obtain the result in $\Theta(\log n)$ time. For more details, see Fibonacci sequence; an implementation of matrix fast exponentiation is available in Matrix-accelerated recurrence.

Multiple permutations

Simply raise the permutation to the $k$-th power and apply it to the sequence. The time complexity is $O(n \log k)$. For more details, see Composition of permutations.

Accelerating operations on point sets in geometry

Referring to Vectors and matrices, each operation can be represented by a transformation matrix, and a sequence of consecutive transformations can be represented by a product of matrices. A Repeat operation is equivalent to raising a matrix to the $k$-th power. Thus, the final matrix of the entire transformation sequence can be computed in $O(m \log k)$ time. Applying it to the $n$ points gives total complexity $O(n + m \log k)$.

Counting paths of fixed length

Raise the adjacency matrix $M$ of the graph to the $k$-th power. Then $M_{i,j}$ is the number of paths of length $k$ from $i$ to $j$. The complexity is $O(n^3 \log k)$. See the Matrix page for details.

Modular integer multiplication

As with binary exponentiation, this time we express one of the factors as a sum of integer powers of 2. When multiplying a number by 2 and taking a modulus, we can use addition and subtraction to avoid integer overflow. This solves the problem in $O (\log m)$ time. The recursive method is:

In practice, however, this method is inefficient because it introduces greater computational complexity. In actual programming, fast multiplication is usually used for multiplication when the modulus range is long long.

High-precision fast exponentiation

Prerequisite: Big-integer multiplication

#include <cmath>
#include <cstring>
#include <iostream>
using namespace std;

const int M = 500;

int a[505], b[505], t[505];

// big-integer multiplication
void mult(int x[], int y[]) {
  memset(t, 0, sizeof(t));
  for (int i = 1; i <= x[0]; i++) {
    for (int j = 1; j <= y[0]; j++) {
      if (i + j - 1 > M) continue;
      t[i + j - 1] += x[i] * y[j];
      t[i + j] += t[i + j - 1] / 10;
      t[i + j - 1] %= 10;
      t[0] = i + j;
    }
  }
  memcpy(b, t, sizeof(b));
}

// fast exponentiation
void binpow(int p) {
  if (p == 1) {
    memcpy(b, a, sizeof(b));
    return;
  }
  binpow(p / 2);  // (2^(p/2))^2=2^p
  mult(b, b);     // square b
  if (p % 2 == 1) mult(b, a);
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int p;
  cin >> p;
  a[0] = 1;  // record the number of digits in array a
  a[1] = 2;  // square 2
  b[0] = 1;  // record the number of digits in array b
  b[1] = 1;  // answer array
  binpow(p);
  cout << (int)(log10(2) * p) + 1 << '\n';
  b[1] -= 1;  // subtract 1 from the last digit
  for (int i = M; i >= 1; i--) {
    cout << b[i];
    if ((i - 1) % 50 == 0) {
      cout << '\n';
    }
  }
}

Precomputed fast exponentiation with a fixed base

When the base $a$ is fixed, we can use block decomposition to preprocess for some amount of time and answer each power query in $O(1)$ time. This algorithm is also called light-speed exponentiation. The process is as follows:

1. Select a number $s$. Precompute the values of $a^0,a^1,\cdots,a^{s-1}$ and $a^0,a^s,\cdots,a^{\lfloor p/s\rfloor s}$ and store them in two arrays.
2. For each query $a^b$, decompose $b$ as $\lfloor b/s\rfloor s+(b\bmod s)$. Then $a^b=a^{\lfloor b/s\rfloor s}\cdot a^{b\bmod s}$, so the answer can be obtained in $O(1)$ time.

Suppose the range of the exponent $b$ is $[0,n]$. The block length $s$ is often chosen as $\sqrt{n}$ or a nearby power of $2$. Choosing $\sqrt{n}$ gives the optimal preprocessing complexity $O(\sqrt{n})$, while choosing a power of $2$ allows bit operations to simplify the computation.

In particular, for modular exponentiation, having the same base $a$ also implies that the modulus $m$ must be the same. By Euler's theorem extension, for any modulus $m$, the upper bound of the preprocessed exponent range is $n = 2\varphi(m)$; for a prime modulus $p$, the upper bound is $n = p - 1$. The preprocessing complexity in both cases is $O(\sqrt{m})$.

int a, mod, pow1[65536], pow2[65536];

void preproc() {
  pow1[0] = pow2[0] = 1;
  for (int i = 1; i < 65536; i++) pow1[i] = 1LL * pow1[i - 1] * a % mod;
  int pow65536 = 1LL * pow1[65535] * a % mod;
  for (int i = 1; i < 65536; i++) pow2[i] = 1LL * pow2[i - 1] * pow65536 % mod;
}

int query(int pows) {
  return 1LL * pow1[pows & 65535] * pow2[pows >> 16] % mod;
}

Exercises

• UVa 1230 - MODEX
• UVa 374 - Big Mod
• UVa 11029 - Leading and Trailing
• Codeforces - Parking Lot
• SPOJ - The last digit
• SPOJ - Locker
• SPOJ - Just add it

Some content on this page is translated from the blog post Binary Exponentiation and its English translation Binary Exponentiation. The Russian version is licensed under Public Domain + Leave a Link; the English version is licensed under CC-BY-SA 4.0.