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.
