bzoj 4128 Matrix

http://www.elijahqi.win/archives/3248
Description

给定矩阵A,B和模数p,求最小的x满足

A^x = B (mod p)
Input

第一行两个整数n和p，表示矩阵的阶和模数,接下来一个n * n的矩阵A.接下来一个n * n的矩阵B
Output

输出一个正整数，表示最小的可能的x，数据保证在p内有解
Sample Input

2 7
1 1
1 0
5 3
3 2
Sample Output

4
HINT

对于100%的数据，n <= 70,p <=19997，p为质数,0<= A_{ij},B_{ij}< p
保证A有逆
为了避免矩阵求逆 所以使用另外一种形式的bsgs 注意这个算需要算到ceil(sqrt(mod)) 并且储存在hash表里的东西需要从大到小储存
即i*m-j=x 那么将这个j挪动到右边就是正的次幂了
然后针对每个矩阵做hash再把hash值存到hash表里即可
证明只需要做到ceil即可
a^(k%(p-1))=a^k%p
k%(p-1)=k-m*(p-1)
原式看成a^k/(a^(p-1)^m) 因为费马小定理 a^(p-1)=1所以 变成了a^k=a^k

#include<cmath>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
#define ll unsigned long long
using namespace std;
inline char gc(){
    static char now[1<<16],*S,*T;
    if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(!isdigit(ch)) {if (ch=='-') f=-1;ch=gc();}
    while(isdigit(ch)) x=x*10+ch-'0',ch=gc();
    return x*f;
}
const int g1=100003;
const int g2=1000003;
const int N=80;
ll p1[N],p2[N];int mod,n;
struct matrix{int f[N][N];}a,b,c,bas;
inline matrix multiply(const matrix &a,const matrix &b){
    matrix c;memset(c.f,0,sizeof(c.f));
    for (int i=1;i<=n;++i)
        for (int k=1;k<=n;++k)
            for (int j=1;j<=n;++j)
                (c.f[i][j]+=a.f[i][k]*b.f[k][j])%=mod;
    return c;
}
const int md=99989;
struct node{
    ll x;int j,next;
    node(){x=next=j=0;}
    node(ll &_x,int &_j,int &_next){x=_x;j=_j;next=_next;}
}data[md];
int h[md],num;
inline void insert1(ll x,int id){
    static int tmp;tmp=x%md;
    for (int i=h[tmp];i;i=data[i].next) if (data[i].x==x) {data[i].j=id;return;}
    data[++num]=node(x,id,h[tmp]);h[tmp]=num;
}
inline int query(ll x){
    static int tmp;tmp=x%md;
    for (int i=h[tmp];i;i=data[i].next) if (data[i].x==x) return data[i].j;
    return -1;
}
inline ll make_hs(const matrix &a){
    static ll hs;hs=0;
    for (int i=1;i<=n;++i)
        for (int j=1;j<=n;++j) hs+=p1[i-1]*p2[j-1]*a.f[i][j];
    return hs;
}
inline void print(const matrix &a){
    for (int i=1;i<=n;++i){
        for (int j=1;j<=n;++j) printf("%d ",a.f[i][j]);puts("");
    }
}
int main(){
    freopen("bzoj4128.in","r",stdin);
    n=read();mod=read();int m=ceil(sqrt(mod));p1[0]=1;p2[0]=1;
    for (int i=1;i<=n;++i) p1[i]=p1[i-1]*g1;
    for (int i=1;i<=n;++i) p2[i]=p2[i-1]*g2;
    for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) a.f[i][j]=read();bas=a;
    for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) b.f[i][j]=read();c=b;
    for (int i=1;i<=m;++i) c=multiply(c,a),insert1(make_hs(c),i);//print(c);
    memset(c.f,0,sizeof(c.f));for (int i=1;i<=n;++i) c.f[i][i]=1;//puts("---");
    for (int t=m;t;bas=multiply(bas,bas),t>>=1) if (t&1) c=multiply(c,bas);bas=c;
    for (int i=1;i<=m;++i){
        int j=query(make_hs(c));//print(c);
        if (j!=-1) {printf("%d\n",i*m-j);return 0;}c=multiply(c,bas);
    }
    return 0;
}