1. 程式人生 > >P4783 【模板】矩陣求逆

P4783 【模板】矩陣求逆

傳送門

線性代數真的好珂怕……以下如果有漏洞歡迎指出

定義矩陣的三種初等行變換:

1.交換某兩行

2.將某一行的所有元素乘上\(k\)(\(k\neq 0\))

3.將某一行的所有元素乘上\(k\)加到另一行去

每一個初等變換都對應一個初等矩陣,即矩陣\(A\)做某一線性變換等價於用一個對應的初等矩陣左乘\(A\)。若有一堆初等變換\(1,2,...l\),對應的初等矩陣分別為\(P_1,P_2,...,P_l\),那麼經過這些線性變換後的矩陣即為\(P_l....P_2P_1A=PA\)\(P\)為之前那堆東西的乘積)

對於一個矩陣\(A\)\(A\)可逆的充分必要條件是\(A\)

經過若干次初等行變換可以變成\(E\)\(E\)即單位矩陣),即存在一個矩陣\(P\)使得\(PA=E\),則\(P=A^{-1}\)

通過初等行變換使得\(A\)變為\(E\)並不困難,可以用高斯消元解決,先消成上三角矩陣,然後再消成對角矩陣

考慮怎麼求出\(P\),因為有\(PA=E,PE=P\),如果我們同時維護兩個矩陣\(A,B\),令\(B\)一開始時等於\(E\),在把\(A\)變為\(E\)的過程中對\(B\)也做相等的初等變換,那麼當\(A\)變為\(E\)時,\(B\)也就變為了\(P\)(因為做初等行變換等價於被對應的初等矩陣左乘)

如果在高斯消元的過程中發現無法將\(A\)

變為\(E\),輸出無解即可

//minamoto
#include<bits/stdc++.h>
using namespace std;
#define getc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
char buf[1<<21],*p1=buf,*p2=buf;
int read(){
    int res,f=1;char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int C=-1,Z=0;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
void write(int x){
    if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]=' ';
}
const int N=405,mod=1e9+7;
int n;
struct Matrix{
    int a[N][N];
    inline void clr(){memset(a,0,sizeof(a));}
    int* operator [](int x){return a[x];}
    void SWAP(int x,int y){for(int i=1;i<=n;++i)swap(a[x][i],a[y][i]);}
    //交換某兩行
    void MUL(int x,int k){for(int i=1;i<=n;++i)a[x][i]=(1ll*a[x][i]*k%mod+mod)%mod;}
    //將某一行的所有元素乘上k
    void MD(int x,int y,int k){for(int i=1;i<=n;++i)a[x][i]=((a[x][i]+(1ll*a[y][i]*k%mod))%mod+mod)%mod;}
    //將某一行的所有元素乘上k加到另一行去
    void print(){
        for(int i=1;i<=n;++i){
            for(int j=1;j<=n;++j)write(a[i][j]);
            sr[++C]='\n';
        }
    }
}A,B;
int ksm(int a,int b=(mod-2)){
    int res=1;
    for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)res=1ll*res*a%mod;
    return res;
}
int main(){
//  freopen("testdata.in","r",stdin);
    n=read();
    for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)A[i][j]=read();
    for(int i=1;i<=n;++i)B[i][i]=1;
    for(int i=1;i<=n;++i){
        //消成上三角矩陣 
        if(!A[i][i]){
            for(int j=i+1;j<=n;++j)if(A[j][i]){
                A.SWAP(i,j),B.SWAP(i,j);break;
            }
        }
        if(!A[i][i])return puts("No Solution"),0;
        //如果消著消著某一列沒有數了,說明無解 
        B.MUL(i,ksm(A[i][i])),A.MUL(i,ksm(A[i][i]));
        for(int j=i+1;j<=n;++j)
        B.MD(j,i,-A[j][i]),A.MD(j,i,-A[j][i]);
    }
    //消成對角矩陣 
    for(int i=n-1;i;--i)for(int j=i+1;j<=n;++j)
    B.MD(i,j,-A[i][j]),A.MD(i,j,-A[i][j]);
    B.print();return Ot(),0;
}