先備:1. 指數原根 2. Baby step-Giant step
這題真不是普通的好玩XDD
總之,照題意,要解x:
先抓p的原根g
然後對兩邊取index:
然後解Linear方程。
/* LANG:C++ ** http://web2.ck.tp.edu.tw/~step5/probdisp.php?pid=0035 ** NUMBER THEORY ** 2013-11-8 */ #include<cstring> #include<stdio.h> #include<stdlib.h> #include<algorithm> #include<math.h> #include<assert.h> #include<vector> #define PMAX 100000 #define Mid(a,b) ((a+b)>>1) #define LL long long int #define PS printf("PS[%d]\n",step++); using namespace std; int step=0; ///Set up Prime Table vector<LL> PRIME; bool seive[PMAX+1]; void BuildPrimeTable() { memset(seive,true,sizeof(seive)); seive[1]=false; PRIME.clear(); for(LL lx=2;lx<=PMAX;lx++) { if(seive[lx]) { PRIME.push_back(lx); for(LL ly=2;ly*lx<=PMAX;ly++) seive[lx*ly]=false; } } } ///Get the primitive root of P LL P; vector<LL> Drc; LL _Pow(LL a,LL n,LL _PP) { if(n==0) return 1; if(n==1) return a%_PP; LL k=_Pow(a%_PP,n/2,_PP); if(n%2) return (((k*k)%_PP)*a)%_PP; else return (k*k)%_PP; } LL _Pow(LL a,LL n) { return _Pow(a,n,P); } //536877277 void DrcPrime() { LL prc=P-1; Drc.clear(); //printf("a1"); for(int lx=0;(PRIME.size()>lx)&&(PRIME[lx]<prc);lx++) { if(prc%PRIME[lx]==0) { Drc.push_back(PRIME[lx]); while(prc%PRIME[lx]==0) prc/=PRIME[lx]; //printf("%I64d\n",PRIME[lx]); } } if(prc>1) Drc.push_back(prc); //printf("a1"); } LL GetPrimitiveRoot() { for(LL lx=2;lx<=P;lx++) { bool isRoot=true; for(LL ly=0;ly<Drc.size();ly++) if(_Pow(lx,(P-1)/Drc[ly])==1){isRoot=false;break;} if(isRoot) return lx; } } ///Babystep prt^a=r class Data { public: LL index; LL cont; Data(){} void set(const LL _i,const LL _c){index=_i,cont=_c;} }; bool operator<(Data _aa,Data _bb){return (_aa.cont<_bb.cont);} Data* brysch(Data *I,LL c,LL m) { //prLLf("SCH:%d\n",c); LL start=0,end=m-1; while((end-start)>1) { //prLLf("END=%d START=%d\n",end,start); if(I[start].cont==c) return &I[start]; else if(I[end].cont==c) return &I[end]; else if(I[Mid(start,end)].cont==c) return &I[Mid(start,end)]; if(c>I[Mid(start,end)].cont) start=Mid(start,end); else end=Mid(start,end); } if(I[start].cont==c) return &I[start]; else if(I[end].cont==c) return &I[end]; else return NULL; } Data sqrtb[PMAX]; LL GetInd(LL prt,LL r) { r=(r%P+P)%P;//prLLf("GetInd(%d,%d)\n",prt,r); LL m=ceil(sqrt(P)); for(LL lx=0,prc=1;lx<m;lx++,prc=(prc*prt)%P) sqrtb[lx].set(lx,prc); sort(sqrtb,sqrtb+m); /*for(LL lx=0;lx<m;lx++) prLLf("%d:%d\n",sqrtb[lx].index,sqrtb[lx].cont);*/ LL base=_Pow(prt,m*(P-2)); //printf("base=%d\n",base); for(LL lx=0,prc=1;lx<m;lx++,prc=(prc*base)%P) { //printf("prc=%d\n",prc); Data *RES=brysch(sqrtb,(r*prc)%P,m); if(RES==NULL) continue; return (*RES).index+m*lx; } return -1; } ///ax=b mod P-1 ,return x LL phi(LL m) { LL prc=P; LL ret=1; for(LL lx=0,prc=m;(lx<PRIME.size())&&(PRIME[lx]<=prc);lx++) { if(prc%PRIME[lx]==0) { LL nn=1;prc/=PRIME[lx]; while(prc%PRIME[lx]==0) prc/=PRIME[lx],nn*=PRIME[lx]; ret*=nn*(PRIME[lx]-1); } } return ret; } LL Linear(LL a,LL b) { LL MM=P-1; LL gg=__gcd(a,MM); //printf("gg=%I64d\n",gg); //assert(gg); MM/=gg; a/=gg; b/=gg; //printf("{%I64d}}\n",phi(MM)); return (b*_Pow(a,phi(MM)-1,MM))%MM; } int main() { BuildPrimeTable(); LL X0,X1,K,N; LL proot; long long int Kqq; //scanf("%d",&P); scanf("%I64d %I64d %I64d %I64d %I64d",&X0,&X1,&Kqq,&N,&P); //PS if(Kqq>=0) K=Kqq; else K=P+Kqq; DrcPrime(); proot=GetPrimitiveRoot(); LL aa,bb;//ax=b; LL qq=0; if(X1>=X0) qq=(X1-X0); else qq=P+X1-X0; aa=N-1;bb=GetInd(proot,(K*_Pow(qq,P-2))%P); //printf("a=%d b=%d\n",aa,bb); //printf("proot=%I64d\n",proot); //PS LL a,b;a=_Pow(proot,Linear(aa,bb)); if(X1>=X0*a) qq=X1-X0*a; else qq=P+X1-X0*a; b=(qq%P+P)%P; //PS printf("%I64d %I64d\n",a,b); for(LL lx=0,prc=X0;lx<=N;lx++,prc=(prc*a+b)%P) printf("%I64d ",prc); printf("\n"); return 0; }
沒有留言:
張貼留言