以前基本没有了解原根相关的一块内容,最近正好碰到了这个题,于是写一篇博客记录一下。
这道题的总体思路就是比较明显,就是先算出 a p = x a^p=x ap=x对于每个 x x x的解的个数,然后NTT算一下即可。
先来讲一下怎么求欧拉函数 ϕ ( x ) \phi (x) ϕ(x),即1到 x − 1 x-1 x−1中与 x x x互素的数的个数,我们一般在线性筛里面通过递推来求出这个函数。也就是我们要找出 ϕ ( i ∗ p [ j ] ) \phi(i*p[j]) ϕ(i∗p[j])和 ϕ ( i ) , ϕ ( p [ j ] ) , i , p [ j ] \phi (i),\phi(p[j]),i,p[j] ϕ(i),ϕ(p[j]),i,p[j]之间的关系。我们把与 ϕ ( i ) \phi(i) ϕ(i)所代表的集合写在一行里面,然后第二行写上前一行的数加上 i i i,一直写 p [ j ] p[j] 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 m o d p [ j ] ≠ 0 i \mod p[j] \ne 0 imodp[j]=0,也就是 ϕ ( i ) \phi(i) ϕ(i)所表示的集合中存在 p [ j ] p[j] p[j]的倍数。然后你可以发现,每一列里面恰好就有一个是 p [ j ] p[j] p[j]的倍数,我们不妨设某一列第一个数为 x x x,且 y = x m o d p [ j ] y=x\mod p[j] y=xmodp[j],那么因为 p [ j ] p[j] p[j]为质数,所以 y m o d p [ j ] , y + i m o d p [ j ] , y + 2 i m o d p [ j ] , . . . , y + ( p [ j ] − 1 ) i m o d p [ j ] y\mod p[j],y+i\mod p[j],y+2i\mod p[j],...,y+(p[j]-1)i\mod p[j] ymodp[j],y+imodp[j],y+2imodp[j],...,y+(p[j]−1)imodp[j]中 0 , 1 , . . . , p [ j ] − 1 0,1,...,p[j]-1 0,1,...,p[j]−1恰好每个一个,因此可以知道 p h i [ i ∗ p [ j ] ] = p h i [ i ] ∗ ( p [ j ] − 1 ) ( i m o d p [ j ] ≠ 0 ) phi[i*p[j]]=phi[i]*(p[j]-1)(i \mod p[j] \ne 0) phi[i∗p[j]]=phi[i]∗(p[j]−1)(imodp[j]=0)。
在例2中,因为 i m o d p [ j ] = 0 i\mod p[j] =0 imodp[j]=0,因此 ϕ ( i ) \phi(i) ϕ(i)所表示的集合里面没有 p [ j ] p[j] p[j]的倍数,因此第一行里面没有被删掉的,因为相邻两行之间的差都是 p [ j ] p[j] p[j]的倍数,所以不会有任何的数为 p [ j ] p[j] p[j]的倍数。
然后就是原根部分的知识。
如果 x x x是模 P P P意义下的原根,那么 1 , 2 , . . . , P − 1 1,2,...,P-1 1,2,...,P−1中每个数都能写成 x t ( 0 < = t < = P − 2 ) x^t(0<=t<=P-2) xt(0<=t<=P−2)的形式,把0的特殊情况处理掉,然后我们就可以把乘法转换成加法,也就是求 i ∗ j = x m o d ( P − 1 ) ( 0 < = i < = P − 2 , 1 < = j < = N ) i*j=x\mod (P-1)(0<=i<=P-2,1<=j<=N) i∗j=xmod(P−1)(0<=i<=P−2,1<=j<=N)每个 x x x的解。
到这一步以后我就卡了很久,因为不知道怎么进一步去求了,虽然和 P − 1 P-1 P−1公因数相同的数会产生一样的循环节,但是如果 N N N不是 P − 1 P-1 P−1的倍数就无从下手了,因此我感觉自己还有什么性质没有发现,然后我就不停的打表观察。
直到我打表直接求出每个 x x x的解的时候,我发现无论 N N N怎么取,和 P − 1 P-1 P−1公因数相同的数的解的个数是一样的,那么我们就可以枚举 P − 1 P-1 P−1的每个因子,然后暴力枚举每个 i i i有多少个解,这样时间复杂度是 O ( P P ) O(P\sqrt P) O(PP),理论上是可以过的,虽然有点不太优秀。然后我就开始思考其背后的原因。
假设 N N N是不断增大的, N N N每增大 1 1 1,这个结果都能保持不变,也就是说对于任意一个 j j j, i ∗ j ( 0 < = i < = P − 2 ) i*j(0<=i<=P-2) i∗j(0<=i<=P−2)都能满足和 P − 1 P-1 P−1公因数相同的数出现个数一样(固定 j j j,变化 i i i)。
然后我们尝试理解一组数乘完某个数以后仍然在同一组。显然 j = 1 j=1 j=1时,即 0 , 1 , 2 , . . . , P − 2 0,1,2,...,P-2 0,1,2,...,P−2种显然满足这个条件,那么我们把相同公因数的数放到同一组里面去,考虑一组里面的数同时乘上同一个数会怎么样。结论是要么这个公因数保持不变,要么翻倍。不妨设某组数和 P − 1 P-1 P−1的最大公因数为 d d d,这一组里面某个数为 a d ad ad,然后乘以一个 k k k。那么 g c d ( k a d m o d P − 1 , P − 1 ) gcd(kad\mod P-1,P-1) gcd(kadmodP−1,P−1)就为 g c d ( k a d , P − 1 ) gcd(kad,P-1) gcd(kad,P−1),这就是辗转相除法,因为 g c d ( a , P − 1 ) = 1 gcd(a,P-1)=1 gcd(a,P−1)=1,所以 g c d ( k a d , P − 1 ) = g c d ( k d , P − 1 ) gcd(kad,P-1)=gcd(kd,P-1) gcd(kad,P−1)=gcd(kd,P−1)。特别地,如果 g c d ( k a d , P − 1 ) = P − 1 gcd(kad,P-1)=P-1 gcd(kad,P−1)=P−1,结果就是0。
然后我们还要弄明白一点,就是同一组数乘完某个数以后虽然仍然在同一组,那么是不是那一组里面每个数的个数都是相同的,答案是肯定的。首先证明,对于全部与 P − 1 P-1 P−1互素的数,如果同时乘上一个与 P − 1 P-1 P−1互素的数,那么乘完以后还是原来的那些数,只是顺序发生了改变。设原来两个数为 x , y ( x ≠ y ) ,那么 k x − k y = k ( x − y ) ≠ 0 ( m o d P − 1 ) x,y(x\ne y),那么kx-ky=k(x-y) \ne 0(mod\ P-1) x,y(x=y),那么kx−ky=k(x−y)=0(mod P−1),也就是两两不同。然后我们先考虑 g c d ( x , P − 1 ) gcd(x,P-1) gcd(x,P−1)为 1 1 1的组,并且乘的 k k k是 P − 1 P-1 P−1的因子的情况,然后就可以根据刚刚证明的性质把全部情况转换成这种情况,也就是先要提最大公因数,转换成互质,然后把乘的数拆成两部分,一部分互质,一部分是因子,前一部分不会改变数的值,只会改变顺序,后一部分就是特殊的情况。现在来考虑这种特殊情况,再次强调 k ∣ P − 1 k|P-1 k∣P−1,那么我们令 t = P − 1 k t=\frac{P-1}{k} t=kP−1,如果两个数 x , y x,y x,y(其中 g c d ( x y , P − 1 ) = 1 ) gcd(xy,P-1)=1) gcd(xy,P−1)=1)满足 k x = k y ( m o d P − 1 ) kx=ky(mod\ P-1) kx=ky(mod P−1),那么有 x = y ( m o d t ) x=y(mod\ t) x=y(mod t),这个问题就能转化成,所有与 P − 1 P-1 P−1互素的数在模 t t t意义下是不是每个数出现次数相同。想证明这个可以用到前面求欧拉函数的方法,我们不妨考虑和 t t t互素的每个数对应多少个和 P − 1 P-1 P−1互素的个数。因为 k = P − 1 t k=\frac{P-1}{t} k=tP−1,然后我们就可以把 k k k拆成若干个质数相乘,每乘一个质数就用一次求欧拉函数的方法(同一列都是对应同一个与 t t t互素的数),就能知道每个与 t t t互素的数在每次操作以后对应的数的个数相同,这样就能证明这个性质了。
在具体实现的时候,我采用了先枚举每个 P − 1 P-1 P−1的因子 i i i,然后遍历整个周期 p e r i o d = P − 1 i period=\frac{P-1}{i} period=iP−1,这样的复杂度是小于 p l o g p \sqrt plogp plogp的(但是要先预处理出每个数和 P − 1 P-1 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 ) = ∏ i = 1 n p i c i + 1 − 1 p i − 1 < = ∏ i = 1 n p i c i + 1 p i − 1 = ∏ i = 1 n p i c i ∏ i = 1 n p i p i − 1 < = ( P − 1 ) ∏ i = 1 n i + 1 i = ( n + 1 ) ( P − 1 ) < = ( P − 1 ) l o g ( P − 1 ) (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}^np_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) (1+p1+p12+...+p1c1)(1+p2+p22+...+p2c2)...(1+pn+pn2+...+pncn)=∏i=1npi−1pici+1−1<=∏i=1npi−1pici+1=∏i=1npici∏i=1npi−1pi<=(P−1)∏i=1nii+1=(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;
}