以前基本沒有了解原根相關的一塊內容，最近正好碰到了這個題，于是寫一篇博客記錄一下。

這道題的總體思路就是比較明顯，就是先算出$a^p=x$對于每個$x$的解的個數，然后NTT算一下即可。

先來講一下怎么求歐拉函數$?(x)$，即1到$x?1$中與$x$互素的數的個數，我們一般在線性篩里面通過遞推來求出這個函數。也就是我們要找出$?(i?p[j])$和$?(i),?(p[j]),i,p[j]$之間的關系。我們把與$?(i)$所代表的集合寫在一行里面，然后第二行寫上前一行的數加上$i$，一直寫$p[j]$行。

現在舉兩個例子。

例1：

i%p[j]!=0
p[j]=3,i=8//雖然實際不會出現這種，但是對于證明來講問題不大
1  3  5  7
9  11 13 15
17 19 21 23然后phi[24]表示的集合為
1       5   711  13
17 19       23  

例2：

i%p[j]==0
p[j]=2,i=10
1  3  7  9
11 13 17 19
phi[i*p[j]]=8

在例1中，因為$i\mathrm{mod}p[j]\ne 0$，也就是$?(i)$所表示的集合中存在$p[j]$的倍數。然后你可以發現，每一列里面恰好就有一個是$p[j]$的倍數，我們不妨設某一列第一個數為$x$，且$y=x\mathrm{mod}p[j]$，那么因為$p[j]$為質數，所以$y\mathrm{mod}p[j],y+i\mathrm{mod}p[j],y+2i\mathrm{mod}p[j],...,y+(p[j]?1)i\mathrm{mod}p[j]$中$0,1,...,p[j]?1$恰好每個一個，因此可以知道$phi[i?p[j]]=phi[i]?(p[j]?1)(i\mathrm{mod}p[j]\ne 0)$。

在例2中，因為$i\mathrm{mod}p[j]=0$，因此$?(i)$所表示的集合里面沒有$p[j]$的倍數，因此第一行里面沒有被刪掉的，因為相鄰兩行之間的差都是$p[j]$的倍數，所以不會有任何的數為$p[j]$的倍數。

然后就是原根部分的知識。

如果$x$是模$P$意義下的原根，那么$1,2,...,P?1$中每個數都能寫成$x^t(0<=t<=P?2)$的形式，把0的特殊情況處理掉，然后我們就可以把乘法轉換成加法，也就是求$i?j=x\mathrm{mod}(P?1)(0<=i<=P?2,1<=j<=N)$每個$x$的解。

到這一步以后我就卡了很久，因為不知道怎么進一步去求了，雖然和$P?1$公因數相同的數會產生一樣的循環節，但是如果$N$不是$P?1$的倍數就無從下手了，因此我感覺自己還有什么性質沒有發現，然后我就不停的打表觀察。

直到我打表直接求出每個$x$的解的時候，我發現無論$N$怎么取，和$P?1$公因數相同的數的解的個數是一樣的，那么我們就可以枚舉$P?1$的每個因子，然后暴力枚舉每個$i$有多少個解，這樣時間復雜度是$O(P\sqrt{P})$，理論上是可以過的，雖然有點不太優秀。然后我就開始思考其背后的原因。

假設$N$是不斷增大的，$N$每增大$1$，這個結果都能保持不變，也就是說對于任意一個$j$,$i?j(0<=i<=P?2)$都能滿足和$P?1$公因數相同的數出現個數一樣(固定$j$，變化$i$)。

然后我們嘗試理解一組數乘完某個數以后仍然在同一組。顯然$j=1$時，即$0,1,2,...,P?2$種顯然滿足這個條件，那么我們把相同公因數的數放到同一組里面去，考慮一組里面的數同時乘上同一個數會怎么樣。結論是要么這個公因數保持不變，要么翻倍。不妨設某組數和$P?1$的最大公因數為$d$，這一組里面某個數為$ad$,然后乘以一個$k$。那么$gcd(kad\mathrm{mod}P?1,P?1)$就為$gcd(kad,P?1)$，這就是輾轉相除法，因為$gcd(a,P?1)=1$，所以$gcd(kad,P?1)=gcd(kd,P?1)$。特別地，如果$gcd(kad,P?1)=P?1$，結果就是0。

然后我們還要弄明白一點，就是同一組數乘完某個數以后雖然仍然在同一組，那么是不是那一組里面每個數的個數都是相同的，答案是肯定的。首先證明，對于全部與$P?1$互素的數，如果同時乘上一個與$P?1$互素的數，那么乘完以后還是原來的那些數，只是順序發生了改變。設原來兩個數為$x,y(x\ne y)，那么kx?ky=k(x?y)\ne 0(modP?1)$，也就是兩兩不同。然后我們先考慮$gcd(x,P?1)$為$1$的組，并且乘的$k$是$P?1$的因子的情況，然后就可以根據剛剛證明的性質把全部情況轉換成這種情況，也就是先要提最大公因數，轉換成互質，然后把乘的數拆成兩部分，一部分互質，一部分是因子，前一部分不會改變數的值，只會改變順序，后一部分就是特殊的情況。現在來考慮這種特殊情況，再次強調$k|P?1$,那么我們令$t=\frac{P?1}{k}$，如果兩個數$x,y$(其中$gcd(xy,P?1)=1)$滿足$kx=ky(modP?1)$，那么有$x=y(modt)$，這個問題就能轉化成，所有與$P?1$互素的數在模$t$意義下是不是每個數出現次數相同。想證明這個可以用到前面求歐拉函數的方法，我們不妨考慮和$t$互素的每個數對應多少個和$P?1$互素的個數。因為$k=\frac{P?1}{t}$，然后我們就可以把$k$拆成若干個質數相乘，每乘一個質數就用一次求歐拉函數的方法（同一列都是對應同一個與$t$互素的數），就能知道每個與$t$互素的數在每次操作以后對應的數的個數相同，這樣就能證明這個性質了。

