【矩阵快速幂】- 洛谷 p3390 模板题 （补一道斐波那契数列）

矩阵快速幂

大意：
就是对矩阵求快速幂，将快速幂中res=1换成单位阵，而取模运算完全在我们定义的乘法运算中进行，即Mul，其他的跟快速幂没有区别，当然我现在只掌握了基础，所以做了一道模板题。
题目描述：
题目链接：洛谷-p3390

题目背景
矩阵快速幂

题目描述
给定n*n的矩阵A，求A^k

输入输出格式
输入格式：
第一行，n,k

第2至n+1行，每行n个数，第i+1行第j个数表示矩阵第i行第j列的元素

输出格式：
输出A^k

共n行，每行n个数，第i行第j个数表示矩阵第i行第j列的元素，每个元素模10^9+7

输入输出样例
输入样例#1：
2 1
1 1
1 1
输出样例#1：
1 1
1 1
说明
n<=100, k<=10^12, |矩阵元素|<=1000 算法：矩阵快速幂

就是个板子，不多解释，直接上板子！！！

完整代码：

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string>
#include<cstring>
#include<cmath>
#include<stack>
using namespace std;
#define maxn 1005
typedef long long ll;
const int p=1e9+7;
ll n,k;
typedef struct mat{
ll m[maxn][maxn];
}Matrix;
Matrix Mul(Matrix a,Matrix b){
Matrix ans;
memset(ans.m,0,sizeof ans.m);
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
for(int k=0;k<n;k++){
ans.m[i][j]+=(a.m[i][k]%p)*(b.m[k][j]%p)%p;
ans.m[i][j]%=p;
}
}
}
return ans;
}
Matrix qsm(Matrix a,ll b){
Matrix res;
memset(res.m,0,sizeof res.m);
for(int i=0;i<n;i++) res.m[i][i]=1;
while(b){
if(b&1) res=Mul(res,a);
a=Mul(a,a);
b>>=1;
}
return res;
}
int main(){
Matrix mat,ans;
cin>>n>>k;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++) scanf("%lld",&mat.m[i][j]);
}
ans=qsm(mat,k);
for(int i=0;i<n;i++){
for(int j=0;j<n;j++) printf("%lld ",ans.m[i][j]);
printf("\n");
}
return 0;
}

应用：
矩阵快速幂可用于递推式求解，比如斐波那契数列，我们可以这样来看：
题目：斐波那契数列f(n),给一个n，求f(n)%10000,n<=1e9;

(这题是可以用fibo的循环节去做的，不过这里不讲，反正是水题)

矩阵快速幂是用来求解递推式的，所以第一步先要列出递推式:

f(n)=f(n-1)+f(n-2)

第二步是建立矩阵递推式，找到转移矩阵:

这里还是说一下构建矩阵递推的大致套路,一般An与A(n-1)都是按照原始递推式来构建的，当然可以先猜一个An,主要是利用矩阵乘法凑出矩阵T,第一行一般就是递推式，后面的行就是不需要的项就让与其的相乘系数为0。矩阵T就叫做转移矩阵(一定要是常数矩阵),它能把A(n-1)转移到A(n);然后这就是个等比数列，直接写出通项:此处A1叫初始矩阵。所以用一下矩阵快速幂然后乘上初始矩阵就能得到An,这里An就两个元素(两个位置)，根据自己设置的A(n)对应位置就是对应的值，按照上面矩阵快速幂写法，res[1][1]=f(n)就是我们要求的。

斐波那契数列：

题目分析：
同样是板子题，只不过更加形象的体现除了矩阵快速幂的作用，你直接循环暴力跑这道题…T成🐶，我们这里开二阶矩阵，因为递推式是这样子的：f(n)=f(n-1)+f(n-2)。即含两项，所以用二阶矩阵。

代码展示：

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<string>
#include<cstring>
#include<cmath>
#include<stack>
using namespace std;
#define maxn 1005
typedef long long ll;
const int p=10000;
int n;
typedef struct ma{
ll m[2][2];
}Matrix;
Matrix Mul(Matrix a,Matrix b){
Matrix temp;
memset(temp.m,0,sizeof temp.m);
for(int i=0;i<2;i++){
for(int j=0;j<2;j++){
for(int k=0;k<2;k++){
temp.m[i][j]+=a.m[i][k]*b.m[k][j]%p;
temp.m[i][j]%=p;
}
}
}
return temp;
}
Matrix qsm(Matrix a,int b){
Matrix res;
for(int i=0;i<2;i++){
for(int j=0;j<2;j++){
if(i==j) res.m[i][j]=1;
else res.m[i][j]=0;
}
}
while(b){
if(b&1) res=Mul(res,a);
a=Mul(a,a);
b>>=1;
}
return res;
}
int main(){
while(~scanf("%d",&n)){
if(n==-1) return 0;
if(n==0){
cout<<"0\n";
}
else{
Matrix a,b;
a.m[0][0]=1,a.m[0][1]=1,a.m[1][0]=1,a.m[1][1]=0;
b=qsm(a,n);
cout<<b.m[1][0]<<endl;
}
}
return 0;
}