參考資料
原根
如果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\n",m); for(std::set<int>::iterator it=ans.begin();it!=ans.end();it++) printf("%d\n",(*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\n",m); for(std::set<int>::iterator it=ans.begin();it!=ans.end();it++) printf("%d\n",(*it)); } int main(){ scanf("%lld%lld%lld",&p,&k,&a); solve(p,k,a); }