题目链接: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<cstdio>
#define For(i,L,R) for(int i=L;i<=R;i++)
#define gc getchar
typedef long long ll;
const int N = 105; //根据题意修改
const int M = 105;
const ll mod = 1e9+7;
inline ll read(){
ll a=0;int f=0;char c=getchar();
while(!('0'<=c&&c<='9')){f|=c=='-';c=getchar();}
while(('0'<=c&&c<='9')){a=(a<<3)+(a<<1)+(c^48);c=getchar();}
return f?-a:a;
}
typedef struct Matrix{
int n,m;
ll mat[N][M];
Matrix( int Tn=0,int Tm=0 ):n(Tn),m(Tm) {}
void Clear(){
For(i,0,n) For(j,0,m)
mat[i][j] = 0;
}
}Mat;
Matrix operator + (const Matrix a,const Matrix b){
Matrix c;
int n = a.n , m = b.m;
For(i,0,n-1) For(j,0,m-1){
c.mat[i][j] = a.mat[i][j]+b.mat[i][j];
if( c.mat[i][j] >= mod )
c.mat[i][j] %=mod;
}
return c;
}
Matrix operator * (const Matrix a,const Matrix b){
int n = a.n;
int m = b.m;
ll tmp ;
Matrix c(n,m);
c.Clear();
For(i,0,n-1) For(k,0,m-1) For(j,0,a.m-1){
tmp = a.mat[i][j]*b.mat[j][k] ;
if( tmp >= mod ) tmp %=mod;
c.mat[i][k] = c.mat[i][k]+tmp ;
if( c.mat[i][k] >= mod || c.mat[i][k] < 0 )
c.mat[i][k] = (c.mat[i][k]+mod)%mod;
}
return c;
}
Matrix operator ^ (Matrix a,ll b){
int n = a.n;
Matrix c(n,n);
For( i,0,n-1) For( j,0,n-1) //单位矩阵,主对角线为1
c.mat[i][j] = (i==j) ;
while(b){
if( b&1 )
c = a*c;
a = a * a;
b>>=1;
}
return c;
}
int main(){
ll n,k;
n = read(),k=read();
Matrix a(n,n);
For(i,0,n-1) For(j,0,n-1)
a.mat[i][j] = read();
Matrix b = a ^ k;
For(i,0,n-1) For(j,0,n-1)
printf("%lld%c",b.mat[i][j],j==n-1?'\n':' ');
return 0;
}
P4783 【模板】矩阵求逆
题解
参考了SarvaTathagata大佬写的代码
①找到当前方阵(假设是第k个)的主元
②将主元所在行、列与第k行、第k列交换,并且记录下交换的行和列
③把第k行第k个元素变成它的逆元(同时判无解:逆元不存在)
④更改当前行(乘上刚才那个逆元)
⑤对矩阵中的其他行做和高斯消元几乎完全一样的操作,只有每行的第k列不一样,具体看代码
⑥最后把第②步中的交换逆序地交换回去。
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const ll mod=1e9+7;
ll a[410][410];
int n,row[410],col[410];
void exgcd(int a,int b,int &x,int &y){
if(!b)return x=1,y=0,void();
exgcd(b,a%b,y,x);y-=x*(a/b);
}
int inv(int p){
int x,y;exgcd(p,mod,x,y);
return (x+mod)%mod;
}
void inv(){
for(int k=1;k<=n;k++){
for(int i=k;i<=n;i++) // 找非零主元
for(int j=k;j<=n;j++)if(a[i][j]){
row[k]=i,col[k]=j;break;
}
for(int i=1;i<=n;i++) // 换主元
swap(a[k][i],a[row[k]][i]);
for(int i=1;i<=n;i++)
swap(a[i][k],a[i][col[k]]);
if(!a[k][k]){ //零行,矩阵的秩<n,无解
puts("No Solution");
exit(0);
}
a[k][k]=inv(a[k][k]); //首行元素,求逆元
for(int j=1;j<=n;j++) if(j!=k) //第k行元素,全部乘以逆元
(a[k][j]*=a[k][k])%=mod;
for(int i=1;i<=n;i++) if(i!=k) //把所有除(*,k),(k,*)元素外
for(int j=1;j<=n;j++) if(j!=k) //减去(i,k)*(k,j)
(a[i][j]+=mod-a[i][k]*a[k][j]%mod)%=mod;
for(int i=1;i<=n;i++) if(i!=k) //第k列元素,全部乘以逆元
a[i][k]=(mod-a[i][k]*a[k][k]%mod)%mod;
}
for(int k=n;k;k--){ // 把以前换出去,换回来
for(int i=1;i<=n;i++)
swap(a[col[k]][i],a[k][i]);
for(int i=1;i<=n;i++)
swap(a[i][row[k]],a[i][k]);
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
scanf("%lld",&a[i][j]);
inv();
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
printf("%lld%c",a[i][j],j==n?'\n':' ');
return 0;
}