矩阵模板（矩阵快速幂算法用法）

矩阵快速幂用法

应用矩阵解线性方程组入门视频

用于计算递归的方程组等：如f(n)=f(n-1)+f(n-2)+f(n-3);

在引入矩阵快速幂之前
先引入问题：
想必大家都不陌生斐波那切数.
数字公式表达：f(n)= f(n-1)+ f(n-2) , f(1)=1; f(0)=0; n>=2;
嗯很美观~ 然后呢？
假如n很大呢， 比如 1e13 呢? 嗯哼? 能在一秒呢计算出吗。
显然不可能呢~

来–上矩阵，解决。
先别急；
来个简单例子 温习下线性方程组：
2x+3y=z
x-y=2z

这里x，y，z均是未知数
那么用矩阵表示:
$$A=\left[\begin{array}{cc}2& 3\\ 1& -1\end{array}\right]\text{\hspace{1em}(2)}$$
$$X=\left[\begin{array}{c}x\\ y\end{array}\right]\text{\hspace{1em}(2)}$$
A*X 得到 $$结果=\left[\begin{array}{c}z\\ 2z\end{array}\right]\text{\hspace{1em}(2)}$$
来—上公式
f(n-1) +f(n-2)=f(n);
f(n-1)=f(n-1);
这里我们暂且可以把 等式右边的f(n-1), f(n-2) 看成未知数.

哇哇哇~ 这不就轻轻松松求出来了， 现在的关键是要知道矩阵A的(n-2)次方，这简单 ，快速幂就好啦.
回到最前面的问题：n= 1e13；log(n)时间内即可得到结果。

上代码：

#include<bits/stdc++.h>
#define ll long long
#define rep(i,a,b) for(int i=a;i<=b;i++)
ll gcd(ll a,ll b){ return b? gcd(b,a%b):a;}
const int N=2e5+10;
const ll mod=1e9+7;
ll read(){
    ll s = 0, f = 1; char ch = getchar();
    while(!isdigit(ch)){
        if(ch == '-') f = -1;
        ch = getchar();
    }
    while(isdigit(ch)) s = (s << 3) + (s << 1) + (ch ^ 48), ch = getchar();
    return s * f;
}
using namespace std;
ll n,m;
struct Matrix{
	ll  mat[100][100];
	int n, m;
	Matrix(){
		n = m = 20;
		memset(mat, 0, sizeof(mat));
	}
	//重新定义矩阵的大小 
	void init(int row,int col){
		n = row;
		m = col;
	}
	//单位矩阵 
	void init_e(){
		rep(i,1,n){
			rep(j,1,m){
				mat[i][j] = (i == j);
			}
		}
	}
	//打印矩阵 
	void print(){
		rep(i,1,n){
			rep(j,1,m){
				cout<<mat[i][j]<<" ";
			}
			cout<<endl;
		}
		cout<<endl;
	}
};
// 矩阵加法 
Matrix operator +(Matrix a,Matrix b){
	Matrix ret;
	ret.init(a.n,a.m);
	rep(i,1,n){
		rep(j,1,m){
			ret.mat[i][j] = a.mat[i][j] + b.mat[i][j];
		}
	}
	return ret;
}
// 矩阵乘法 
Matrix  operator *(Matrix a,Matrix b){
	Matrix ret;
	ret.init(a.n, b.m);
	rep(i,1,a.n){
		rep(j,1,b.m){
			rep(k,1,a.m){
				ret.mat[i][j] = (ret.mat[i][j]+a.mat[i][k] * b.mat[k][j])%mod;
			}
		}
	}
	return ret;
}
// 矩阵快速幂  求递归方程 
Matrix operator ^ (Matrix a,ll b){
	// n X n
	Matrix sum = a;
	//sum.init(a.n, a.m);
	sum.init_e();
	//a=a*a;
	//return a;
	while(b){
		if(b&1){
			sum = sum * a;
		}
		a = a * a;
		b = b >> 1;
	}
	return sum;
}
void solve(){
	n=read();
	if(n==1||n==2){
		cout<<n-1<<endl;
		return ;
	}
	Matrix a;
	a.init(2,2);
	a.mat[1][1]=1; a.mat[1][2]=1;
	a.mat[2][1]=1; a.mat[2][2]=0;
	a=a^(n-2);
	//cout<<"dsa"<<endl;
	Matrix b;
	b.init(2,1); b.mat[1][1]=1; b.mat[2][1]=0;
	a=a*b;
	cout<<a.mat[1][1]<<endl;
}
int main (){
  freopen("in.txt","r",stdin);
  freopen("out.txt","w",stdout);
  int T; cin>>T;
  while(T--){
  solve();
  }
  return 0;
}
/*
in:
5
1
3
10
100
1000000000000000000
out:
0
1
34
94208912
470273943
*/

