基础算法—快速幂详解
幂运算是非常常见的一种运算，求取$a^n$,最容易想到的方法便是通过循环逐个累乘，其复杂度为$O(n)$，这在很多时候是不够快的，所以我们需要一种算法来优化幂运算的过程。
一、快速幂——反复平方法
该怎样去加速幂运算的过程呢？既然我们觉得将幂运算分为n步进行太慢，那我们就要想办法减少步骤，把其中的某一部分合成一步来进行。
比如，如果$n$能被2整除，那我们可以先计算一半，得到$a^{n/2}$的值，再把这个值平方得出结果。这样做虽然有优化，但优化的程度很小，仍是线性的复杂度。
再比如，如果我们能找到$2^k = n$，那我们就能把原来的运算优化成$((a^2)^2)^2...$，只需要$k$次运算就可以完成，效率大大提升。可惜的是，这种条件显然太苛刻了，适用范围很小。不过这给了我们一种思路，虽然我们很难找到$2^k = n$，但我们能够找到$2^{k_1} + 2^{k_2} + 2^{k_3} +......+ 2^{k_m} = n$。这样，我们可以通过递推，在很短的时间内求出各个项的值。
我们都学习过进制与进制的转换，知道一个$b$进制数的值可以表示为各个数位的值与权值之积的总和。比如，2进制数$1001$，它的值可以表示为10进制的$1\times2^3 + 0\times2^2 + 0\times2^1 + 1\times2^0$，即$9$。这完美地符合了上面的要求。可以通过2进制来把$n$转化成$2^{k_m}$的序列之和，而2进制中第$i$位（从右边开始计数，值为$1$或是$0$）则标记了对应的$2^{i - 1}$是否存在于序列之中。譬如，$13$为二进制的$1101$，他可以表示为$2^3 + 2^2 + 2^0$，其中由于第二位为$0$，$2^1$项被舍去。
如此一来，我们只需要计算$a、a^2、a^4、a^8......a^{2^{k_m}}$的值（这个序列中的项不一定都存在，由$n$的二进制决定）并把它们乘起来即可完成整个幂运算。借助位运算的操作，可以很方便地实现这一算法，其复杂度为$O(\log n)$。
typedef long long ll; ll mod; ll qpow(ll a, ll n)//计算a^n % mod { ll re = 1; while(n) { if(n & 1)//判断n的最后一位是否为1 re = (re * a) % mod; n >>= 1;//舍去n的最后一位 a = (a * a) % mod;//将a平方 } return re % mod; }
取模运算一般情况下是需要的，当然也可以省去。
二、矩阵快速幂
快速幂只是通过二进制拆分$n$来加速幂运算的手段，当然并不只适用于求取数字的幂次，对于矩阵的$n$次方，也可以用同样的手段求取。除了乘法的规则与上面的快速幂不同之外，其他方面并没有太大的差别。
不过这有什么意义呢？利用矩阵的幂次，我们可以快速地完成递推。
比如，在POJ3070 Fibonacci中，就需要我们快速地求取斐波那契数列的第$n$项（取模），对于$n=10^{18}$，一步步推过去显然太慢了，那么我们可以考虑构造矩阵来帮我们完成递推。
首先复习一下矩阵的乘法:
对于$i \ge 3$
有$Fib_i = Fib_{i-1} + Fib_{i-2}$
可以构造出矩阵递推式$$\left[\begin{array}{ll}{1} & {1} \\ {1} & {0}\end{array}\right]\left[\begin{array}{c}{F i b_{i}} \\ {F i b_{i-1}}\end{array}\right]=\left[\begin{array}{c}{F i b_{i+1}} \\ {F i b_{i}}\end{array}\right]$$
我们就可以利用矩阵快速幂以$O(\log n)$求取斐波那契数列的第$n$项了。
const ll mod = 10000; const int maxv = 2; struct Matrix { ll a[maxv][maxv]; //矩阵 Matrix operator*(const Matrix &b) const& { //矩阵乘法，复杂度O(maxv^3)，也可看作常数，但maxv较大（大于5）时会使运算时间提高好几个数量级 Matrix ans; for (int i = 0; i < maxv; ++i) { for (int j = 0; j < maxv; ++j) { ans.a[i][j] = 0; for (int k = 0; k < maxv; ++k) { ans.a[i][j] += a[i][k] * b.a[k][j] % mod; ans.a[i][j] %= mod; } } } return ans; } static Matrix qpow(Matrix x, ll n) {//矩阵快速幂，将乘法复杂度看作常数则复杂度为O(log n) Matrix ans; for (int i = 0; i < maxv; ++i) { for (int j = 0; j < maxv; ++j) { if (i == j) ans.a[i][j] = 1; else ans.a[i][j] = 0; } }//初始化为单位矩阵，参考普通数字的快速幂这里初始化为1 while (n) {//其余细节基本相同 if (n & 1) ans = ans * x; x = x * x; n >>= 1; } return ans; } Matrix(ll temp[maxv][maxv]) {//构造方法 for (int i = 0; i < maxv; ++i) { for (int j = 0; j < maxv; ++j) { a[i][j] = temp[i][j]; } } } Matrix() { } }; ll fib(ll n) {//求取斐波那契数列第n项（本题取模） if (n == 0) return 0; if (n <= 2) return 1; ll temp[maxv][maxv] = { 1, 1, 1, 0 }; Matrix m(temp); m = Matrix::qpow(m, n - 2); return (m.a[0][1] + m.a[0][0]) % mod; }
当然了，构造矩阵的方法是多种多样的，例如上面的题目中就给出了另一种构造出斐波那契数列的方法。
矩阵快速幂能解决的递推也远不止斐波那契数列，下面根据我个人的习惯列举几个比较常见、比较“套路”的情况：
（1）对于$a_n = x\cdot a_{n-1} + y\cdot a_{n-2} + c^n$，可以得到
$$
\left[\begin{array}{ccc}{x} & {y} & {c} \\ {1} & {0} & {0} \\ {0} & {0} & {c}\end{array}\right]\left[\begin{array}{c}{a_{i}} \\ {a_{i-1}} \\ {c^{i}}\end{array}\right]=\left[\begin{array}{c}{a_{i+1}} \\ {a_{i}} \\ {c^{i+1}}\end{array}\right]
$$
（2）对于$a_n = a_{n-1} + a_{n-2} + n^2$，可以得到
$$
\left[\begin{array}{ccccc}{1} & {1} & {1} & {2} & {1} \\ {1} & {0} & {0} & {0} & {0} \\ {0} & {0} & {1} & {2} & {1} \\ {0} & {0} & {0} & {1} & {1} \\ {0} & {0} & {0} & {0} & {1}\end{array}\right]\left[\begin{array}{c}{a_{i}} \\ {a_{i-1}} \\ {i^{2}} \\ {i} \\ {1}\end{array}\right]=\left[\begin{array}{c}{a_{i+1}} \\ {a_{i}} \\ {(i+1)^{2}} \\ {i+1} \\ {1}\end{array}\right]
$$
（3）对于$a_n=x\cdot a_{n-1} + y\cdot a_{n-2}$，求$a_n$的前缀和$s_n$
$$
\left[\begin{array}{lll}{1} & {x} & {y} \\ {0} & {x} & {y} \\ {0} & {1} & {0}\end{array}\right]\left[\begin{array}{c}{s_{i}} \\ {a_{i}} \\ {a_{i-1}}\end{array}\right]=\left[\begin{array}{c}{s_{i+1}} \\ {a_{i+1}} \\ {a_{i}}\end{array}\right]
$$
