BZOJ 4128

矩阵BSGS，同样处理即可。

有一个优化可以更快，判断矩阵相等可以$O(n)$。列矩阵的使用。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=71;
int n,m,P;
struct matrix{
    int a[N][N];
    int t;
    matrix(){}
    int* operator [] (int x){return a[x];}
    matrix operator * (matrix x){
        matrix tmp ;
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++){
                tmp[i][j]=0;
                for (int k=0;k<n;k++)
                    tmp[i][j]=(tmp[i][j]+a[i][k]*x[k][j])%P;
            }
        return tmp;
    }
    void operator *= (matrix x){*this = *this * x;}
    matrix operator ^ (ll x){
        matrix ret;
        for (int i=0;i<n;i++)ret[i][i]=1;
        matrix tmp = *this;
        for (;x;x>>=1,tmp*=tmp){if(x&1)ret *=tmp;}
        return ret;
    }
    void print(){
        for (int i=0;i<n;i++){
            for (int j=0;j<n;j++)
                printf("%d ",a[i][j]);
            puts("");
        }
    }
    bool operator == (const matrix x) const{
    	for(int i=0;i<n;i++)
    		for(int j=0;j<n;j++)
    			if(a[i][j]!=x.a[i][j]) return false;
    	return true;
    }
};
struct Hash{
	#define Maxn 10007
	vector<matrix> mp[Maxn];
	void init(){
		for(int i=0;i<Maxn;i++)
			mp[i].clear();
	}
	int getkey(matrix t){
		ll ret=0;
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++)
				ret=(ret*817+t[i][j])%Maxn;
		return ret;
	}
	void insert(matrix tmp,ll i){
		int key=getkey(tmp);
		tmp.t=i;
		mp[key].push_back(tmp);
	}
	int Search(matrix t){
		int key=getkey(t);
		int mx=-1;
		// for(auto E:mp[key]){
		// 	if(E==t){
		// 		mx=max(mx,E.t);
		// 	}
		// } 
		for(int i=0;i<(int)mp[key].size();i++){
			if(mp[key][i]==t){
				mx=max(mx,mp[key][i].t);
			}
		}
		return mx;
	}
}T;

inline ll BSGS(matrix a,matrix b,int p){//a^x%p=b
	ll l=1,r=p;
	ll m;
	while(l<=r){
		ll mid=(l+r)>>1;
		if(mid*mid>=p)
			r=mid-1,m=mid;
		else l=mid+1;
	}
	matrix tmp=b;
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			tmp[i][j]%=p;
	T.init();
	for(int i=0;i<=m;i++){
		T.insert(tmp,i);
		tmp=tmp*a;
	}
	tmp=a^m;
	matrix ret;
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			if(i==j) ret[i][j]=1;
			else ret[i][j]=0;
	for(int i=1;i<=m;i++){
		ret=ret*tmp;
		int id=T.Search(ret);
		if(id!=-1){
			return m*i-id;
		}
	}
	return -1;
}

int main(){
	matrix a,b;
	scanf("%d%d",&n,&P);
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			scanf("%d",&a[i][j]);
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			scanf("%d",&b[i][j]);
	cout<<BSGS(a,b,P);
	return 0;
}