A Simple Math Problem
传送门
题意很简单：
就是求 $f(m)$;
然后给出了：
$f(x)$ 关系表达式：
$ifx<10f(x)=x$.
$ifx>=10f(x)=a0\ast f(x-1)+a1\ast f(x-2)+a2\ast f(x-3)+\dots \dots +a9\ast f(x-10);$
组成是个方程组就OK了；
$f(x)=a0\ast f(x-1)+a1\ast f(x-2)+a2\ast f(x-3)+\dots \dots +a9\ast f(x-10)$
$f(x-1)=1\ast f(x-1)$
$f(x-2)=1\ast f(x-2)$
$f(x-3)=1\ast f(x-3)$
$f(x-4)=1\ast f(x-4)$
$f(x-5)=1\ast f(x-5)$
$f(x-6)=1\ast f(x-6)$
$f(x-7)=1\ast f(x-7)$
$f(x-8)=1\ast f(x-8)$
$f(x-9)=1\ast f(x-9)$
系数矩阵A表示：
$$X=\left[\begin{array}{cccccccccc}a0& a1& a2& a3& a4& a5& a6& a7& a8& a9\\ 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\end{array}\right]\text{\hspace{1em}(2)}$$
那就很清楚了
代码如下：

#include<bits/stdc++.h>
#define ll long long
#define rep(i,a,b) for(int i=a;i<=b;i++)
ll gcd(ll a,ll b){ return b? gcd(b,a%b):a;}
const int N=2e5+10;
ll read(){
    ll s = 0, f = 1; char ch = getchar();
    while(!isdigit(ch)){
        if(ch == '-') f = -1;
        ch = getchar();
    }
    while(isdigit(ch)) s = (s << 3) + (s << 1) + (ch ^ 48), ch = getchar();
    return s * f;
}
using namespace std;
ll n,m;
ll mod;
struct Matrix{
	ll  mat[100][100];
	int n, m;
	Matrix(){
		n = m = 20;
		memset(mat, 0, sizeof(mat));
	}
	void init(int row,int col){
		n = row;
		m = col;
	}
	void init_e(){
		rep(i,1,n){
			rep(j,1,m){
				mat[i][j] = (i == j);
			}
		}
	}
	void print(){
		rep(i,1,n){
			rep(j,1,m){
				cout<<mat[i][j]<<" ";
			}
			cout<<endl;
		}
		cout<<endl;
	}
};
Matrix operator +(Matrix a,Matrix b){
	Matrix ret;
	ret.init(a.n,a.m);
	rep(i,1,n){
		rep(j,1,m){
			ret.mat[i][j] = a.mat[i][j] + b.mat[i][j];
		}
	}
	return ret;
}
Matrix  operator *(Matrix a,Matrix b){
	Matrix ret;
	ret.init(a.n, b.m);
	rep(i,1,a.n){
		rep(j,1,b.m){
			rep(k,1,a.m){
				ret.mat[i][j] = (ret.mat[i][j]+a.mat[i][k] * b.mat[k][j])%mod;
			}
		}
	}
	return ret;
}
Matrix operator ^ (Matrix a,int b){
	// n X n
	Matrix sum = a;
	//sum.init(a.n, a.m);
	sum.init_e();
	//a=a*a;
	//return a;
	while(b){
		if(b&1){
			sum = sum * a;
		}
		a = a * a;
		b = b >> 1;
	}
	return sum;
}
 ll k;
void solve(){
	Matrix a;
//    int n,m;
    //ll k;
    //k=read(); mod=read();
    if(k<=9){
        ll ans;
    	for(int i=1;i<=m;i++) ans=read();
    	cout<<k%mod<<endl;
	}else{
		Matrix a;
		a.init(10,10);
	    rep(i,1,10) a.mat[1][i]=read();
	    rep(i,2,10) a.mat[i][i-1]=1;
	    a=a^(k-9);
	    //a.print();
	    Matrix b;
	    b.init(10,1);
	    for(int i=1;i<=10;i++) b.mat[i][1]=10-i;
	    a=a*b;
	    cout<<a.mat[1][1]<<endl;
	}
}
int main (){
  //freopen("in.txt","r",stdin);
  //freopen("out.txt","w",stdout);
  //int T; cin>>T;
  while(scanf("%d %d",&k,&mod)!=EOF){
  solve();
  }
  return 0;
}

又一道经典的矩阵快速幂题：

PS： 其实也是动态规划题。
552. 学生出勤记录 II
--------------------------------- 先鸽着…以后补doge---------------------------------------------------
参考：
起风了_唯有努力生存
疯狂奔跑的少年