Skip to content

Fast Exponentiation

Introduction

Fast exponentiation, also known as binary exponentiation or exponentiation by squaring, is a technique for computing \(a^n\) in \(\Theta(\log n)\) time, whereas the naive approach requires \(\Theta(n)\) time.

This technique can be applied in any scenario where multiplication of \(a\) is associative, such as modular exponentiation, matrix exponentiation, etc. For more details, see the Applications section.

Process

Computing \(a\) raised to the \(n\)-th power means multiplying \(n\) copies of \(a\) together: \(a^{n} = \underbrace{a \times a \cdots \times a}_{n\text{ copies of a}}\). However, when \(n\) is too large or each multiplication is too expensive, this method becomes impractical. The idea of binary exponentiation is to break down the exponentiation task into smaller tasks based on the binary representation of the exponent.

Example

Suppose we want to compute \(3^{13}\). Expanding it as a product requires \(13-1=12\) multiplications. However, since

\[ 3^{13} = 3^{(1101)_2} = 3^8 \times 3^4 \times 3^1, \]

we can compute \(3^{13}\) with only \(2\) multiplications if we can quickly compute \(3^{1},3^{2},3^{4},3^{8}\). Thus, we only need a fast way to compute the sequence of \(3\) raised to powers of \(2^k\). This is easy because each element in the sequence (except the first) is the square of the previous element.

Based on this analysis, the computation of \(3^{13}\) proceeds as follows:

\[ \begin{aligned} 3^1 &= 3, \\ 3^2 &= \left(3^1\right)^2 = 3^2 = 9, \\ 3^4 &= \left(3^2\right)^2 = 9^2 = 81, \\ 3^8 &= \left(3^4\right)^2 = 81^2 = 6561, \\ 3^{13} &= 6561 \times 81 \times 3 = 1594323. \end{aligned} \]

Only \(5\) multiplication operations were performed.

This is the basic idea of fast exponentiation. As for specific implementations, there are two common versions.

Iterative Version

Let the binary representation of \(n\) be \((n_tn_{t-1}\cdots n_1n_0)_2\), i.e.,

\[ n = n_t2^t + n_{t-1}2^{t-1} + \cdots + n_12^1 + n_02^0, \]

where \(n_i\in\{0,1\}\). Then we have

\[ \begin{aligned} a^n & = a^{n_t2^t + n_{t-1}2^{t-1} + \cdots + n_12^1 + n_02^0}\\ & = a^{n_0 2^0} \times a^{n_1 2^1}\times \cdots \times a^{n_{t-1}2^{t-1}} \times a^{n_t2^t}. \end{aligned} \]

Note that only terms where \(n_i=1\) actually contribute to the product.

Based on this expression, we can first compute the \(\Theta(\log n)\) powers \(a^{2^k}\) in \(\Theta(\log n)\) time, then spend \(\Theta(\log n)\) time to multiply those corresponding to binary digits equal to \(1\) into the final result. This is the iterative version of fast exponentiation.

Pseudocode:

\[ \begin{array}{l} \textbf{Algorithm }\text{FastPow}(a, n): \\ \textbf{Input. }\text{Base }a\text{ and exponent }n.\\ \textbf{Output. }\text{Power }a^n.\\ \textbf{Method.}\\ \begin{array}{ll} 1 & \textit{result}\gets\mathrm{Id}\\ 2 & \textbf{while }n > 0\textbf{ do}\\ 3 & \qquad \textbf{if }n \bmod 2 = 1\textbf{ then}\\ 4 & \qquad \qquad \textit{result} \gets \textit{result}\cdot a\\ 5 & \qquad \textbf{end if}\\ 6 & \qquad a \gets a \cdot a\\ 7 & \qquad n \gets n / 2\\ 8 & \textbf{end while}\\ 9 & \textbf{return }\textit{result} \end{array} \end{array} \]

Using this method to compute fast exponentiation requires \(\Theta(\log n)\) multiplication operations.

Recursive Version

This process can also be implemented recursively. Note that the binary expansion of exponent \(n\) can be written recursively as

\[ (n_tn_{t-1}\cdots n_1n_0)_2 = 2 \times (n_tn_{t-1}\cdots n_1)_2 + n_0. \]

Therefore, the power \(a^n\) can be computed recursively as

\[ a^n = \begin{cases} 1, & n = 0,\\ (a^{\lfloor n/2\rfloor})^2, & n > 0 \text{ and }n\text{ is even},\\ (a^{\lfloor n/2\rfloor})^2\cdot a, & n > 0 \text{ and }n\text{ is odd}.\\ \end{cases} \]

This is the recursive version of fast exponentiation.

Pseudocode:

\[ \begin{array}{l} \textbf{Algorithm }\text{FastPow}(a, n): \\ \textbf{Input. }\text{Base }a\text{ and exponent }n.\\ \textbf{Output. }\text{Power }a^n.\\ \textbf{Method.}\\ \begin{array}{ll} 1 & \textbf{if }n = 0\textbf{ then}\\ 2 & \qquad \textbf{return }\mathrm{Id}\\ 3 & \textbf{end if}\\ 4 & \textit{result} \gets \text{FastPow}(a, n / 2) \\ 5 & \textbf{if }n\bmod 2 = 0\textbf{ then}\\ 6 & \qquad \textbf{return }\textit{result}\cdot\textit{result}\\ 7 & \textbf{else}\\ 8 & \qquad \textbf{return }\textit{result}\cdot\textit{result}\cdot a\\ 9 & \textbf{end if} \end{array} \end{array} \]

