【BZOJ3122】随机数生成器（SDOI2013）-BSGS+矩阵乘法

测试地址：随机数生成器
做法：本题需要用到BSGS+矩阵乘法。
首先看到线性递推式，就想到用矩阵乘法来做，那么题目要求的就是关于x$x$的方程：

(a0b1)x(x11)=(t1)$${\left(\begin{array}{cc}a& b\\ 0& 1\end{array}\right)}^x\left(\begin{array}{c}{x}_1\\ 1\end{array}\right)=\left(\begin{array}{c}t\\ 1\end{array}\right)$$

的最小非负整数解（以上等式表示在模p$p$意义下等），那么x+1$x+1$就是我们所求的答案。
然而注意到求这个方程比较困难，因为x$x$在指数的位置。我们联想到，如果式子中不是矩阵而是整数的话，我们可以用BSGS在O(p–√)$O(\sqrt{p})$的时间内求出解，但是这里换成了矩阵我们要怎么做呢？
其实非常相似，令上式中以x$x$为指数的矩阵为A$A$，我们构造一个矩阵M$M$，使得MA=E$MA=E$，其中
E$E$为单位矩阵，显然：

M=(1a0−ba1)$$M=\left(\begin{array}{cc}\frac{1}{a}& -\frac{b}{a}\\ 0& 1\end{array}\right)$$

那么我们仿照标准的BSGS算法，令x=im+j$x=im+j$，那么上式就会变成：

Aim(x11)=Mj(t1)$$A^{im}\left(\begin{array}{c}{x}_1\\ 1\end{array}\right)=M^j\left(\begin{array}{c}t\\ 1\end{array}\right)$$

于是我们先对右边求出所有可能的取值，存在哈希表中，然后对于左边可能的取值，在哈希表中查询有没有相同的元素存在，这其实和标准的BSGS算法完全一样。那么我们只需要令m$m$为p–√$\sqrt{p}$级别的一个数，这个算法的复杂度就是O(p–√)$O(\sqrt{p})$的了。
还有一点要注意，当a=0$a=0$时，M$M$不存在，所以特判一下即可。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll hashmod=1000003;
int T,q[1000010],top=0,blocksize;
ll p,a,b,x1,t,hashlist[1000010],val[1000010]={0};
struct matrix
{
    ll mat[2][2];
}E,A,S,M;

void hashinsert(ll x,ll v)
{
    int now=x%hashmod;
    while (hashlist[now]!=-1&&hashlist[now]!=x) now=(now+1)%hashmod;
    if (hashlist[now]==-1) hashlist[now]=x,val[now]=v,q[++top]=now;
}

ll hashfind(ll x)
{
    int now=x%hashmod;
    while (hashlist[now]!=-1&&hashlist[now]!=x) now=(now+1)%hashmod;
    if (hashlist[now]==x) return val[now];
    else return -1;
}

ll power(ll a,ll b)
{
    ll s=1,ss=a;
    while(b)
    {
        if (b&1) s=(s*ss)%p;
        ss=(ss*ss)%p;b>>=1;
    }
    return s;
}

matrix mult(matrix A,matrix B)
{
    matrix S;
    memset(S.mat,0,sizeof(S.mat));
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
            for(int k=0;k<2;k++)
                S.mat[i][j]=(S.mat[i][j]+A.mat[i][k]*B.mat[k][j])%p;
    return S;
}

void init()
{
    for(int i=1;i<=top;i++) hashlist[q[i]]=-1;
    scanf("%lld%lld%lld%lld%lld",&p,&a,&b,&x1,&t);
    blocksize=(int)(sqrt(p)+1);
    A.mat[0][0]=a,A.mat[0][1]=b,A.mat[1][0]=0,A.mat[1][1]=1;
    S=E;
    for(int i=1;i<=blocksize;i++) S=mult(S,A);
    M.mat[0][0]=power(a,p-2),M.mat[0][1]=((-b*power(a,p-2))%p+p)%p;
    M.mat[1][0]=0,M.mat[1][1]=1;
}

void BSGS()
{
    matrix N=E;
    top=0;
    for(int i=0;i<blocksize;i++,N=mult(N,M))
        hashinsert((N.mat[0][0]*t+N.mat[0][1])%p,i);
    N=E;
    for(int i=0;i<=blocksize;i++,N=mult(N,S))
    {
        ll s=hashfind((N.mat[0][0]*x1+N.mat[0][1])%p);
        if (s!=-1)
        {
            printf("%lld
",i*blocksize+s+1);
            return;
        }
    }
    printf("-1
");
}

int main()
{
    scanf("%d",&T);
    E.mat[0][0]=E.mat[1][1]=1;
    E.mat[0][1]=E.mat[1][0]=0;
    for(int i=0;i<=hashmod;i++) hashlist[i]=-1;
    while(T--)
    {
        init();
        if (a) BSGS();
        else
        {
            if (t==x1) printf("1
");
            else if (t==b) printf("2
");
                 else printf("-1
");
        }
    }

    return 0;
}