本文共 7816 字，大约阅读时间需要 26 分钟。

---

1. 快速幂概念

快速幂 Exponentiation By Squaring (平方求幂)是一种简单有效的算法，可以以 $O(\log n)$ 的时间复杂度计算幂次。它的应用场景比较常见，题型众多，而且有许多算法都以它为基础。

幂运算 $a^n$ 即 $n$ 个 $a$ 相乘，快速幂就是快速计算 $a^n$ 。如 $n=10^9$ 时，计算 $a^n$ 这样大的数即使用Java也不能处理：数字太大，计算时间太长。我们这里只考虑如何缩短计算时间，以 $9^9$ 次方为例。

最朴素的做法：逐个做乘法，$9\times 9=81,81\times 9=\dots$ ，这样一共进行 $8$ 次乘法。时间复杂度是 $O(n)$ ，即使能够算出来，也会超时。

改进一下：先算 $9^4=9\times 9\times 9\times 9=6561$ ，用到 $3$ 次乘法，再计算它的平方，$(9^4{)}^2=43046721$ ，最后再乘以 $9$ 。总共用到 $5$ 次乘法。但是这也不是最优解。

最优解：先计算 $9^2=81$ ，再计算 $(9^2{)}^2=9^4$ ，再计算 $(9^4{)}^2=9^8$ ，最后计算 $9^8\ast 9$ 。只用到 $4$ 次乘法。

这就是分治法的思路，复杂度为 $O(\log n)$ 。

(1) 快速幂递归实现

我们可以得到下面的递归方程：

