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解即可。

const maxn=1000008;
maxh=1000007;
var a,b,c,rt,t1,t2,x,y,d:int64;
i:longint;
ans,pm,pri:array[0..maxn*2] of int64;
pd:array[0..maxn*2] of boolean;
cnt,nm:longint;
h:array[0..maxh,1..2] of int64;
procedure init;
var i,j:longint;
begin
fillchar(pd,sizeof(pd),false);
i:=2; nm:=0;
while i<=maxn do
begin
inc(nm);
pm[nm]:=i;
j:=i;
while j<=maxn do
begin
pd[j]:=true;
j:=j+i;
end;
while pd[i] do inc(i);
end;
end;
function pow(x,y,p:int64):int64;
var sum:int64;
begin
x:=x mod p;
sum:=1;
while y>0 do
begin
if y and 1 =1 then sum:=sum*x mod p;
x:=x*x mod p;
y:=y >> 1;
end;
exit(sum);
end;
procedure divide(n:int64);
var i:longint;
begin
cnt:=0;
i:=1;
while pm[i]*pm[i]<=n do
begin
if n mod pm[i]=0 then
begin
inc(cnt);
pri[cnt]:=pm[i];
while n mod pm[i]=0 do n:=n div pm[i];
end;
inc(i);
end;
if n>1 then
begin
inc(cnt);
pri[cnt]:=n;
end;
end;
function findrt(p:int64):int64;
var g,t:int64;
flag:boolean;
begin
divide(p-1);
g:=2;
while true do
begin
flag:=true;
for i:=1 to cnt do
begin
t:=(p-1) div pri[i];
if pow(g,t,p)=1 then
begin
flag:=false;
break;
end;
end;
if flag then exit(g);
inc(g);
end;
end;
procedure insert(x,y:int64); inline;
var hash:int64;
begin
hash:=x mod maxh;
while (h[hash,1]<>x) and (h[hash,1]<>0) do hash:=(hash+1) mod maxh;
h[hash,1]:=x;
h[hash,2]:=y;
end;
function find(x:int64):int64; inline;
var hash:int64;
begin
hash:=x mod maxh;
while (h[hash,1]<>x) and (h[hash,1]<>0) do hash:=(hash+1) mod maxh;
if h[hash,1]=0 then exit(-1) else exit(h[hash,2]);
end;
function work(a,b,p:int64):int64;
var j,m,x,cnt,ans,t:int64;
i:longint;
begin
ans:=2000000000;
m:=trunc(sqrt(p))+1;
x:=pow(a,m,p);
j:=1;
for i:=1 to m do
begin
j:=j*x mod p;
if find(j)=-1 then insert(j,i);
end;
j:=1;
for i:=0 to m-1 do
begin
t:=find(j*b mod p);
if t<>-1 then
begin
cnt:=m*t-i;
if cnt<ans then ans:=cnt;
end;
j:=j*a mod p;
end;
exit(ans);
end;
function gcd(x,y:int64):int64;
begin
if y=0 then exit(x) else exit(gcd(y,x mod y));
end;
procedure exgcd(a,b:int64;var x,y:int64);
var t:int64;
begin
if b=0 then
begin
x:=1;
y:=0;
exit;
end;
exgcd(b,a mod b,x,y);
t:=x;
x:=y;
y:=t-a div b*y;
end;
procedure swap(var a,b:int64); inline;
var c:longint;
begin
c:=a; a:=b; b:=c;
end;
procedure sort(l,r:int64);
var i,j,x:int64;
begin
i:=l; j:=r; x:=ans[(l+r) div 2];
while i<=j do
begin
while ans[i]<x do inc(i);
while x<ans[j] do dec(j);
if i<=j then
begin
swap(ans[i],ans[j]);
inc(i); dec(j);
end;
end;
if l<j then sort(l,j);
if i<r then sort(i,r);
end;
begin
init;
readln(b,a,c);
rt:=findrt(c);
t1:=work(rt,b,c);
t2:=c-1;
d:=gcd(a,t2);
if t1 mod d<>0 then
begin
writeln(0);
exit;
end;
exgcd(a,t2,x,y);
t1:=t1 div d;
t2:=t2 div d;
ans[1]:=((x*t1 mod t2)+ t2) mod t2;
for i:=2 to d do ans[i]:=ans[i-1]+t2;
for i:=1 to d do ans[i]:=pow(rt,ans[i],c);
sort(1,d);
writeln(d);
for i:=1 to d-1 do write(ans[i],' ');
writeln(ans[d]);
end.