求：$a^{bx \%p}\equiv 1(\mod p)$ 的一个可行的 $x$.
根据欧拉定理，我们知道 $a^{\phi(p)}\equiv 1(\mod p)$
而在 $a^x\equiv 1(\mod p)$ 这个式子中 $x$ 是存在很多个解的.
这些解之间存在着循环节，使得任意解 $x$ 可以被表示成循环节的倍数.
我们设这个循环节为 $cir$.
由于已知 $\phi(p)$ 一定是一个可行解，所以最小循环节一定是 $\phi(p)$ 的约数.
然后我们就可以对 $\phi(p)$ 进行质因数分解来求这个最小循环节 $cir$.
求出来 $cir$ 后，可得 $bx\%p=cir\times y$
你发现因为有这个 $\%p$ 操作，所以会导致这个方程解不出来.
但是好在我们发现，左面那个式子可以变为 $bx+pk=l$，而这个 $l$ 可以被表示为 $i\times gcd(b,p)$
故我们可以将式子变为 $i\times gcd(b,p)=y\times cir$.
然后我们可以对这个求最小通解（因为如果最小解小于 $p$，则一定可以被 $bx\%p$ 表示出来，而且满足 $x\leqslant p-1$ )
这个最小通解是 $c=\frac{cir\times gcd(b,p)}{gcd(gcd(b,p),cir)}$
然后用 exgcd 求一下 $bx+pk=c$ 的 $x$ 的最小正整数解就可以了.
这里可以证明一下为什么只要存在 $bx+pk=c$ 就能保证 $x<p$：
我们可以将 $x$ 表示成 $p+d$ 的形式，那么原式为 $b(p+d)+pk=c$
$\Rightarrow bp+bd+pk=c$
$\Rightarrow b\times d+p\times (k+d)=c$
所以，一旦 $x>p$，我们就可以将一些部分导到 $p\times k$ 那里，以此来实现 $x<p$
#include <string>
#include <ctime>
#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
#define N 10000060
#define ll long long
using namespace std;
void setIO(string s) {
string in=s+".in";
string out=s+".out";
freopen(in.c_str(),"r",stdin);
freopen(out.c_str(),"w",stdout);
}
int H=0;
int fr[N];
int prime[N],vis[N],phi[N];
int answer[N];
int MN[N];
struct data {
int a,b,p,id;
data(int a=0,int b=0,int p=0,int id=0):a(a),b(b),p(p),id(id){}
};
vector<data>G[N];
inline int qpow(int x,int y,int mod) {
int tmp=1;
while(y) {
if(y&1) {
tmp=(ll)tmp*x%mod;
}
x=(ll)x*x%mod;
y>>=1;
}
return tmp;
}
int exgcd(int a,int b,int &x,int &y) {
if(!b) {
x=1,y=0;
return a;
}
int gcd=exgcd(b,a%b,x,y);
int tmp=x;
x=y,y=tmp-(a/b)*y;
return gcd;
}
void Linear_shaker() {
int i,j,cnt=0;
fr[1]=1;
for(i=2;i<N;++i) {
if(!vis[i]) {
prime[++cnt]=i;
phi[i]=i-1;
fr[i]=1;
}
for(j=1;j<=cnt&&prime[j]*i<N;++j) {
vis[i*prime[j]]=1;
if(i%prime[j]) {
phi[i*prime[j]]=phi[i]*(prime[j]-1);
fr[i*prime[j]]=i;
}
else {
fr[i*prime[j]]=i;
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
}
}
}
int main() {
// setIO("mod");
Linear_shaker();
int i=0,j=0,t1,t2,tp=0;
int a,b,p,Mx=0,k;
while(scanf("%d",&a)!=EOF) {
scanf("%d%d",&b,&p);
++i;
++tp;
if(__gcd(a,p)!=1) {
answer[i]=-1;
}
else {
Mx=max(Mx,phi[p]);
G[phi[p]].push_back(data(a,b,p,i));
MN[i]=phi[p];
int w=phi[p],tmp,c=1;
while(w!=1) {
while(w!=1&&qpow(a,fr[w]*c,p)==1) {
w=fr[w];
}
tmp=w/fr[w];
while(w%tmp==0&&w!=1) {
c*=tmp;
w/=tmp;
}
}
MN[i]=c;
}
}
for(i=1;i<=Mx;++i) {
for(j=0;j<G[i].size();++j) {
a=G[i][j].a;
b=G[i][j].b;
p=G[i][j].p;
int x=0,y=0;
int id=G[i][j].id;
int delta=MN[id];
int gcd=__gcd(p,b);
int tmp=(1ll*gcd*delta)/(__gcd(gcd,delta));
if(tmp>=p) {
answer[id]=-1;
}
else {
gcd=exgcd(b,p,x,y);
x=(1ll*x*(tmp/gcd)%(p/gcd)+(p/gcd))%(p/gcd);
answer[id]=x;
}
}
}
for(i=1;i<=tp;++i) printf("%d\n",answer[i]);
return 0;
}
