bzoj 3122 随机数生成器 （BSGS）

3122: [Sdoi2013]随机数生成器

Time Limit: 10 Sec   Memory Limit: 256 MB
Submit: 1049   Solved: 430
[ Submit][ Status][ Discuss]

Description

Input

输入含有多组数据，第一行一个正整数T，表示这个测试点内的数据组数。  
 
接下来T行，每行有五个整数p，a，b，X1，t，表示一组数据。保证X1和t都是合法的页码。 

注意：P一定为质数

Output

共T行，每行一个整数表示他最早读到第t页是哪一天。如果他永远不会读到第t页，输出-1。 

Sample Input

3
7 1 1 3 3
7 2 2 2 0
7 2 2 2 1

Sample Output

1
3
-1

HINT

0<=a<=P-1,0<=b<=P-1,2<=P<=10^9

题解：这是一道不错的数论题，困扰了我好久，最终AC了。

如果问我做这道题最大的感慨是什么的话

1.数论太神啦

2.没事就取模，内存够果断开long long （因为炸飞调了好久）

这道题要分情况讨论，所以有很多看似很简单，但很不容易想全的特判

1.x1==t ,直接输出1

2.a==0 (1)b==t 输出2

             (2)否则输出-1

3.a==1

  式子就可以化简成

 xn=x1+(ans-1)b mod p

整理得   xn-x1=(ans-1)b mod p

    令  z=xn-x1, x=ans-1

      bx=z mod p

即  bx+mp=z  其中b,p,z 为已知，就可以用扩张欧几里德先求bx+mp=1 的解再乘z

当然也可以 ans-1 =(xn-x1)*inv(b)  ,inv(b)表示b 的逆元，在模意义下，除以b等于乘b的逆元。

根据费马小定理  a,p 互质，p为质数  a^p-1=1 mod p  ,ax=1 mod p a的逆元x=a^p-2;

4. 将原式化简   

x[i+1]=ax[i]+b(mod p)

x[i+1]+b/(a-1)=a(x[i]+b/(a-1))] (mod p)  把 x[i+1]=ax[i]+b mod p代入得

x[n]+b/(a-1)=a^(n-1)(x1+b/(a-1)) (mod p)  x[n]=t

(t+b*inv(a-1))*inv(x1+b*inv(a-1))=a^(n-1)(mod p)  inv

整理可得  a^x=z mod p  x=n-1, z=(t+b*inv(a-1))*inv(x1+b*inv(a-1))

这样就可以用bsgs求解

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<map>
using namespace std;
long long  a,b,p,t,x,n,z,ax,ay;
map<long long,long long> mp;
void exgcd(long long  a,long long b)
{
  if (b==0)
   {
   	 ax=1; ay=0; return;
   }
  else
   {
   	  exgcd(b,a%b);
   	  long long t=ay;
   	  ay=ax-(a/b)*ay;
   	  ax=t;
   }
}
long long gcd(long long  x,long long y)
{
	long long r=0;
	while(y!=0)
	 {
	   r=x%y;
	   x=y;
	   y=r;
	 }
	return x;
}
long long  quickpow(long long  num,long long x)
{
	long long k=num%p; long long ans=1;
	while (x>0)
	 {
	 	if (x&1)
	 	 ans=(ans*k)%p;
	 	x=x>>1;
	 	k=(k*k)%p;
	 }
	return ans;
}
int main()
{
  freopen("input.in","r",stdin);
  freopen("output.out","w",stdout);
  scanf("%d",&n);
  for (int l=1;l<=n;l++)
   {
   	 bool pd=false;
   	 scanf("%I64d%I64d%I64d%I64d%I64d",&p,&a,&b,&x,&t);
   	 if(x==t)
   	  {
   	  	printf("1\n");
   	  	continue;
   	  }
   	 if (a==0)
   	  if (b==t)
   	   {
   	   	 printf("2\n");
   	   	 continue;
   	   }
   	 else
   	  {
   	  	printf("-1\n");
   	  	continue;
   	  }
   	 if (a==1)
   	  {
   	  	z=((t-x)%p+p)%p;
   	  	if (z%(gcd(b,p))!=0)
   	  	 {
   	  	 	printf("-1\n");
   	  	 	continue;
   	  	 }
   	  	exgcd(b,p);
   	  	ax=(ax%p*z%p)%p;
   	  	printf("%d\n",(ax%p+p)%p+1);
   	  	continue;
   	  }
   	 long long inva=quickpow(a-1,p-2)%p;
   	 long long invb=quickpow(x+b*inva,p-2)%p;
   	 long long nk=(long long)((t+b*inva)%p*invb%p)%p;
   	 mp.clear();
   	 long long ans=nk;
   	 long long m=ceil(sqrt(p));
   	 mp[ans%p]=0;
   	 for (long long i=1;i<=m;i++)
   	  {
   	  	ans=(long long)(ans%p*a%p)%p;
		mp[ans]=i;
   	  }
   	 long long an=quickpow(a,m); ans=1;
   	 for(long long i=1;i<=m;i++)
   	  {
   	  	ans=(long long)(ans%p*an%p)%p;
   	  	if (mp[ans])
   	  	 {
   	  	 	printf("%d\n",((i*m-mp[ans])%p+p)%p+1);
   	  	 	pd=true;
   	  	 	break;
   	  	 }
   	  }
   	 if (!pd) printf("-1\n");
   }
}