HDU3930 （原根）

给定方程 X^A = B (mol C)  ，求 在[0,C) 中所有的解 ， 并且C为质数。

设 rt 为 C 的原根 , 则 X = rt^x  （这里相当于求 A^x =B (mol C) 用大步小步算法即可）

那么 ( rt^x ) ^ A = b (mol C)

　　 rt^Ax = b (mol C)

由费马小定理， 设 Ax = (C-1)*y +t1   ---------------- ( * ) 

可得  rt^t1 =b ( mod C)

这里运用大步小步算法可以计算出 t1 。

得到 t1 后反代会 （*）式 ， 利用扩展欧几里得求出符合条件的x解。

由于此方程相当于解 Ax mod (C-1) = t1 , 共用 gcd ( a , C-1 ) 组解。

最后用快速幂计算出所有的X解即可。 

  1 const maxn=1000008;
  2       maxh=1000007;
  3 var a,b,c,rt,t1,t2,x,y,d:int64;
  4     i:longint;
  5     ans,pm,pri:array[0..maxn*2] of int64;
  6     pd:array[0..maxn*2] of boolean;
  7     cnt,nm:longint;
  8     h:array[0..maxh,1..2] of int64;
  9 procedure init;
 10 var i,j:longint;
 11 begin
 12     fillchar(pd,sizeof(pd),false);
 13     i:=2; nm:=0;
 14     while i<=maxn do
 15     begin
 16         inc(nm);
 17         pm[nm]:=i;
 18         j:=i;
 19         while j<=maxn do
 20         begin
 21             pd[j]:=true;
 22             j:=j+i;
 23         end;
 24         while pd[i] do inc(i);
 25     end;
 26 end;
 27 function pow(x,y,p:int64):int64;
 28 var sum:int64;
 29 begin
 30     x:=x mod p;
 31     sum:=1;
 32     while y>0 do
 33     begin
 34         if y and 1 =1 then sum:=sum*x mod p;
 35         x:=x*x mod p;
 36         y:=y >> 1;
 37     end;
 38     exit(sum);
 39 end;
 40 procedure divide(n:int64);
 41 var i:longint;
 42 begin
 43     cnt:=0;
 44     i:=1;
 45     while pm[i]*pm[i]<=n do
 46     begin
 47         if n mod pm[i]=0 then
 48         begin
 49             inc(cnt);
 50             pri[cnt]:=pm[i];
 51             while n mod pm[i]=0 do n:=n div pm[i];
 52         end;
 53         inc(i);
 54     end;
 55     if n>1 then
 56     begin
 57         inc(cnt);
 58         pri[cnt]:=n;
 59     end;
 60 end;
 61 function findrt(p:int64):int64;
 62 var g,t:int64;
 63     flag:boolean;
 64 begin
 65     divide(p-1);
 66     g:=2;
 67     while true do
 68     begin
 69         flag:=true;
 70         for i:=1 to cnt do
 71         begin
 72             t:=(p-1) div pri[i];
 73             if pow(g,t,p)=1 then
 74             begin
 75                 flag:=false;
 76                 break;
 77             end;
 78         end;
 79         if flag then exit(g);
 80         inc(g);
 81     end;
 82 end;
 83 procedure insert(x,y:int64); inline;
 84 var hash:int64;
 85 begin
 86     hash:=x mod maxh;
 87     while (h[hash,1]<>x) and (h[hash,1]<>0) do hash:=(hash+1) mod maxh;
 88     h[hash,1]:=x;
 89     h[hash,2]:=y;
 90 end;
 91 function find(x:int64):int64; inline;
 92 var hash:int64;
 93 begin
 94     hash:=x mod maxh;
 95     while (h[hash,1]<>x) and (h[hash,1]<>0) do hash:=(hash+1) mod maxh;
 96     if h[hash,1]=0 then exit(-1) else exit(h[hash,2]);
 97 end;
 98 function work(a,b,p:int64):int64;
 99 var j,m,x,cnt,ans,t:int64;
100     i:longint;
101 begin
102     ans:=2000000000;
103     m:=trunc(sqrt(p))+1;
104     x:=pow(a,m,p);
105     j:=1;
106     for i:=1 to m do
107     begin
108         j:=j*x mod p;
109         if find(j)=-1 then insert(j,i);
110     end;
111     j:=1;
112     for i:=0 to m-1 do
113     begin
114         t:=find(j*b mod p);
115         if t<>-1 then
116         begin
117             cnt:=m*t-i;
118             if cnt<ans then ans:=cnt;
119         end;
120         j:=j*a mod p;
121     end;
122     exit(ans);
123 end;
124 function gcd(x,y:int64):int64;
125 begin
126     if y=0 then exit(x) else exit(gcd(y,x mod y));
127 end;
128 procedure exgcd(a,b:int64;var x,y:int64);
129 var t:int64;
130 begin
131     if b=0 then
132     begin
133         x:=1;
134         y:=0;
135         exit;
136     end;
137     exgcd(b,a mod b,x,y);
138     t:=x;
139     x:=y;
140     y:=t-a div b*y;
141 end;
142 procedure swap(var a,b:int64); inline;
143 var c:longint;
144 begin
145     c:=a; a:=b; b:=c;
146 end;
147 procedure sort(l,r:int64);
148 var i,j,x:int64;
149 begin
150     i:=l; j:=r; x:=ans[(l+r) div 2];
151     while i<=j do
152     begin
153         while ans[i]<x do inc(i);
154         while x<ans[j] do dec(j);
155         if i<=j then
156         begin
157             swap(ans[i],ans[j]);
158             inc(i); dec(j);
159         end;
160     end;
161     if l<j then sort(l,j);
162     if i<r then sort(i,r);
163 end;
164 begin
165     init;
166     readln(b,a,c);
167     rt:=findrt(c); 
168     t1:=work(rt,b,c);
169     t2:=c-1;
170     d:=gcd(a,t2);
171     if t1 mod d<>0 then
172     begin
173         writeln(0);
174         exit;
175     end;
176     exgcd(a,t2,x,y);
177     t1:=t1 div d;
178     t2:=t2 div d;
179     ans[1]:=((x*t1 mod t2)+ t2) mod t2;
180     for i:=2 to d do ans[i]:=ans[i-1]+t2;
181     for i:=1 to d do ans[i]:=pow(rt,ans[i],c);
182     sort(1,d);
183     writeln(d);
184     for i:=1 to d-1 do write(ans[i],' ');
185     writeln(ans[d]);
186 end.

 

转载于:https://www.cnblogs.com/rpSebastian/p/4571677.html