【洛谷P3390】 【模板】矩阵快速幂

题目链接：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;
}