题解太高深，自己整理一份。
阶的定义： 若 \(\gcd(a,m)=1\)，则使 \(a^l\equiv1\pmod{m}\) 的最小的 \(l\) 称为 \(a\) 关于模 \(m\) 的阶，记为 \(\operatorname{ord}_m a\)。
阶的性质：
根据欧拉定理，\(a^{\varphi(m)}\equiv 1\pmod{m}\)，所以 \(\operatorname{ord}_m a\mid \varphi(m)\)。
进一步的，若 \(a^x\equiv1\pmod{m}\)，则 \(\operatorname{ord}_m a\mid x\)。
由此可见，阶相当于最小循环节。
证明： \(\operatorname{ord}_m a^t=\frac{\operatorname{ord}_m a}{\gcd(\operatorname{ord}_m a,t)}\)。
设 \(\operatorname{ord}_m a^t=l\)，则 \(a^{t\cdot l}\equiv1\pmod{m}\)。
因为 \(\operatorname{ord}_m a\mid t\cdot l\)，\(t\mid t\cdot l\)，所以 \(t\cdot l=\operatorname{lcm}(\operatorname{ord}_m a,t)\)。
所以 \(l=\frac{\operatorname{lcm}(\operatorname{ord}_m a,t)}{t}=\frac{\operatorname{ord}_m a\cdot t}{\gcd(\operatorname{ord}_m a,t)\cdot t}=\frac{\operatorname{ord}_m a}{\gcd(\operatorname{ord}_m a,t)}\)。
原根的定义： 当 \(\operatorname{ord}_m a=\varphi(m)\) 时，称 \(a\) 为关于 \(m\) 的原根。
原根的性质：
只有 \(2,4,p_i,2\cdot p_i\) 有原根，其中 \(p_i\) 为奇素数。
当 \(g\) 为原根时 \(g,g^2,g^3,...,g^{\varphi(m)}\) 构成模 \(m\) 意义下的既约剩余系，而因为 \(\gcd(g,m)=1\)，所以 \(\gcd(g^x,m)=1\)，且剩余系中的各个数互不相等。
考虑反证法。若有两个数 \(g^x,g^y\) 相等，则说明 \(g^{y-x}\equiv1\pmod{m}\)，由于 \(y-x<\varphi(m)\)，所以 \(g\) 不为原根，矛盾。
那么，可能的原根只能在 \(g,g^2,g^3,...,g^{\varphi(m)}\) 中。
如果 \(g^t\) 满足 \(\operatorname{ord}_m g^t=\varphi(m)\)，则 \(\frac{\operatorname{ord}_m g}{\gcd(\operatorname{ord}_m g,t)}=\varphi(m)\)，即 \(\frac{\varphi(m)}{\gcd(\varphi(m),t)}=\varphi(m)\)。所以当且 \(\gcd(t,\varphi(m))=1\) 时，\(g^t\) 为 \(m\) 的原根。
如何求 \(g\):
显然，\(g^{\varphi(m)}\equiv1\pmod{m}\) 且不存在 \(x<\varphi(m)\) 满足 \(g^x\equiv1\pmod{m}\)。
而如果存在这样的 \(x\)，则必然有 \(x\mid \varphi(m)\)。
所以我们可以通过质因数分解找到 \(\varphi(m)\) 的所有质因子 \(p_i\)，枚举每一个 \(p_i\) 判断 \(g^{\frac{\varphi(m)}{p_i}}\equiv1\pmod{m}\) 即可。
具体实现：
-
线性筛素数，预处理 \(\varphi\) 函数，\(O(n)\)。
-
分解 \(\varphi(n)\)，\(O(\sqrt n)\)。
-
枚举找到 \(g\)，\(O(n^{0.25}\log^2 n)\)。
-
找出所有原根并排序，\(O(n\log n)\)。
综上，复杂度 \(O(Tn\log n)\)。
code:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+5;
int T,n,d,p[N],phi[N],cnt=0,tot=0,ans[N];
bool vis[N],rt[N];
vector<int>yin;
void init(){
phi[1]=1;
for(ll i=2;i<N;++i){
if(!vis[i]){
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*p[j]<N;++j){
vis[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
rt[2]=rt[4]=1;
for(int i=2;i<=cnt;++i){
for(ll j=p[i];j<N;j*=p[i])rt[j]=1;
for(ll j=p[i]*2;j<N;j*=p[i])rt[j]=1;
}
}
void fen(int x){
for(ll i=2;i*i<=x;++i){
if(x%i==0){
yin.push_back(i);
while(x%i==0)x/=i;
}
}
if(x>1)yin.push_back(x);
}
int ksm(ll x,ll y){
ll res=1;
while(y){
if(y&1)res=res*x%n;
y>>=1;
x=x*x%n;
}
return res;
}
int getg(){
for(int i=1;i<=n;++i){
if(ksm(i,phi[n])!=1)continue;
bool flag=1;
for(int j=0;j<yin.size();++j){
if(ksm(i,phi[n]/yin[j])==1){
flag=0;
break;
}
}
if(flag)return i;
}
}
int gcd(int x,int y){return y?gcd(y,x%y):x;}
int main(){
cin>>T;
init();
while(T--){
scanf("%d%d",&n,&d);yin.clear();tot=0;
if(rt[n]){
fen(phi[n]);
ll g=getg();
ll fac=1;
for(int i=1;i<=phi[n];++i){
fac=fac*g%n;
if(gcd(i,phi[n])==1)ans[++tot]=fac;
}
sort(ans+1,ans+tot+1);
}
printf("%d\n",tot);
for(int i=d;i<=tot;i+=d)printf("%d ",ans[i]);
cout<<endl;
}
return 0;
}