$$a^n=\{\begin{array}{ll}{a}^{\frac{n-1}{2}}\cdot a^{\frac{n-1}{2}}\cdot a,& \text{if }n\text{ is odd}\\ a^{\frac{n}{2}}\cdot a^{\frac{n}{2}},& \text{if }n\text{ is even but not 0}\\ 1,& \text{if }n=0\end{array}$$ 计算 $a^n$ ，如果 $n=0$ ，递归基准，返回 $a^0=1$ ；如果 $n$ 是不为 $0$ 的偶数，则先计算 $a^{\frac{n}{2}}$ 再平方；如果 $n$ 是奇数，则先计算 $a^{n-1}$ 再乘以 $a$ 。整个过程十分自然，代码也很简单，基本就是递归方程的翻译：

int fastPow(int a, int n) {
   	if (n == 0) return 1; //0^0也定义为1	int temp = fastPow(a, n / 2);	if (n & 1) //奇数		return temp * temp * a;	else //偶数		return temp * temp;}

注意：$temp$ 变量的存在是必须的，如果不记录下来，就会计算两次 $a^{\frac{n}{2}}$ ，整个算法退化为 $O(n)$ 。另外，这个程序中的递归层数只有 ${\log }_{2}n$ ，不会溢出。

(2) 快速幂迭代实现

递归虽然简洁，但是存在额外的开销，把递归改写为循环，就可以避免栈空间的大量占有和函数调用带来的时间成本。为此，我们使用位运算进行快速幂，时间复杂度也是 $O(\log n)$ 。下面以 $a^{11}$ 为例说明算法的原理。

我们可以把 $a^{11}$ 分解成 $a^8,a^2,a^1$ 的乘积，即 $a^{11}=a^{8+2+1}=a^8\times a^2\times a^1$。如何求 $a^8,a^2,a^1$ 的值？似乎不需要分开计算，我们可以发现：$a^1\times a^1=a^2,a^2\times a^2=a^4,a^4\times a^4=a^8$ 等等，$1,2,4,8$ 则是 $2$ 的幂数，产生的 $a^i$ 都是倍乘关系，因此逐渐递推就可以了。在程序中，这会用 base *= base 实现。

那么如何把 $n$ 分解成 $11=8+2+1$ 这样的可以倍乘递推的关系呢？根据数学常识，每一个正整数可以唯一表示为若干个指数不重复的 $2$ 的次幂的和。即 $n$ 在二进制表示下有 $k$ 位，其中第 $i(0\le i<k)$ 位的数字是 $c_i(c_i=0,1)$ ，那么把 $n$ 转为二进制形式：$$n=c_{k-1}{c}_{k-2}\dots c_0$$

即为：$$n=c_{k-1}{2}^{k-1}+c_{k-2}{2}^{k-2}+\cdots +c_0{2}^0$$

于是：$$a^n=a^{c_{k-1}\ast 2^{k-1}}\ast a^{c_{k-2}\ast 2^{k-2}}\ast \cdots \ast a^{c_0\ast 2^0}$$

因为 $k=\lceil {\log }_2(n+1)\rceil$ ，所以上式中的乘积项的数量不多于 $\lceil {\log }_2(n+1)\rceil$ 。同时，二进制数的每一位的权值都是低一位的两倍，对应的 $a^{2^i}$ 则是倍乘的关系：$$a^{2^i}=(a^{2^{i-1}}{)}^2$$ 所以很容易通过 $k$ 次递推求出每个乘积项。如果还不理解，举一个二进制表示的例子。例如 $n=11_{10}=1011_2=2^3+2^1+2^0=8+2+1$ ，所以只需要把 $n$ 按照二进制处理即可。遍历 $n$ 在二进制表示下的所有数位 $c_i$ ，整个算法的时间复杂度为 $O(\log n)$ 。

还有一个问题，如何跳过那些不需要的 $a^{2^i}$ ？例如求 $a^{11}$ ，因为 $11=8+2+1$ ，需要跳过 $a^{2^2}=a^4$ 。这很简单，$1011_2$ 中的 $0$ 就是要跳过的，用位运算 $n\&1$ 做个判断即可：

• $n\&1$ ，取 $n$ 的最后一位，并且判断这一位是否要跳过，$c_i=0$ 时跳过， $c_i=1$ 时把该乘积项累积到答案中；
• $n>>=1$ ，将 $n$ 右移一位，把处理过的 $n$ 的最右一位丢弃。

代码如下，先看代码然后举例说明：

int fastPow(int a, int n) {
   	int base = a;	//不定义base,直接用a也可以	int res = 1;	//返回的结果	while (n) {
   		if (n & 1)  //如果n的最后一位是1,表示这个地方需要乘以a^i			res *= base;		base *= base; //推算乘积,a^2->a^4->a^8->a^16...		n >>= 1;    //n右移一位,把刚处理过的n的最后一位丢弃	}	return res;}

对照这一程序，$n=11$ 时执行过程如下：

	$n$	res(res *= base)	base(base *= base)
第一轮	$1011_2$	$a^1$	$a^2$
第二轮	$101_2$	$a^1\times a^2$	$a^4$
第三轮	$10_2$	是 $0$ ，res 不变	$a^8$
第四轮	$1_2$	$a^1\times a^2\times a^8$	$a^{16}$
结束	$0$

如果当时的二进制位为 $1$ ，就说明结果需要乘以 base 。res 也可能不变，而 base 总是在变化，一直进行倍乘。

---

2. 快速幂取模

由于幂运算的结果非常大，常常会超过变量类型的最大值，因此涉及到快速幂的题目，往往需要对大素数做取模操作，缩小结果。根据模运算的性质，在快速幂中做取模操作，对 $a^n$ 取模，和先对 $a$ 取模再做幂运算的结果是一样的。即：

$$a^n\mathrm{m}\mathrm{o}\mathrm{d}m=(a\mathrm{m}\mathrm{o}\mathrm{d}m{)}^n\mathrm{m}\mathrm{o}\mathrm{d}m$$

原则是步步取模。据此修改快速幂函数，加上取模操作。以HDU 2817为例，取模操作如下：

//如果mod较大,应该全部开long longconst int mod = 200907;int fastPow(int a, int n) {
   	if (n == 0) return 1; 	if (n & 1) //奇数		return fastPow(a, n - 1) * a % mod;	else {
   	   //偶数		int temp = fastPow(a, n / 2) % mod;		return temp * temp % mod;	}}//------------------------------------------int fastPow(int a, int n) {
    	int base = a;	int res = 1;	while (n) {
   		if (n & 1) //修改: res = (long long)res * base % mod;			res = (res * base) % mod; //步步取模		base = (base * base) % mod;   //步步取模		n >>= 1;	}	return res;}

上面给出的代码可能会出现错误，因为C++中两个数值进行算术运算时，以参与运算的最高数值类型为基准，与保存结果的变量类型无关。即两个32位整数的乘积可能超过 int 类型的表示范围，但是CPU只会用一个32位寄存器保存结果，造成越界现象，导致类似 res * base, base * base 这样的乘积结果出错，然后取模的结果也会出错。为此可以将其中一个数强制转换为64位整数类型 long long 参与运算，从而得到正确的乘积，对 mod 取模后执行赋值操作，隐式转换为 int 存储到 res 中。

更直接的方法是所有数值都使用 long long 类型，不过C++内置的最高整数类型只有64位……如果运算中的 an % mod 中的三个变量 a, n, mod 都在 1018 级别，则不会有一个可供使用的128位整数类型，需要使用一些特殊的处理方法。

---

3. 快速幂拓展

上述都是整数的快速幂。但是在算 $a^n$ 时，只要 $a$ 的数据类型支持乘法和满足结合律，快速幂的算法都是有效的。比如浮点数：

double fastPow(double base, int exponent) {
          if (exponent < 0) {
    //exponent为负数     	base = 1.0 / base;     	exponent = -exponent;     }     double res = 1; //结果      while (exponent) {
            if (exponent & 1)          	res *= base;         base *= base;         exponent >>= 1;     }     return res;}

一个模板如下，矩阵、高精度整数，都可以使用这个思路：

template 
   
    T fastPow(T a, long long n) {
    	T res = 1; //赋值为乘法单位元,可能需要根据T的构造函数进行修改	while (n) [		if (n & 1) 			ans = ans * a; //重载*即可; 用*=需要重载*=		a = a * a;		n >>= 1;	}	return res;}
   

注意，复杂类型的快速幂，时间复杂度需要具体分析，它依赖于底数的乘法操作，而不再是 $O(\log n)$ 。

---

4. 矩阵快速幂

矩阵快速幂是快速幂的一个拓展，给定一个 $m\times n$ 的矩阵 $A$ ，求它的 $n$ 次幂 $A^n$ 。这也是常见的操作，做法是把矩阵当做变量来操作，程序和上面的很相似。

首先定义矩阵的结构体，并且定义矩阵相乘的操作，注意：矩阵相乘也需要取模。这里则仅仅展示朴素版本的矩阵快速幂：

const int MAXN = 2;		//根据题目要求定义矩阵的阶 const int MOD = 1e4;    //根据题目要求定义模struct Matrix {
   	int m[MAXN][MAXN];	Matrix() {
   		memset(m, 0, sizeof(m));	}}Matrix Multi(const Matrix &a, const Matrix &b) {
    //矩阵乘法	Matrix res;	for (int i = 0; i < MAXN; ++i)		for (int j = 0; j < MAXN; ++j)			for (int k = 0; k < MAXN; ++k)				res.m[i][j] = (res.m[i][j] + a.m[i][k] * b.m[k][j]) % MOD; }Matrix fastPow(const Matrix &a, int n) {
   	Matrix res;	for (int i = 0; i < MAXN; ++i) //初始化为单位矩阵,相当于int res = 1;		res.m[i][i] = 1;	while (n) {
   		if (n & 1)			res = Multi(res, a);		a = Multi(a, a);		n >>= 1;	}	return res;}

朴素的矩阵乘法，得到的矩阵快速幂的复杂度：上面求 $A^n$ ，$A$ 是 $m\times n$ 的方阵，其中矩阵乘法的复杂度是 $O(m^3)$ ，快速幂的复杂度是 $O({\log }_{2}n)$ ，合起来就是 $O(m^3{\log }_{2}n)$ 。应用矩阵快速幂的难点在于，如何把递推关系转换为矩阵。

---

5. 应用题目

求 $n^n$ 的末尾数字，$n\le 10^9$ 快速幂取模、LCM 把递推关系转换为矩阵 有难度的矩阵快速幂 数位DP、矩阵快速幂 AC自动机、矩阵快速幂 和下面一题一样，矩阵快速幂的经典题目，算 Fibonacci 数列。求第 $10^9$ 个 Fibonacci 数，因为直接递推无法完成，所以先用矩阵表示 Fibonacci 数列的递推关系，然后转换为求这个矩阵的 $10^9$ 次幂；

转载地址：https://memcpy0.blog.csdn.net/article/details/108410952 如侵犯您的版权，请留言回复原文章的地址，我们会给您删除此文章，给您带来不便请您谅解！