何为快速幂
快速幂（Exponentiation by squaring，平方求幂）是一种简单而有效的小算法，它可以在 $O(log_{2}{n})$ 的时间复杂度计算乘方。快速幂不仅本身非常常见，而且后续很多算法也都会用到快速幂。
让我们来思考一个问题：$7^{10}$ 怎么运算比较快？
- 方法1：最朴素的想法，7×7=49，49×7=434，……一步步算，这样一共进行了9次乘法
这样算无疑太慢了，对于高效的CPU来说，每次运算都只乘一个个位数，很浪费性能，时间复杂度为O(n)。这时我们可以想，也许可以拆分问题
- 方法2：先算 $7^5$ ，即 7×7×7×7×7 ，再算它的平方，这样一共进行了5次乘法
但这样也并不是最优解，因为显然可以发现，对于 $7^5$ ，我们仍然可以拆分问题
- 方法3：先算 7×7 = 49 ，则 $7^5$ = 49×49×7，再算它的平方，这样一共进行了4次乘法
模仿这样的拆分过程，我们就得到了一个在 $O(log_{2}{n})$ 的时间内算出幂的方法，也就是快速幂的根源
快速幂实现
一、递归快速幂
根据上文提出的问题解决方式，我们可以发现这无非是一个二分的思路，因此我们很自然地可以得到一个递归方程：
计算 $a^n$ ，如果 n 是偶数（且不为0），那么就先计算 $a^{\frac{n}{2}}$ ，然后平方；如果 n 是奇数，那么就先计算$a^{n-1}$，再乘上 a ；递归出口是 $a^0 = 1$ 。可以看到思路很自然，代码也很好写
// 递归快速幂
int qpow(int a, int n)
{
if (n == 0)
return 1;
else if (n % 2 == 1)
return qpow(a, n - 1) * a;
else
{
int tmp = qpow(a, n / 2);
return tmp * tmp;
}
}注意，这个 tmp 变量是必要的，因为如果不把 $a^{\frac{n}{2}}$ 记录下来，直接写成 qpow(a, n /2)*qpow(a, n /2) ，那会计算两次 $a^{\frac{n}{2}}$ ，整个算法就退化为 O(n)复杂度
在实际解决问题中，题目常常会要求对一个大素数取模，这是因为计算结果可能很大，但是没必要在快速幂中考察高精度，此时我们的快速幂也应当进行取模。要注意，原则是步步取模，如果 mod 较大，应该把变量改为 long long
//递归快速幂（对大素数取模）
typedef long long ll;
const ll mod = 1e9+7;
ll qpow(ll a, ll n)
{
if (n == 0)
return 1;
else if (n % 2 == 1)
return qpow(a, n - 1) * a % mod;
else
{
ll tmp = qpow(a, n / 2) % mod;
return tmp * tmp % mod;
}
}虽然递归的代码简单易懂，但是我们也知道递归会产生额外的空间开销，递归层数过多容易导致爆栈。我们可以把递归改写为循环，来避免对栈空间的大量占用，也就是非递归快速幂。
二、非递归快速幂
根据知识我们可知，计算机是通过二进制进行存储数据的，对二进制数据可以直接进行操作，因此我们换一个角度来引入非递归的快速幂。还是 $7^{10}$ ，但这次我们把 10 写成二进制的形式，也就是 $(1010)_2$
现在我们要计算 $7^{(1010)_2}$ ，可以怎么做？我们很自然地会想到把它拆分为 $7^{(1000)_2}$ × $7^{(10)_2}$ ，即 $7^8$ × $7^2$ = $7^{8+2}$ = $7^{10}$ 。同样我们也可以发现，对于任意的整数，我们都可以把它拆成若干个 $7^{(10……)_2}$ 的形式相乘。而这些 $7^{(10……)_2}$ ，恰好就是 $7^1$ 、$7^2$、$7^4$……我们只需不断把底数平方就可以算出它们。
下面给出代码，然后再来一步步解释意思：
// 非递归快速幂
int qpow(int a, int n)
{
int ans = 1; // 存储答案
int tmp = a; // 存储底数,用于自乘(也可直接用a)
while(n) {
if(n&1) ans = ans * tmp; // 如果n的当前末位为1,ans乘上当前的tmp
tmp = tmp * tmp; // 每次循环,底数都要平方,即tmp自乘
n >>= 1; // n往右移一位
}
return ans;
} 最初设置 ans = 1，然后我们对指数10的二进制 $(1010)_2$ 来一位一位算：
- 1010的最后一位是0，所以 $a^1$ 不需要计入ans，跳过，然后 1010 右移一位变为 101，a 自乘变为 $a^2$。
- 101的最后一位是1，所以 $a^2$ 需要计入ans，ans乘上 $a^2$ ，然后 101 右移一位变为 10，$a^2$ 自乘变为 $a^4$。
- 10的最后一位是0，所以 $a^4$ 不需要计入ans，跳过，然后 10 右移一位变为 1， $a^4$ 自乘变为 $a^8$。
- 1的最后一位是1，所以 $a^8$ 需要计入ans，ans乘上 $a^8$ 。循环结束，返回结果 ans 。
可以看出，非递归快速幂的时间复杂度即为指数 n 的二进制表示的数位长度，即 $O(log_{2}{n})$ 。同样的，非递归快速幂也有取模的版本，代码如下：
// 非递归快速幂（对大素数取模）
typedef long long ll;
const ll mod = 1e9+7;
ll qpow(ll a, ll n)
{
ll ans = 1; // 存储答案
ll tmp = a % mod; // 存储底数,用于自乘(也可直接用a)
while(n) {
if(n&1) ans = (ans * tmp) % mod; // 如果n的当前末位为1,ans乘上当前的tmp
tmp = (tmp * tmp) % mod; // 每次循环,底数都要平方,即tmp自乘
n >>= 1; // n往右移一位
}
return ans;
} 代码中的 >> 是右移位运算符，表示把二进制数往右移一位，相当于/2；&是按位与，&1可以理解为判断二进制数最后一位是否为1，相当于%2==1。可以看到虽然非递归快速幂因为牵扯到二进制理解起来稍微复杂一点，但基本思路其实和递归快速幂没有太大的出入。
快速幂拓展
一、广义快速幂
上面提及的都是整数的快速幂，但其实，在算 $a^n$ 时，只要 a 的数据类型支持乘法(或者说有某种二元运算)且满足结合律，快速幂的算法都是有效的，这个就是广义快速幂的思维。以下给出广义快速幂的一些定义（涉及一点离散数学的知识）。
已知群有4大性质：1.运算封闭性，2.满足结合律，3.有幺元，4.有逆元，设：○ 为一种二元运算且与集合 V 构成群，a∈V，e 为 ○ 运算的幺元。幺元 e 满足对于任意的 a ,有 e ○ a = a ○ e = a ，我们可以记 $a^0$ = e，$a^n$ = $a^{n-1}$ ○ a ，则有以下性质 $a^{n+m}$ = $a^n$ ○ $a^m$ 。则此时计算 a 关于 ○ 运算的 n 次幂的快速幂可以这样写：
//泛型的非递归快速幂
template <typename T>
T qpow(T a, ll n)
{
T ans = 1; // 赋值为运算单位元，可能要根据构造函数修改
T tmp = a;
while (n)
{
if (n & 1)
ans = ans * tmp; // 这里就最好别用自乘了，不然重载完*还要重载*=，有点麻烦。
tmp = tmp * tmp;
n >>= 1;
}
return ans;
}例如 $a^n$ = a×a×a×a×a… 是关于乘法的幂运算，又因为 1×a = a×1 = a,所以乘法幺元 e 就是 1 ，带入上面的程序就可以得到最常见的乘法快速幂。同理，a×b = a+a+a+a……，其实就是关于加法的幂运算，至于加法的幺元，因为 a+0 = 0+a = a，所以加法幺元就是0，带入上面程序便可得到快速乘法运算。
一般来说此时的泛型 T 是我们自己定义的一些数据类型，因此广义快速幂代码中的乘法并不是普通四则运算中的乘法，而是自定义数据类型中重载的乘法运算，所以一般广义快速幂还会自己定义数据类型，代码结构如下（可自行修改）
// 自定义数据类型
struct DataType
{
int e; // 运算单位元
DataType(int e) : e(e) {}
// 重载 “乘法” 运算
DataType operator*(const DataType &y)
{
DataType ans;
/*
自定义运算内容
*/
return ans;
}
};要注意的是，较复杂类型的快速幂的时间复杂度并不一定是 $O(log_{2}{n})$ ，而是取决于底数的数据类型的 “乘法” 运算的时间复杂度。
二、矩阵快速幂
矩阵快速幂就是一种典型的广义快速幂，它可以加快矩阵的幂运算，同理高精度整数也有快速幂，这里先介绍矩阵快速幂。矩阵快速幂最为经典的用法就是求斐波那契数列的第 n 项（当n特别大的时候）
题目背景
大家都知道，斐波那契数列是满足如下性质的一个数列：
当 n <= 2时，$F_n$ = 1 ; 当 n >= 3 时，$F_n$ = $F_{n-1}$ + $F_{n-2}$
题目描述
请你求出 $F_n$ mod ${10}^{9}$+7 的值
数据范围
对于 100% 的数据，1 <= n <= $2^{63}$
很明显，题目给的 n 过于大，哪怕是递推去求也会时间溢出，而且这么大的 n 也直接杜绝了打表水过去的可能性，因此要用矩阵快速幂去优化。矩阵快速幂适用于类似斐波那契数列的问题，它们一般都有个特点：各项之间的关系由一个递推式来确定；这种题的难点就在于找出递推式并把它转换成对应的变换矩阵。
由上面的过程可知，我们需要做的是找到如 F(1) = 1, F(2) = 1 这样的初始情况，构成初始矩阵；然后根据递推式得出变换矩阵 A ，用于求幂，最后我们算出 $A^{n}$ ，然后用初始矩阵乘这个结果就可以得出最终答案。这样我们就把原本较为复杂的问题转换成了求某个矩阵的幂的问题了。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9+7; // 取模
const int len = 2; // 矩阵维度
// 自定义矩阵
struct Matrix
{
// 矩阵
ll a[len][len];
// 构造函数
Matrix(ll a11 = 0, ll a12 = 0, ll a21 = 0, ll a22 = 0) {
a[0][0] = a11; a[0][1] = a12;
a[1][0] = a21; a[1][1] = a22;
}
// 重载乘法运算(在这里进行取模操作)
Matrix operator*(const Matrix &B) {
Matrix ans;
for(int i = 0; i < len; ++i) {
for(int j = 0; j < len; ++j) {
for(int k = 0; k < len; ++k) {
// 步步取模
ans.a[i][j] = (ans.a[i][j] + ((a[i][k] * B.a[k][j])%mod))%mod;
}
}
}
return ans;
}
};
// 矩阵快速幂
Matrix qpow(Matrix a, ll n)
{
Matrix ans(1,0,0,1); // 单位矩阵
while(n) {
if(n&1) ans = ans * a;
a = a * a;
n >>= 1;
}
return ans;
}
int main()
{
Matrix A(0,1,1,1); // 递推式的变换矩阵A
ll n;
scanf("%lld", &n);
Matrix ans = qpow(A, n-1); // 如果初始条件用F(1)=1,F(2)=1 则只需要计算到n-1次幂
printf("%lld\n", (ans.a[0][0] + ans.a[0][1])%mod);
return 0;
}参考资料
版权属于：PCsky
本文链接：http://www.hyouka.club/index.php/archives/137/
转载时须注明出处及本声明