发现可以用中国剩余定理对 M = p α , A i ∈ [ 1 , p α ] , N = N m o d p α M=p^{\alpha},A_i\in [1,p^\alpha],N=N\bmod p^\alpha M=pα,Ai∈[1,pα],N=Nmodpα求出答案,然后再乘起来。对于每个 A i A_i Ai通过中国剩余定理合并出来的答案也是 唯一 的。
设 N = x × p t N=x\times p^t N=x×pt,其中 x ∈ [ 1 , p α − t ] x\in [1,p^{\alpha-t}] x∈[1,pα−t]且与 p p p互质。这里假设 N > 0 N>0 N>0。
显然,此时 ∏ i = 1 n A i \prod_{i=1}^nA_i ∏i=1nAi一定恰好包含了 t t t个 p p p(否则无法得到 N N N),假设 A i A_i Ai包含了 a i a_i ai个 p p p,那么 A i A_i Ai的方案数为 p α − a i − 1 ( p − 1 ) p^{\alpha-a_i-1}(p-1) pα−ai−1(p−1)。
根据简约剩余系的理论,可以证明 x x x取到任意值的方案数是相等的。证明可以考虑 A n A_n An不为 p p p倍数的部分一定包含了所有取值。
因此方案数为 ( n + t − 1 n − 1 ) × p n ( α − 1 ) − t × ( p − 1 ) n × 1 p α − t − 1 ( p − 1 ) = ( n + t − 1 n − 1 ) × p n α − n − α + 1 × ( p − 1 ) n − 1 \binom{n+t-1}{n-1}\times p^{n(\alpha-1)-t}\times (p-1)^n\times \frac{1}{p^{\alpha-t-1}(p-1)}=\binom{n+t-1}{n-1}\times p^{n\alpha-n-\alpha+1}\times (p-1)^{n-1} (n−1n+t−1)×pn(α−1)−t×(p−1)n×pα−t−1(p−1)1=(n−1n+t−1)×pnα−n−α+1×(p−1)n−1
对于 N = 0 N=0 N=0的情况简单容斥即可。注意要算上 x x x的方案数。
upd:这题有一些求逆元的做法好像会寄掉。但是题解的做法没用到逆元。其实应该在大样例里面加一组这种数据的。这个确实没考虑到。
#include<bits/stdc++.h>
#define ll long long
#define pb push_back
#define fi first
#define se second
#define db double
using namespace std;
const int mod=998244353;
ll n,N,M,res=1;
ll fpow(ll x,ll y=mod-2){
if(y<0)return 0;
x%=mod;ll z(1);
for(;y;y>>=1){
if(y&1)z=z*x%mod;
x=x*x%mod;
}return z;
}
ll binom(ll x,ll y){
ll res=1;
for(int i=1;i<=y;i++)res=res*fpow(i)%mod;
for(int i=0;i<y;i++)res=res*((x-i)%mod)%mod;
return res;
}
void add(ll &x,ll y){
x=(x+y)%mod;
}
ll solve(ll p,ll N,int k){
if(N){
int t=0;
while(N%p==0)N/=p,t++;
return binom(n+t-1,t)*fpow(p,k*n-t-n+1-(k-t))%mod*fpow(p-1,n-1)%mod;
}
else{
ll res=fpow(p,k*n);
for(int t=0;t<k;t++){
add(res,-fpow(p,k-t-1)*((p-1)%mod)%mod*binom(n+t-1,t)%mod*fpow(p,k*n-t-n+1-(k-t))%mod*fpow(p-1,n-1)%mod);
}
return res;
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin>>n>>N>>M;
for(ll i=2;i<=M/i;i++){
if(M%i==0){
ll P=1,k=0;
while(M%i==0)M/=i,P*=i,k++;
ll n2=N%P;
res=res*solve(i,n2,k)%mod;
}
}
if(M>1){
res=res*solve(M,N%M,1)%mod;
}
cout<<(res+mod)%mod;
}