BZOJ4128: Matrix BSGS+矩阵求逆+线性筛

题解：

首先考虑如何对一个可逆矩阵求逆。

据称有一个理论就是任意一个可逆矩阵都是单位矩阵通过基本行变换以及基本列变换形成的。

毕竟我们这样搞肯定搞不出不可逆矩阵的吧。

所以我们对矩阵进行消元变成单位矩阵，假设这个过程相当于乘上了一个矩阵$B$。

也就是说：$A\times{B}=I$。

事实上说明了$B=A^{-1}$。

所以我们让一个单位矩阵也和矩阵$A$进行相同的基本变换，最终就能得到$A$的逆矩阵了。

再加之BSGS那套理论就行了。

这样求逆元和乘法的复杂度都是$O(n^3)$。

这样复杂度就是$O(n^3\sqrt{p})$吧。

有一些细节，比如利用hash表，以及预处理逆元。

这里利用到了线性筛逆元大法好！

比如我们想求$a$的逆元$a^{-1}$，并保证在此之前的逆元都已经被求出。

不妨令$p=r\times{a}+q$,其中$0<q<a$（因为必然不整除），那么也就是说：

$r\times{a}+q\equiv{0}\pmod{p}$

两边同时乘以$q^{-1}$得到：

$r\times{a}\times{q^{-1}}+1\equiv{0}\pmod{p}$

移项得到：

$-r\times{q^{-1}}\times{a}\equiv{1}\pmod{p}$

于是得到：

$a^{-1}\equiv{-r\times{q^{-1}}}\pmod{p}$

而显然$q^{-1}$我们之前已经求出，于是这样就能做到$O(n)$求出所有逆元了。

代码：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
using namespace std;

typedef unsigned long long llu;
static const int seed=131;


int p,inv[20010];

#define N 71
int n;
struct Matrix{
	int A[N][N];
	Matrix(){
		memset(A,0,sizeof A);
	}
	inline void read(){
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				scanf("%d",&A[i][j]);
	}
	inline llu hash(){
		llu s=1,re=0;
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				re^=(s*A[i][j]),s*=seed;
		return re;
	}
	inline Matrix get_inv(){
		static int B[N][N<<1];
		int i,j,k,t;
		memset(B,0,sizeof B);
		for(i=1;i<=n;++i)
			for(j=1;j<=n;++j)
				B[i][j]=A[i][j];
		for(i=1;i<=n;++i)
			B[i][i+n]=1;
		for(i=1;i<=n;++i){
			for(j=i+1;j<=2*n;++j)
				(B[i][j]*=inv[B[i][i]])%=p;
			B[i][i]=1;
			for(j=i+1;j<=n;++j)
				for(t=B[j][i],k=i;k<=2*n;++k)
					B[j][k]=((B[j][k]-t*B[i][k])%p+p)%p;
		}
		for(i=n;i>=1;--i)
			for(j=i+1;j<=n;++j)
				for(t=B[i][j],k=j;k<=2*n;++k)
					B[i][k]=((B[i][k]-t*B[j][k])%p+p)%p;
		Matrix re;
		for(i=1;i<=n;++i)
			for(j=1;j<=n;++j)
				re.A[i][j]=B[i][j+n];
		return re;
	}
	inline Matrix operator*(const Matrix&B)const{
		Matrix c;
		int i,j,k;
		for(i=1;i<=n;++i)
			for(j=1;j<=n;++j)
				for(k=1;k<=n;++k)
					c.A[i][j]=(c.A[i][j]+A[i][k]*B.A[k][j])%p;
		return c;
	}
	inline void operator*=(const Matrix B){
		Matrix c=*this;
		*this=c*B;
	}
};

struct Hashset{
	static const int mod=131313;
	int head[mod],next[20010],v[20010],ind;
	llu f[20010];
	inline void reset(){
		ind=0;
		memset(head,-1,sizeof head);
	}
	inline int operator[](const llu x){
		int ins=x%mod;
		for(int j=head[ins];j!=-1;j=next[j])
			if(f[j]==x)
				return v[j];
		return -1;
	}
	inline void insert(llu x,int _v){
		int ins=x%mod;
		for(int j=head[ins];j!=-1;j=next[j])
			if(f[j]==x){
				v[j]=min(v[j],_v);
				return;
			}
		f[ind]=x;
		v[ind]=_v;
		next[ind]=head[ins];
		head[ins]=ind++;
	}
}M;

int main(){
#ifndef ONLINE_JUDGE
	freopen("tt.in","r",stdin);
#endif
	int i;
	scanf("%d%d",&n,&p);
	for(inv[1]=1,i=2;i<p;++i)
		inv[i]=p-(p/i)*inv[p%i]%p;
	Matrix A,inv_A,B,C;
	A.read();
	B.read();
	int m=ceil(sqrt(p));
	for(i=1;i<=n;++i)
		C.A[i][i]=1;

	M.reset();
	for(i=0;i<m;++i,C*=A)
		M.insert(C.hash(),i);
	
	C=C.get_inv();
	int re;
	for(i=0;i*m<p;++i,B*=C){
		re=M[B.hash()];
		if(re!=-1){
			cout<<i*m+re<<endl;
			break;
		}
	}
	return 0;
}