顺便学了一下矩阵求逆(其实就是一个高斯消元),不过后来没有用。
采用BSGS的思想,实际上令m=sqrt(p),那么就是求最小的i*m-j使得:
(A^m)^i=A^jB,其中A^j一共m个,可以预处理。
因此关键是快速判断两个矩阵是否相等。可以将矩阵看成一个大的p进制数然后将其取模后得到的值设为其hash值,判断为O(1);或者用另外一种方法,即将矩阵去乘一个(n*1)的矩阵,得到一个(n*1)的矩阵,然后就可以O(N)判断了。后面一种实现更简单,而且更快(用小号刷到rk2辣)。
时间复杂度O(P^0.5N^2+N^3+PN)。
AC代码如下:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#define N 75
using namespace std;
int n,mod,num[N];
struct matrix{
int p[N][N];
void init(){
int i,j;
for (i=1; i<=n; i++)
for (j=1; j<=n; j++) scanf("%d",&p[i][j]);
}
}a,b,c,hash[N<<1];
matrix tms(matrix x,matrix y){
matrix z; int i,j,k;
for (i=1; i<=n; i++)
for (j=1; j<=n; j++){
z.p[i][j]=0;
for (k=1; k<=n; k++) z.p[i][j]=(z.p[i][j]+x.p[i][k]*y.p[k][j])%mod;
}
return z;
}
matrix mul(matrix x,matrix y){
matrix z; int i,j;
for (i=1; i<=n; i++)
for (j=1; j<=n; j++){
z.p[i][j]=0;
z.p[i][1]=(z.p[i][1]+x.p[i][j]*y.p[j][1])%mod;
}
return z;
}
matrix ksm(matrix x,int y){
matrix z=x;
for (; y; y>>=1,x=tms(x,x)) if (y&1) z=tms(z,x);
return z;
}
bool ok(matrix x,matrix y){
int i; for (i=1; i<=n; i++) if (x.p[i][1]!=y.p[i][1]) return 0;
return 1;
}
int main(){
scanf("%d%d",&n,&mod); int m=(int)sqrt(mod)+1,i,j;
srand(20160307); a.init(); b.init();
for (i=1; i<=n; i++) c.p[i][1]=num[i]=rand()%mod;
for (i=0; i<m; i++){
hash[i]=mul(b,c); c=mul(a,c);
}
a=ksm(a,m-1);
for (i=1; i<=n; i++) c.p[i][1]=num[i];
for (i=1; i<=m; i++){
c=mul(a,c);
for (j=m-1; j>=0; j--)
if (ok(c,hash[j])){ printf("%d\n",i*m-j); return 0; }
}
return 0;
}
by lych
2016.3.7