Public model for matrix

   以下是可以加减乘除(就是乘逆矩阵啦)以及求若干次幂、行列式和逆的矩阵模板。

欢迎大家指正其中可能存在的错误(只验证了求逆的正确性)。

   顺便提一下这种复杂度低于定义式求逆的方法，来自于我的高等代数书，思想就是对分块矩阵(A E)进行行变换从而得到(E A^-1)，复杂度与消元一样，都是 O(N^3)的。

 

const int N=405,ha=1e9+7;

inline int add(int x,int y){ x+=y; return x>=ha?x-ha:x;}
inline void ADD(int &x,int y){ x+=y; if(x>=ha) x-=ha;}

inline int ksm(int x,int y){
	int an=1;
	for(;y;y>>=1,x=x*(ll)x%ha) if(y&1) an=an*(ll)x%ha;
	return an;
}
inline int Get_inv(int x){ return ksm(x,ha-2);}

struct matrix{
	int a[N][N],n;
	
	inline void clear(){ memset(a,0,sizeof(a));}
	inline void Base(){ clear(); for(int i=1;i<=n;i++) a[i][i]=1;}
	
	inline void input(){
		scanf("%d",&n);
		for(int i=1;i<=n;i++)
		    for(int j=1;j<=n;j++) scanf("%d",a[i]+j);
	}
	inline void output(){
		for(int i=1;i<=n;i++){
			for(int j=1;j<=n;j++) printf("%d ",a[i][j]);
			puts("");
		}
	}
	
	inline void swap(int x,int y){
		for(int i=1;i<=n;i++) swap(a[x][i],a[y][i]);
	}
	inline void add(int from,int tmp,int to){
		for(int i=1;i<=n;i++) ADD(a[to][i],a[from][i]*(ll)tmp%ha);
	}
	inline void mul(int x,int ml){
		for(int i=1;i<=n;i++) a[x][i]=a[x][i]*(ll)ml%ha;
	}
	
	inline int Determinant(){
		int an=1;
		matrix b=*this;
		
		for(int i=1,o,inv;i<=n;i++){
			for(o=i;o<=n;o++) if(b.a[o][i]) break;
			if(o>n) return 0;
			if(o>i){ an=ha-an; b.swap(o,i);}
			
			inv=Get_inv(b.a[i][i]);
			for(int j=i+1,tmp;j<=n;j++) if(b.a[j][i])
				tmp=ha-inv*(ll)b.a[j][i]%ha,b.add(i,tmp,j);
			
			an=an*(ll)b.a[i][i]%ha;
		}
		
		return an;
	}
	
	inline matrix ni(){
		matrix b; b.n=n,b.Base();
		
		if(!Determinant()) b.clear();
		else{
			
			//先化成对角线元素都是1的上三角矩阵 
			
    		for(int i=1,o,inv;i<=n;i++){
    			for(o=i;o<=n;o++) if(a[o][i]) break;
    			if(o>i) swap(o,i),b.swap(o,i);
    			
    			inv=Get_inv(a[i][i]),mul(i,inv),b.mul(i,inv);
    			
    			for(int j=i+1,tmp;j<=n;j++) if(a[j][i])
    				tmp=ha-a[j][i],add(i,tmp,j),b.add(i,tmp,j);
    		}
    		
    		//然后再把A消成单位矩阵 
    		
    		for(int i=n-1;i;i--)
    		    for(int j=i+1;j<=n;j++) if(a[i][j]) b.add(j,ha-a[i][j],i),add(j,ha-a[i][j],i);
		}
		
		return b;
	}
	
	matrix operator +(const matrix &u){
		for(int i=1;i<=n;i++)
		    for(int j=1;j<=n;j++) ADD(a[i][j],u.a[i][j]);
		return *this;
	}
	
	matrix operator -(const matrix &u){
		for(int i=1;i<=n;i++)
		    for(int j=1;j<=n;j++) ADD(a[i][j],ha-u.a[i][j]);
		return *this;
	}
	
	matrix operator *(const matrix &u)const{
		matrix r; r.clear(),r.n=n;
		for(int k=1;k<=n;k++)
		    for(int i=1;i<=n;i++)
		        for(int j=1;j<=n;j++) ADD(r.a[i][j],a[i][k]*(ll)u.a[k][j]%ha);
		return r;
	}
	
	matrix operator /(matrix &u)const{
		return *this*u.ni();
	}
	
	matrix operator ^(int &u){
		matrix b,c=*this; b.n=n,b.Base();
		for(;u;u>>=1,c=c*c) if(u&1) b=b*c;
		return b;
	}
}JHY;