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.