BZOJ 3122 [SDOI2013]随机数生成器 (BSGS)
观察式子$x_{i+1}=a\cdot x_{i}+b (mod\;p)$
我们能用欧几里得算法求解$ax=c(mod\;b)$的线性同余方程问题
还能用$BSGS$算法求解$a^{x}=c(mod\;b)$的指数同余方程问题
而式子里的$b$这一项显得十分棘手
。
。
。
还记得等比数列吗,是数学老师教的
考虑构造一个"等比数列”方便我们递推
$(x_{i+1}+K)=a\cdot (x_{i}+K) (mod\;p)$
解得$K=\frac{b}{a-1}$
这样,原式化成了$A^{X}=B(mod\;p)$的形式
好像求出指数$x$就行了?
于是乎这道题就变成了一个裸的$BSGS$问题?
$a==0$的情况特判
$a==1?$原式变成了等差数列,用$exgcd$求解线性同余方程即可
1 #include <cmath> 2 #include <cstdio> 3 #include <cstring> 4 #include <algorithm> 5 #define N1 400010 6 #define M1 40010 7 #define ll long long 8 #define dd double 9 #define cint const int 10 #define cll const long long 11 #define inf 23333333333333333ll 12 using namespace std; 13 14 int gcd(int a,int b){ if(!b) return a; return gcd(b,a%b); } 15 void exgcd(ll a,ll b,ll &x,ll &y) 16 { 17 if(!b){ x=1; y=0; return; } 18 exgcd(b,a%b,x,y); ll t=x; x=y; y=t-a/b*y; 19 } 20 int P,A,B,X1,T,K; 21 namespace S0{ 22 int solve() 23 { 24 int ans,g; ll inv,invy; T=(T-X1+P)%P; 25 g=gcd(B,P); if(T%g!=0) return -1; 26 T/=g; B/=g; P/=g; 27 exgcd(B,P,inv,invy); inv=(inv%P+P)%P; 28 ans=inv*T%P; P*=g; ans=1ll*ans*g%P; 29 return ans+1; 30 } 31 }; 32 namespace S1{ 33 #define maxn 400000 34 struct Hsh{ 35 int head[N1],to[M1],nxt[M1],val[M1],cte; 36 void ins(int x,int w) 37 { 38 int u=x%maxn,j,v; 39 for(j=head[u];j;j=nxt[j]) 40 { 41 v=to[j]; 42 if(v==x) return; 43 } 44 cte++; to[cte]=x; nxt[cte]=head[u]; 45 head[u]=cte; val[cte]=w; 46 } 47 int find(int x) 48 { 49 int u=x%maxn,j,v; 50 for(j=head[u];j;j=nxt[j]) 51 { 52 v=to[j]; 53 if(v==x) return val[j]; 54 } 55 return -1; 56 } 57 }h; 58 ll qpow(ll x,ll y,cll &P) 59 { 60 ll ans=1; 61 while(y){ 62 if(y&1) ans=ans*x%P; 63 x=x*x%P; y>>=1; 64 }return ans; 65 } 66 #undef maxn 67 int solve(int p) 68 { 69 int i,j,sq;ll K,pw,now,ans,inv,invy,Y,tmp; cll P=p; 70 exgcd(A-1,P,inv,invy); inv=(inv%P+P)%P; K=inv*B%P; 71 T=(T+K)%P; X1=(X1+K)%P; 72 exgcd(X1,P,inv,invy); inv=(inv%P+P)%P; T=T*inv%P; 73 sq=sqrt(P); 74 for(pw=qpow(A,sq,P),now=1,i=1;(i-1)*sq<P;i++) 75 { 76 now=1ll*now*pw%P; 77 h.ins(now,i); 78 } 79 for(ans=inf,now=T,i=0;i<sq;i++) 80 { 81 tmp=h.find(now); 82 if(tmp!=-1) ans=min(ans,tmp*sq-i); 83 now=now*A%P; 84 } 85 memset(&h,0,sizeof(h)); 86 if(ans==inf) return -1; 87 return ans+1; 88 } 89 }; 90 int Case; 91 92 int main() 93 { 94 scanf("%d",&Case); 95 while(Case--) 96 { 97 scanf("%d%d%d%d%d",&P,&A,&B,&X1,&T); 98 if(X1==T){ puts("1"); continue; } 99 if(A==0) 100 { 101 if(T==B) puts("2"); 102 else puts("-1"); 103 continue; 104 } 105 if(A==1) printf("%d\n",S0::solve()); 106 else printf("%d\n",S1::solve(P)); 107 } 108 return 0; 109 }