Using this method to compute fast exponentiation requires \(\Theta(\log n)\) recursive calls and also \(\Theta(\log n)\) multiplication operations. Although the complexity is the same, the iterative version is faster in practice due to the overhead of recursion.

Applications

Modular Exponentiation

Luogu P1226【Template】Fast Exponentiation

Given three integers \(a,b,p\), compute \(a^b\bmod p\), where \(p\ge 2\).

This is a very common application, for example, it can be used to compute modular multiplicative inverses. Since we know that the modulo operation does not interfere with multiplication, we only need to take modulo during the computation.

First, we can directly implement the recursive method described above:

Reference Implementation 
1
2
3
4
5
6
7
8
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;
}

1
2
3
4
5
6
7
8
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. During the loop, it accumulates the powers corresponding to binary digits that are 1 into the answer. Although the theoretical complexity of both is the same, the second one is faster in practice because recursion incurs some overhead.

Reference Implementation 
1
2
3
4
5
6
7
8
9
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;
}

1
2
3
4
5
6
7
8
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
Note
  • The modulus is usually greater than \(1\). In very special cases, the modulus \(p\) may equal \(1\), in which case the situation \(b=0\) needs special consideration.
  • When the exponent is very large, it is necessary to reduce the exponent using the Extended Euler's Theorem before computation.

Computing Fibonacci Numbers

According to the recurrence relation of the Fibonacci sequence \(F_n = F_{n-1} + F_{n-2}\), 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 use the idea of fast exponentiation to compute the result in \(\Theta(\log n)\) time. For more details, see Fibonacci Numbers. The implementation of matrix fast exponentiation can be found in Matrix Acceleration for Recurrences.

Multiple Permutations

Problem Description

Given a sequence of length \(n\) and a permutation, apply this permutation to the sequence \(k\) times.

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.

Note

Construct a graph from this permutation, then raise each cycle to the \(k\)-th power (which is equivalent to \(k\) modulo the cycle length) separately. This can solve the problem in \(O(n)\) time.

Accelerating Operations on Point Sets in Geometry

HDU 4087 A Letter to Programmers

Given \(n\) points \(p_i\) in three-dimensional space, apply \(m\) operations to these points. The operations include:

  1. Shift: Move points along a certain vector.
  2. Scale: Scale the coordinates of points by a certain factor.
  3. Rotate: Rotate points around a certain line.

There is also a special operation: repeat a sequence of operations \(k\) times (Repeat), and Repeat operations can be nested. Output the coordinates of each point after all operations.

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

Fixed-Length Path Counting

Problem Description

Given a directed graph (edge weights are 1), count the number of paths of length \(k\) from any vertex \(u\) to any vertex \(v\).

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

Integer Multiplication Under Modulo

Problem Description

Given non-negative integers \(a,b\) and a positive integer \(m\), compute \(a\times b\bmod m\), where \(a,b\le m\le 10^{18}\).

Similar to the idea of binary exponentiation, this time we express one of the multiplicands as a sum of powers of 2. Because when multiplying a number by 2 and taking modulo, we can convert it to addition and subtraction to prevent integer overflow. This solves the problem in \(O (\log m)\) time. The recursive method is as follows:

\[ a \cdot b = \begin{cases} 0 &\text{if }a = 0 \\ 2 \cdot \frac{a}{2} \cdot b &\text{if }a > 0 \text{ and }a \text{ even} \\ 2 \cdot \frac{a-1}{2} \cdot b + b &\text{if }a > 0 \text{ and }a \text{ odd} \end{cases} \]

However, in practice, this method introduces greater computational complexity and is not time-efficient. In actual programming, fast multiplication is often used for multiplication operations when the modulus range is within long long.

High-Precision Fast Exponentiation

Prerequisite skill: Big Integer Multiplication

Luogu P1045 [NOIP 2003 Popularization Group] Mersenne Number

Given an integer \(P\) (\(1000 < P < 3100000\)), compute the number of digits and the last \(500\) digits (in decimal representation) of \(2^P−1\), padding with leading zeros if fewer than \(500\) digits.

Code Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include <cmath>
#include <cstring>
#include <iostream>
using namespace std;

const int M = 500;

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

// Large 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 of 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 Fixed Base

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

  1. Choose a number \(s\), precompute \(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\), split \(b\) into \(\lfloor b/s\rfloor s+(b\bmod s)\), then \(a^b=a^{\lfloor b/s\rfloor s}\cdot a^{b\bmod s}\), which can be answered in \(O(1)\).

Assuming the exponent \(b\) ranges in \([0,n]\), the block size \(s\) is often chosen as \(\sqrt{n}\) or a nearby power of \(2\). Choosing \(\sqrt{n}\) yields the optimal preprocessing complexity \(O(\sqrt{n})\), while choosing a power of \(2\) simplifies computation using bit operations.

In particular, for modular exponentiation, the same base \(a\) implies the same modulus \(m\). Due to the Extended Euler's Theorem, for any modulus \(m\), the upper bound for the exponent range in preprocessing is \(n = 2\varphi(m)\); for prime modulus \(p\), the upper bound is \(n = p - 1\). In both cases, the preprocessing complexity is \(O(\sqrt{m})\).

Reference Code 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
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

Parts of this page are translated from the blog post Бинарное возведение в степень 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.

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用