参考资料：OI-Wiki
\(\text{BSGS}\) (\(\text{Big-Step-Giant-Step}\))
最基本的线性同余方程：
\[ax\equiv b\pmod p \]可以转化成不定方程：
\[ax+py=b \]使用扩展欧几里得算法解出 \(x\) 值得解。
而面对形如：
\[a^x\equiv b\pmod p \]的求指数方程，在保证 \(\gcd(a,p)=1\) 的情况下，需要用到 \(\text{BSGS}\)。
思考暴力的做法，可以直接暴力枚举左侧的指数，而右侧没有任何变化。
使用类似 \(\text{meet-in-the-middle}\) 的优化策略，将 \(x\) 表示为 \(A\left\lfloor\sqrt{p}\right\rfloor-B\)，这样得到新的形式：
\[a^{A\left\lfloor\sqrt{p}\right\rfloor}\equiv ba^{B}\pmod p \]左侧指数与右侧分别只有 \(O(\sqrt{p})\) 种取值，可以预处理出右侧的所有取值，再枚举左侧取值并套用 unordered_map 判断是否存在。
注意到这里移项成立的充要条件是 \(a^B\) 存在模 \(p\) 意义下的逆元，可得条件应为 \(\gcd(a,p)=1\)。
点击查看代码
int a,b,p;
unordered_map<int,int> mp;
inline void BSGS(){
if(1%p==b%p) return printf("0\n"),void();
int B=sqrt(p),cntB=(p-1)/B+1;
int bs=1,now=b;
for(int i=0;i<B;++i){
mp[now]=i;
bs=1ll*bs*a%p,now=1ll*now*a%p;
}
now=bs;
for(int i=1;i<=cntB;++i){
if(mp.count(now)){
int ans=i*B-mp[now];
printf("%d\n",ans);
return;
}
now=1ll*now*bs%p;
}
printf("no solution\n");
}
int main(){
p=read(),a=read(),b=read();
BSGS();
return 0;
}
\(\text{exBSGS}\)
扩展算法就是考虑不互质的情况。
不互质时，设 \(\gcd(a,p)=d\)，我们同时消去一个 \(d\)，得到：
\[\dfrac{a}{d}\times a^{x-1}\equiv \dfrac{b}{d}\pmod {\dfrac{p}{d}} \]当然这步成立的前提是 \(d\mid b\)，这个方程的底数仍为 \(a\)，模数发生了改变，但仍不保证二者互质，继续进行这个操作，直到 \(\gcd\left(a,\dfrac{p}{\prod_{i=1}^k d_i}\right)=1\)，方程就变成：
\[\dfrac{a^k}{\prod_{i=1}^k d_i}\times a^{x-k}\equiv \dfrac{b}{\prod_{i=1}^k d_i}\pmod {\dfrac{p}{\prod_{i=1}^k d_i}} \]此时左侧的系数显然和模数互质，是存在逆元的，移项过去之后只需解：
\[a^{x-k}\equiv \dfrac{b}{a^k}\pmod {\dfrac{p}{\prod_{i=1}^k d_i}} \]套用 \(\text{BSGS}\) 即可。
注意到原方程的解是由此方程的解加上 \(k\) 而得，因此 \([0,k)\) 部分无法计算到，只需要在刚开始每次消去 \(d\) 时暴力判断一下。
点击查看代码
int a,b,p;
inline int q_pow(int A,int B,int MOD){
int res=1;
while(B){
if(B&1) res=1ll*res*A%MOD;
A=1ll*A*A%MOD;
B>>=1;
}
return res;
}
int gcd(int A,int B){
if(!B) return A;
return gcd(B,A%B);
}
int exgcd(int A,int B,int &X,int &Y){
if(!B){
X=1,Y=0;
return A;
}
int GCD=exgcd(B,A%B,Y,X);
Y-=A/B*X;
return GCD;
}
unordered_map<int,int> mp;
inline void BSGS(int A,int B,int P,int K){
if(1%P==B%P) return printf("%d\n",K),void();
mp.clear();
int siz=sqrt(P),cnt=(P-1)/siz+1;
int bs=1,now=B;
for(int i=0;i<siz;++i){
mp[now]=i;
bs=1ll*bs*A%P,now=1ll*now*A%P;
}
now=bs;
for(int i=1;i<=cnt;++i){
if(mp.count(now)){
int ans=i*siz-mp[now]+K;
printf("%d\n",ans);
return;
}
now=1ll*now*bs%P;
}
printf("No Solution\n");
}
inline void exBSGS(int A,int B,int P){
int D=1,K=0;
int now=1;
while(1){
int d=gcd(A,P);
if(d==1) break;
if(now==b) return printf("%d\n",K),void();
now=1ll*now*a%p;
++K;
D=D*d,P/=d;
if(B%D) return printf("No Solution\n"),void();
}
int pw=q_pow(A,K,P),inv,Y;
exgcd(pw,P,inv,Y);
inv=(inv%P+P)%P;
BSGS(A,1ll*B*inv%P,P,K);
}
int main(){
// freopen("exbsgs.in","r",stdin);
// freopen("exbsgs.out","w",stdout);
while(1){
a=read(),p=read(),b=read();
if(!a&&!p&&!b) break;
b%=p;
exBSGS(a,b,p);
}
return 0;
}
阶与原根
定义和性质
关于阶的性质在此不多赘述，只需知道 \(n\) 为 \(a\) 模 \(m\) 的阶，表示 \(n\) 为最小正整数使得 \(a^n\equiv 1\pmod m\)，记作 \(\delta_{m}(a)\)。
定义原根 \(g\) 为满足 \(\gcd(g,m)=1\) 且 \(\delta_{m}(g)=\varphi(m)\) 的正整数。
原根最有用的性质是：当 \(m\) 为质数时，对于 \(i\in [0,m)\)，\(g^i\bmod m\) 两两不同且取遍 \([0,m)\)。
原根存在定理
一个数 \(m\) 存在原根，当且仅当 \(m=2,4,p^{\alpha},2p^{\alpha}\)，其中 \(p\) 是不为 \(2\) 的素数。
证明略。
原根的个数
由于 \(g^k\) 在模 \(m\) 下可以表示任意一个数，于是有：
\[\delta_{m}(g^k)=\dfrac{\delta_{m}(g)}{\gcd(\delta_{m}(g),k)}=\dfrac{\varphi(m)}{\gcd(\varphi(m),k)} \]（这里运用了阶的性质。）
于是当且仅当 \(\gcd(\varphi(m),k)=1\) 时，\(g^k\) 同为原根。
由欧拉定理知，\(g\) 的 \([1,\varphi(m)]\) 次幂在模 \(m\) 下形成循环节，于是有意义的原根均在 \([1,\varphi(m)]\) 中，这之间与 \(\varphi(m)\) 互质的 \(k\)，自然有 \(\varphi(\varphi(m))\) 个。
原根判定定理
用通俗易懂的话来描述原根：
由于 \(\gcd(g,m)=1\)，因此 \(g^{\varphi(m)}\equiv 1\pmod m\) 一定成立，如果 \(\varphi(m)\) 同时是最小的满足此条件的正整数（也就是阶），那么 \(g\) 就是原根。
所以判定原根实际上是在保证互质的基础上判定最小，也就是任意 \(\varphi(m)\) 的一个真因子 \(d\) 都不能使 \(g^{d}\equiv 1\pmod m\)。
一个简洁而快速的检查方法是，令 \(d\) 取遍 \(\varphi(m)\) 除去一个质因子后的值，这样即可覆盖所有的真因子。
求底数的同余方程
形如：
\[x^a\equiv b\pmod p \]其中 \(p\) 是质数的同余方程。
因为 \(p\) 是质数也就一定有原根 \(g\)，也就可以表示任意一个 \(x\)，即 \(x\equiv g^k\pmod p\)。
于是有：
\[\begin{aligned} (g^k)^a&\equiv b\pmod p\\ (g^a)^k&\equiv b\pmod p \end{aligned}\]这样求底数的问题就被转化成求指数了，并且模数是质数，可以直接 \(\text{BSGS}\) 计算。
点击查看代码
int a,b,p;
int phi,g;
vector<int> pr;
inline void get_pr(){
int tmp=phi;
for(int i=2;i*i<=tmp;++i){
if(tmp%i==0){
pr.push_back(i);
while(tmp%i==0) tmp/=i;
}
}
if(tmp!=1) pr.push_back(tmp);
}
inline int q_pow(int A,int B,int P){
int res=1;
while(B){
if(B&1) res=1ll*res*A%P;
A=1ll*A*A%P;
B>>=1;
}
return res;
}
unordered_map<int,int> mp;
vector<int> ans;
inline void BSGS(int A,int B,int P){
if(1%P==B%P) ans.push_back(0);
int siz=sqrt(P),cnt=(p-1)/siz+1;
int bs=1,now=B;
for(int i=0;i<siz;++i){
mp[now]=i;
bs=1ll*bs*A%P,now=1ll*now*A%P;
}
now=bs;
for(int i=1;i<=cnt;++i){
if(mp.count(now)) ans.push_back(i*siz-mp[now]);
now=1ll*now*bs%P;
}
}
int main(){
p=read(),a=read(),b=read();
phi=p-1;
get_pr();
for(g=2;;++g){
bool pd=1;
for(int i:pr){
if(q_pow(g,phi/i,p)==1){
pd=0;
break;
}
}
if(pd) break;
}
BSGS(q_pow(g,a,p),b,p);
for(int i=0;i<ans.size();++i){
int now=ans[i];
ans[i]=q_pow(g,now,p);
}
sort(ans.begin(),ans.end());
ans.erase(unique(ans.begin(),ans.end()),ans.end());
printf("%d\n",(int)ans.size());
for(int i:ans) printf("%d ",i);
printf("\n");
return 0;
}