在具體實現的時候，我采用了先枚舉每個$P?1$的因子$i$，然后遍歷整個周期$period=\frac{P?1}{i}$，這樣的復雜度是小于$\sqrt{p}logp$的（但是要先預處理出每個數和$P?1$的最大公因數，否則要多一個log）。因為這是全部因子之和，而我們知道因子和的計算方式為$(1+p_1+p_1^2+...+p_1^{c_1})(1+p_2+p_2^2+...+p_2^{c_2})...(1+p_n+p_n^2+...+p_n^{c_n})={\prod }_{i=1}^n\frac{p_i^{c_i+1}?1}{p_i?1}<={\prod }_{i=1}^n\frac{p_i^{c_i+1}}{p_i?1}={\prod }_{i=1}^n{p}_i^{c_i}{\prod }_{i=1}^n\frac{p_i}{p_i?1}<=(P?1){\prod }_{i=1}^n\frac{i+1}{i}=(n+1)(P?1)<=(P?1)log(P?1)$

然后我們可以找出一個原根，算出每個數是原根的幾次方，然后就能得到每個數出現的次數。然后用ntt卷積一次即可。

#include<bits/stdc++.h>
#define rep(i,x,y) for(int i=x;i<=y;i++)
#define dwn(i,x,y) for(int i=x;i>=y;i--)
#define ll long long
using namespace std;
template<typename T>inline void qr(T &x){x=0;int f=0;char s=getchar();while(!isdigit(s))f|=s=='-',s=getchar();while(isdigit(s))x=x*10+s-48,s=getchar();x=f?-x:x;
}
int cc=0,buf[31];
template<typename T>inline void qw(T x){if(x<0)putchar('-'),x=-x;do{buf[++cc]=int(x%10);x/=10;}while(x);while(cc)putchar(buf[cc--]+'0');
}
const int N=1e5+10,mod=998244353;
namespace NTT{const int G=3,Gi=332748118;int n,m;int lim=1,cnt;int pos[N<<2];int a[N<<2],b[N<<2];int power(int x,int y){int ret=1;while(y){if(y&1)ret=1ll*ret*x%mod;x=1ll*x*x%mod;y>>=1;}return ret;}void ntt(int *A,int type){rep(i,0,lim-1)if(i<pos[i])swap(A[i],A[pos[i]]);for(int mid=1;mid<lim;mid<<=1){int Wn=power(type==1?G:Gi,(mod-1)/(mid<<1));for(int R=mid<<1,j=0;j<lim;j+=R){int w=1;for(int k=0;k<mid;k++,w=1ll*w*Wn%mod){int x=A[j+k],y=1ll*w*A[j+mid+k]%mod;A[j+k]=(x+y)%mod;A[j+k+mid]=(x-y+mod)%mod;}}}}void work(){while(lim<=n+m)lim<<=1,cnt++;rep(i,0,lim-1)pos[i]=(pos[i>>1]>>1)|((i&1)<<(cnt-1));ntt(a,1),ntt(b,1);rep(i,0,lim)a[i]=1ll*a[i]*b[i]%mod;ntt(a,-1);int inv=power(lim,mod-2);rep(i,0,n+m)a[i]=1ll*a[i]*inv%mod;}
}
int P,n;
int cnt,p[N],phi[N];bool v[N];
int sum[N],s[N];
int pre_gcd[N];
int pos[N];
int num[N];
int gcd(int x,int y){return !y?x:gcd(y,x%y);}
void solve(){qr(P),qr(n);rep(i,2,100000){if(!v[i])p[++cnt]=i,phi[i]=i-1;for(int j=1;j<=cnt&&i*p[j]<=100000;j++){v[i*p[j]]=1;if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}else phi[i*p[j]]=phi[i]*(p[j]-1);}}rep(i,2,P-1){int now=1;bool bk=1;rep(j,1,P-1){now=1ll*now*i%P;if(now==1&&j!=P-1){bk=0;break;}pos[now]=j;}if(bk)break;}rep(i,1,P-2){if(gcd(i,P-1)==1)s[i]++;s[i]+=s[i-1];}rep(i,0,P-1)pre_gcd[i]=gcd(i,P-1);sum[P-1]=n%mod;//gcd(x,P-1)=irep(i,1,P-2)if((P-1)%i==0){int period=(P-1)/i;int group_siz=phi[(P-1)/i];rep(j,1,period){if(j==period){(sum[P-1]+=1ll*group_siz*(n/period)%mod)%=mod;break;}int now=pre_gcd[i*j];int p1=n/period%mod,p2=n%period;int group_siz2=phi[(P-1)/now];(sum[now]+=1ll*(group_siz/group_siz2)*(p1+(j<=p2))%mod)%=mod;}}num[0]=n%mod;rep(i,1,P-1)num[i]=sum[gcd(pos[i],P-1)];NTT::n=NTT::m=P-1;rep(i,0,P-1)NTT::a[i]=NTT::b[i]=num[i];NTT::work();int ans=0;rep(i,0,2*(P-1))(ans+=1ll*NTT::a[i]*num[i%P]%mod)%=mod;qw(ans);
}
int main(){int tt;tt=1;while(tt--)solve();return 0;
}