参考资料
原根
如果g是m的原根，对于任意一个数x(x<m)，都可以找到一个I(x) 小于等于 φ(m),使得 gI(x) = x ,I(x)称为x的指标。也就是说在整数模m乘法群
(所有与 m互素的正整数构成的等价类构成的乘法群)中，g是这个乘法群的生成元。它的生成元的个数就是它的可逆元个数，即φ(φ(m))个。也就是说，m有φ(φ(m))个原根
除了直接运算以外，至今还没有一个办法可以找到模特定m时的原根，但假如已知模m有一个原根，则可找出它其他的原根。
离散对数
从表面通俗地理解，和普通对数一样，两边同时取log，这样就可以将幂次变成多项式，在离散对数中，做了这样的变换后,模m则变成了模φ(m)。于是我们就可以将一些乘法操作变成加法操作，来简化或者解决一些运算问题.
求质数ｐ的最小原根
step1 : 筛素数
step2 : 分解φ(p)
step3 : 按定义判断φ(p)每个质因数d是不是原根
bool not_prime[N+5]; int prime[N+5],tot_prime,pr[N+5],tot_pr; ll ksm(ll a,ll b,ll mod){ ll res=1; while(b){ if(b&1) res = (res*a)%mod; a = (a*a)%mod; b>>=1; } return res; } void sss(){ // 线性筛 for(int i=2;i<N;i++){ if(!not_prime[i]) prime[++tot_prime]=i; for(int j=1;j<=tot_prime;j++){ int t = i*prime[j];if(t>=N) break;not_prime[t]=true; if(i%prime[j]==0)break; } } } void get(int p){ // 分解phi(p)的质因数存到pr[]中 int temp=p-1; for(int i=1;i<=tot_prime&&prime[i]*prime[i]<=temp;i++) if(temp%prime[i]==0){ while(temp%prime[i]==0) temp /= prime[i]; pr[++tot_pr] = prime[i]; } if(temp>1) pr[++tot_pr]=temp; } bool check(int d,int p){ //判断d是不是p的原根 if (std::__gcd(d,p)!=1)return 0; for(int i=1;i<=tot_pr;i++)if(ksm(d,(p-1)/pr[i],p)==1)return 0; return 1; } int get_g(int p) {// 求p的原根,p是质数 int g=1; sss();get(p); for(int i=2;i<p;i++) if(check(i,p)) {g=i;break;} return g; }
BSGS算法
给定a,b,p, 求解 
,p是个质数
时间复杂度 O(√p)
由上面的原根概念可知，这里需要求b的指标I(b)
const int MOD=76543; ll hs[MOD],head[MOD],next[MOD],id[MOD],top; void insert(ll x,int y){ ll k=x%MOD; hs[++top]=x,id[top]=y,next[top]=head[k],head[k]=top; } ll find(int x){ ll k=x%MOD; for(int i=head[k];i!=0;i=next[i]) if(hs[i]==x)return id[i]; return -1; } ll bsgs(ll a,ll n,ll b) { if(b==1)return 0; int m=sqrt(1.0*n),j; ll x=1,p=1,i; for(i=0;i<m;i++,p=p*a%n)insert(p*b%n,i); for(i=m;;i+=m){ if((j=find(x=x*p%n))!=-1) return i-j; if(i>n) break; } return -1; }
扩展BSGS算法
给定a,b,p, 求解 但是不保证p是质数
#include<bits/stdc++.h> typedef long long ll; const int N=1e5; ll p,k,a; bool not_prime[N+5]; int prime[N+5],tot_prime,pr[N+5],tot_pr; const int MOD=76543; ll hs[MOD],head[MOD],next[MOD],id[MOD],top; ll ksm(ll a,ll b,ll mod){ ll res=1; while(b){ if(b&1) res = (res*a)%mod; a = (a*a)%mod; b>>=1; } return res; } void insert(ll x,int y){ ll k=x%MOD; hs[++top]=x,id[top]=y,next[top]=head[k],head[k]=top; } ll find(int x){ ll k=x%MOD; for(int i=head[k];i!=0;i=next[i]) if(hs[i]==x)return id[i]; return -1; } ll bsgs(ll a,ll n,ll b) { if(b==1)return 0; int m=sqrt(1.0*n),j; ll x=1,p=1,i; for(i=0;i<m;i++,p=p*a%n)insert(p*b%n,i); for(i=m;;i+=m){ if((j=find(x=x*p%n))!=-1) return i-j; if(i>n) break; } return -1; } void sss(){ // 线性筛 for(int i=2;i<N;i++){ if(!not_prime[i]) prime[++tot_prime]=i; for(int j=1;j<=tot_prime;j++){ int t = i*prime[j];if(t>=N) break;not_prime[t]=true; if(i%prime[j]==0)break; } } } void get(){ // 分解phi(p)的质因数存到pr[]中 int temp=p-1; for(int i=1;i<=tot_prime&&prime[i]*prime[i]<=temp;i++) if(temp%prime[i]==0){ while(temp%prime[i]==0) temp /= prime[i]; pr[++tot_pr] = prime[i]; } if(temp>1) pr[++tot_pr]=temp; } bool check(ll d){ if (std::__gcd(d,p)!=1)return 0; for(int i=1;i<=tot_pr;i++)if(ksm(d,(p-1)/pr[i],p)==1)return 0; return 1; } ll exgcd(ll m,ll &x,ll n,ll &y){ ll x1=0,y1=1,x0=1,y0=0; ll r=(m%n+n)%n;ll q=(m-r)/n; x=0;y=1; while(r){ x=x0-q*x1;y=y0-q*y1; x0=x1;y0=y1;x1=x;y1=y; m=n;n=r;r=m%n; q=(m-r)/n; } return n; } void solve(ll p,ll k,ll a){ std::set<int>ans; int g=1; sss();get(); for(int i=2;i<p;i++)if(check(i)){g=i;break;} ll Ia=bsgs(g,p,a); ll ar,br; ll gcd=exgcd(k,ar,p-1,br); if(Ia%gcd==0){ ll sp=(p-1)/gcd; ar=((ar*Ia/gcd%sp)%sp+sp)%sp; if(ar==0)ar=sp; for(int i=ar;i<p;i+=sp)ans.insert(ksm(g,i,p)); } int m = ans.size(); printf("%d ",m); for(std::set<int>::iterator it=ans.begin();it!=ans.end();it++) printf("%d ",(*it)); } int main(){ scanf("%lld%lld%lld",&p,&k,&a); solve(p,k,a); }
BSGS算法进阶
求解xk ≡ a(mod p) p是一个质数
该模型可以通过一系列的转化为成BSGS算法的常规模型
gI(a) ≡ a (mod p) BSGS解得I(a)
gI(x) ≡ x (mod p)
原方程变形为 (gI(x))k ≡ gI(a) (mod p)
两边取离散对数 得到 k* I(x) ≡ I(a) (mod φ(p))
即 k*I(x) + φ(p) * yy = I(a) 用扩展欧几里得算法解得一个I(x)
如何获得所有解呢？
对于关于x，y的方程a * x + b * y = gcd 来说，让x增加b/gcd，让y减少a/gcd，等式两边还相等。
对于方程a * x + b * y = c 来说，如果c % gcd(a,b) != 0，则无解。
如果c % gcd(a,b) == 0 且 c / gcd(a,b) = t，那么求出方程 a * x + b * y = gcd(a,b)的所有解x,y，将x,y乘上t，对应的x’,y’即是方程a * x + b * y = c的解
每得到一个I(x) 计算出x ，I(x)≤φ(p)
#include<bits/stdc++.h> typedef long long ll; const int N=1e5; ll p,k,a; bool not_prime[N+5]; int prime[N+5],tot_prime,pr[N+5],tot_pr; const int MOD=76543; ll hs[MOD],head[MOD],next[MOD],id[MOD],top; ll ksm(ll a,ll b,ll mod){ ll res=1; while(b){ if(b&1) res = (res*a)%mod; a = (a*a)%mod; b>>=1; } return res; } void insert(ll x,int y){ ll k=x%MOD; hs[++top]=x,id[top]=y,next[top]=head[k],head[k]=top; } ll find(int x){ ll k=x%MOD; for(int i=head[k];i!=0;i=next[i]) if(hs[i]==x)return id[i]; return -1; } ll bsgs(ll a,ll n,ll b) { if(b==1)return 0; int m=sqrt(1.0*n),j; ll x=1,p=1,i; for(i=0;i<m;i++,p=p*a%n)insert(p*b%n,i); for(i=m;;i+=m){ if((j=find(x=x*p%n))!=-1) return i-j; if(i>n) break; } return -1; } void sss(){ // 线性筛 for(int i=2;i<N;i++){ if(!not_prime[i]) prime[++tot_prime]=i; for(int j=1;j<=tot_prime;j++){ int t = i*prime[j];if(t>=N) break;not_prime[t]=true; if(i%prime[j]==0)break; } } } void get(){ // 分解phi(p)的质因数存到pr[]中 int temp=p-1; for(int i=1;i<=tot_prime&&prime[i]*prime[i]<=temp;i++) if(temp%prime[i]==0){ while(temp%prime[i]==0) temp /= prime[i]; pr[++tot_pr] = prime[i]; } if(temp>1) pr[++tot_pr]=temp; } bool check(ll d){ if (std::__gcd(d,p)!=1)return 0; for(int i=1;i<=tot_pr;i++)if(ksm(d,(p-1)/pr[i],p)==1)return 0; return 1; } ll exgcd(ll m,ll &x,ll n,ll &y){ ll x1,y1,x0,y0; x0=1;y0=0; x1=0;y1=1; ll r=(m%n+n)%n; ll q=(m-r)/n; x=0;y=1; while(r){ x=x0-q*x1; y=y0-q*y1; x0=x1;y0=y1; x1=x;y1=y; m=n;n=r;r=m%n; q=(m-r)/n; } return n; } void solve(ll p,ll k,ll a){ std::set<int>ans; int g=1; sss();get(); for(int i=2;i<p;i++)if(check(i)){g=i;break;} ll Ia=bsgs(g,p,a); ll ar,br; ll gcd=exgcd(k,ar,p-1,br); if(Ia%gcd==0){ ll sp=(p-1)/gcd; ar=((ar*Ia/gcd%sp)%sp+sp)%sp; if(ar==0)ar=sp; for(int i=ar;i<p;i+=sp)ans.insert(ksm(g,i,p)); } int m = ans.size(); printf("%d ",m); for(std::set<int>::iterator it=ans.begin();it!=ans.end();it++) printf("%d ",(*it)); } int main(){ scanf("%lld%lld%lld",&p,&k,&a); solve(p,k,a); }