[学习笔记]矩阵求逆

以前只会高斯消元

高斯消元可以O(n^4)解决

但是可以做到O(n^3)

考虑要求

P*A=E

有

E*A=A

发现，高斯消元就是不断进行“交换两行，某一行乘上k，某一行乘上k加到另一行上去”

这三种变换，都可以通过左乘一个唯一的矩阵等价得到！（而且手玩发现，如果开始另一个矩阵是单位矩阵，就是单位矩阵做同样的变化乘过去就是了）

所以，开始让P=E

然后A不断高斯消元向E靠拢，P同时做完全相同的操作

最后A变成E，

那么P乘A，就是E了。

如果中途存在一个A[i][i]等于0，那么意味着无论如何不能削成单位矩阵。一定无解。（因为变换是可逆的，不可能一个A既能削成单位矩阵又不能）

否则由于可以构造出P，一定有解。

一种镜像映射的感觉

#include<bits/stdc++.h>
#define reg register int
#define il inline
#define fi first
#define se second
#define mk(a,b) make_pair(a,b)
#define numb (ch^'0')
using namespace std;
typedef long long ll;
template<class T>il void rd(T &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');}
template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');}
template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('
');}

namespace Miracle{
const int N=404;
const int mod=1e9+7;
int n;
void _swap(int &x,int &y){
    x^=y^=x^=y;
}
struct tr{
    int a[N][N];
    int* operator [](int x){return a[x];}
    void swap(int x,int y){
        for(reg i=1;i<=n;++i) _swap(a[x][i],a[y][i]);
    }
    void mul(int x,int k){
        for(reg i=1;i<=n;++i) a[x][i]=(ll)a[x][i]*k%mod;
    }
    void mov(int x,int y,int k){//x+=y*k
        for(reg i=1;i<=n;++i) a[x][i]=((ll)a[x][i]+(ll)a[y][i]*k%mod)%mod;
    }
}A,B;
int qm(int x,int y){
    int ret=1;
    while(y){
        if(y&1) ret=(ll)ret*x%mod;x=(ll)x*x%mod;y>>=1; 
    }
    return ret;
}
int main(){
    rd(n);
    for(reg i=1;i<=n;++i){
        for(reg j=1;j<=n;++j){
            rd(A[i][j]);
        }
    }
    for(reg i=1;i<=n;++i) B[i][i]=1;
    bool fl=true;
    for(reg i=1;i<=n;++i){
        if(A[i][i]==0){
            for(reg j=i+1;j<=n;++j){
                if(A[j][i]){
                    A.swap(i,j);B.swap(i,j);break;
                }
            }
        }
        if(!A[i][i]) {
            fl=false;break;
        }
        B.mul(i,qm(A[i][i],mod-2));A.mul(i,qm(A[i][i],mod-2));
        for(reg j=i+1;j<=n;++j){
//            if(A[j][i]){
                B.mov(j,i,mod-A[j][i]);A.mov(j,i,mod-A[j][i]);
//            }
        }
    }
//    for(reg i=1;i<=n;++i){
//        for(reg j=1;j<=n;++j){
//            cout<<A[i][j]<<" ";
//        }cout<<endl;
//    }cout<<endl;
    if(!fl) puts("No Solution");
    else{
        for(reg i=n;i>=2;--i){
            for(reg j=i-1;j>=1;--j){
                if(A[j][i]){
                    B.mov(j,i,mod-A[j][i]);A.mov(j,i,mod-A[j][i]);
                }
            }
//            for(reg k=1;k<=n;++k){
//                for(reg j=1;j<=n;++j){
//                    cout<<A[k][j]<<" ";
//                }cout<<endl;
//            }cout<<endl;
        }
        for(reg i=1;i<=n;++i){
            prt(B[i],1,n);
        }
    }
    
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2019/3/31 10:25:05
*/