第？课——基于矩阵快速幂的递推解法
由于中间的数论部分我自己学的很差，没有办法写出清晰的博客来，所以这里跳过了数论部分的博客，来到矩阵快速幂。
递推
递推是一个非常常用的工具。比如经典的斐波那契数列：
\[f(x)= \left\{ \begin{array}{**lr**} 1 &, 0\leq x\leq 1 \\ f(x-1)+f(x-2)&, 2 \leq x \\ \end{array} \right. \]很明显，\(f(n)\)依赖于\(f(n-1)\)和\(f(n-2)\)，所以我们需要先计算\(f(n-1)f(n-2)\)才能计算\(f(n)\)。
假设我们现在需要求f(1e9+7)，你很快发现了，这个数字非常大。所以我们要求只需要结果对\(MOD\)取模就好了，而\(MOD=1e9+7\)。问题是：我们的迭代算法是\(O(n)\)的。那如何快速的求解这样一个递推问题呢？
矩阵与递推的联系
让我们站在巨人的肩膀上来看这个递推问题的第二项以及之后：
\[\begin{pmatrix}f(n)\\? \\\end{pmatrix}= \begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\ \begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix} \]数学家们把这个问题用矩阵的形式表现了出来。但是矩阵上还有一行是空的。打个比方，我们在求\(f(2)\)的时候，矩阵可以写成：
\[\begin{pmatrix}f(2)\\? \\\end{pmatrix}= \begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\ \begin{pmatrix}f(1)\\f(0)\end{pmatrix} \]那如果我们需要继续计算\(f(3)\)呢？
\[\begin{pmatrix}f(3)\\? \\\end{pmatrix}= \begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\ \begin{pmatrix}f(2)\\f(1)\end{pmatrix} \]我们需要知道
\[\begin{pmatrix}f(2)\\f(1)\end{pmatrix} \]可是\(f(2)f(1)\)哪里来呢？不妨把我们的通项改写一下，在计算\(f(n)\)的同时顺带计算\(f(n-1)\)？
\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}= \begin{pmatrix}1&1\\1&0\\\end{pmatrix}\dot\\ \begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix} \]！！wow，我们得到了一个强大的新递推式子。这时候有人要不乐意了，这不是复杂化了么，我们本来只要两个数加一加，现在还要算一个2$\times$4的矩阵乘法？可是，我们的矩阵是常数。 再稍微改写一下这个式子：
\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}= \begin{pmatrix}1&1\\1&0\\\end{pmatrix}^{n-1} \\ \dot\\ \begin{pmatrix}f(1)\\f(0)\end{pmatrix} \]相信你看到这里已经顿悟了，因为这个矩阵的高次幂我们可以大做文章！
快速幂与矩阵快速幂
快速幂
快速幂很简单，这里直接给出代码：
int quickpow(int a,int b,int x) {
while(x>0) {
if(x&1) a=a*b;
b=b*b;
x>>=1;
}
return a;
}
快速幂原理给一个友情链接：
快速幂总结
矩阵乘法
先放一个市面上常见的通用矩乘：
for(int i=1; i<=n; i++) {
for(int j=1; j<=n; j++) {
for(int k=1; k<=n; k++) {
tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod; \\拖慢速度
tmp.mat[i][j]%=mod;
}
}
}
这个矩阵乘非常符合人的想法，但是有一个点让它比较慢，我们矩阵快速幂只需要进行幂次乘法，可以做一些优化：
for(int i=1; i<=n; i++) {
for(int j=1; j<=n; j++) {
for(int k=1; k<=n; k++) {
tmp.mat[i][k]+=(a.mat[i][j]%mod*b.mat[j][k]%mod)%mod;
tmp.mat[i][k]%=mod;
}
}
}
看出不一样的地方了吗？唯一的区别就在下标处。把k放在第一维会大大增加寻址的计算量，所以需要把k放在第二维上就能减少非常多的计算量。其实矩阵小的时候也没什么区别
矩阵快速幂
其实和普通快速幂没有多大区别，只需要重载一个乘法运算符就好了。这里给一份结构体加快速幂的代码：
#include <iostream>
#include <string.h>
using namespace std;
const int mod = (int)1e9+7;
struct MyMat{
MyMat(int n_=24,int m_=24):n(n_),m(m_){
memset(mat,0,sizeof(mat));
}
int mat[25][25];
int m,n;
friend MyMat operator*(const MyMat&a,const MyMat&b){
MyMat tmp;
if(a.m!=b.n)throw("Matrix size mismatch");
for(int i=1; i<=a.n; i++) {
for(int j=1; j<=b.m; j++) {
for(int k=1; k<=a.m; k++) {
tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod;
tmp.mat[i][j]%=mod;
}
}
}
tmp.n = a.n,tmp.m=b.m;
return tmp;
}
friend istream& operator>>(istream&is,MyMat& M){
for(int i=1;i<=M.n;i++){
for(int j=1;j<=M.m;j++){
is>>M.mat[i][j];
}
}
return is;
}
friend ostream& operator<<(ostream&os,const MyMat& M){
for(int i=1;i<=M.n;i++){
for(int j=1;j<=M.m;j++){
if(j!=1)os<<" ";
os<<M.mat[i][j];
}
os<<"\n";
}
return os;
}
};
MyMat quickpow(MyMat a,MyMat b,int x) {
while(x>0) {
if(x&1) a=b*a;
b=b*b;
x>>=1;
}
return a;
}
int main(){
MyMat a(2,3),b(3,2);
cin>>a>>b;
auto ans = a*b;
cout<<"a:\n"<<a<<"*\nb:\n"<<b<<"=\n"<<ans;
ans = quickpow(a,ans,5);
cout<<"a*ans^5=\n"<<ans;
}
感兴趣的可以当一个参考。下面是本地运行测试输出
g++ -o test test.cpp
./test
1 2 1 2 1 2
1 2 3 1 2 3
a:
1 2 1
2 1 2
*
b:
1 2
3 1
2 3
=
9 7
9 11
a*ans^5=
2480048 2480080 2480048
3188656 3188624 3188656
解题
斐波那契数列
让我们回到经典的斐波那契数列，写出这个矩阵递推式后，加上矩阵快速幂
\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}= \begin{pmatrix}1&1\\1&0\\\end{pmatrix}^{n-1} \\ \dot\\ \begin{pmatrix}f(1)\\f(0)\end{pmatrix} \]我们就能快速地写出斐波那契数列的代码了，和上面不同的地方只有quickpow函数的x改为了long long：
MyMat quickpow(MyMat a,MyMat b,long long x) {
while(x>0) {
if(x&1) a=b*a;
b=b*b;
x>>=1;
}
return a;
}
int main(){
long long n;
cin>>n;
if(n<2)cout<<"1"<<endl;
else {
MyMat A(2,2);
A.mat[1][1]=1;
A.mat[1][2]=1;
A.mat[2][1]=1;
MyMat x(2,1);
x.mat[1][1]=1;
x.mat[2][1]=1;
auto ans = quickpow(x,A,n);
cout<<ans.mat[1][1]<<endl;
}
return 0;
}
A Simple Math Problem
重点：
这个矩阵我就不写了，自己动手尝试一下吧～
Count
重点：
//AC代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=123456789;
int t;
ll n;
struct Matrix {
ll mat[25][25];
} pre,A;;
Matrix operator * (Matrix a,Matrix b) {
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
for(int i=1; i<=6; i++) {
for(int j=1; j<=6; j++) {
for(int k=1; k<=6; k++) {
tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod;
tmp.mat[i][j]%=mod;
}
}
}
return tmp;
}
Matrix quickpow(Matrix a,Matrix b,ll x) { //A*B^x
while(x>0) {
if(x&1) a=a*b;
b=b*b;
x>>=1;
}
return a;
}
void init() {
memset(pre.mat,0,sizeof(pre.mat));
memset(A.mat,0,sizeof(A.mat));
//1 2 27 9 3 1
pre.mat[1][1]=1;
pre.mat[1][2]=2;
pre.mat[1][3]=27;
pre.mat[1][4]=9;
pre.mat[1][5]=3;
pre.mat[1][6]=1;
A.mat[1][2]=2;
A.mat[2][1]=1;
A.mat[2][2]=1;
A.mat[3][2]=1;
A.mat[3][3]=1;
A.mat[4][3]=3;
A.mat[4][4]=1;
A.mat[5][3]=3;
A.mat[5][4]=2;
A.mat[5][5]=1;
A.mat[6][3]=1;
A.mat[6][4]=1;
A.mat[6][5]=1;
A.mat[6][6]=1;
}
int main() {
init();
scanf("%d",&t);
while(t--) {
scanf("%lld",&n);
Matrix ans=quickpow(pre,A,n-1);
printf("%lld\n",ans.mat[1][1]);
}
return 0;
}