题意:有$n$堆石子,每堆石子数量相同,以质因数分解给出,不停地从$1$到$n$依次拿石子,使得取完后石子个数为原来的因数(不能不取),当一堆只剩$1$个时结束,问在每堆石子结束的方案数
记石子个数为$\prod\limits_{i=1}^mp_i^{e_i}$,拿石头就是选一些$e_i$把它们减小
令一堆石子取$x$次取完的方案数为$f_x$,取$x-1$次还未取完的方案数为$g_x$,那么$g_x=f_x$,因为对于每个取$x-1$次未取完的方案,把它取完就对应一个$x$次取完的方案,对于每个$x$次取完的方案,它的前$x-1$次操作一定不完全相同,所以两种方案数相同
记$s=\sum e_i$,那么取$j$次取完并在$i$结束的方案数为$g_{j+1}^{i-1}f_jg_j^{n-i}=f_{j+1}^{i-1}f_j^{n-i+1}$,所以答案为$ans_i=\sum\limits_{j=1}^sf_{j+1}^{i-1}f_j^{n-i+1}$
现在我们来算$f$,直接算将$e_i$分成$x$份的方案数就是$h_x=\prod\limits_{i=1}^m\binom{e_i+x-1}{x-1}$,但$h_x\ne f_x$,因为$h_x$中可能包含那些某一次一个石子都没取的方案,所以要容斥
取$x$次中,至少有$y$次一个石子都没取的方案数为$\binom xyh_{y-x}$,那么$f_x=\sum\limits_{y=0}^x(-1)^y\binom xyh_{x-y}$,也就是说求出$h$后我们就可以$O(s\log s)$求$f$了
原题中$h$可以暴力求,但毒瘤代爷稍微加强了一下数据,注意到因为$\sum e_i=s$,所以不同的$e_i$最多会有$O(\sqrt s)$个,所以我们对相同的$e_i$一起处理,就可以在$O(s\sqrt s\log s)$的时间内递推求得$h$
总时间复杂度$O(s\sqrt s\log s)$
#include<stdio.h> #include<algorithm> using namespace std; typedef long long ll; const int mod=998244353; int mul(int a,int b){return(ll)a*b%mod;} int ad(int a,int b){return(a+b)%mod;} int de(int a,int b){return(a-b)%mod;} void inc(int&a,int b){(a+=b)%=mod;} int pow(int a,int b){ int s=1; while(b){ if(b&1)s=mul(s,a); a=mul(a,a); b>>=1; } return s; } int rev[1048576],N,iN; void pre(int n){ int i,k=0; for(N=1,k=0;N<=n;N<<=1)k++; for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1)); iN=pow(N,mod-2); } void ntt(int*a,int on){ int i,j,k,t,w,wn; for(i=0;i<N;i++){ if(i<rev[i])swap(a[i],a[rev[i]]); } for(i=2;i<=N;i<<=1){ wn=pow(3,on==1?(mod-1)/i:mod-1-(mod-1)/i); for(j=0;j<N;j+=i){ w=1; for(k=0;k<i>>1;k++){ t=mul(a[i/2+j+k],w); a[i/2+j+k]=de(a[j+k],t); a[j+k]=ad(a[j+k],t); w=mul(w,wn); } } } if(on==-1){ for(i=0;i<N;i++)a[i]=mul(a[i],iN); } } int e[10010],inv[300010],fac[300010],rfac[300010],h[300010],f[300010],c[300010],pe[300010],pc[300010]; int C(int n,int k){return mul(fac[n],mul(rfac[k],rfac[n-k]));} int a[1048576],b[1048576]; int main(){ int m,n,s,i,j,t,M; scanf("%d%d",&m,&n); s=0; for(i=1;i<=m;i++){ scanf("%d%d",&t,e+i); s+=e[i]; c[e[i]]++; } M=0; for(i=1;i<=s;i++){ if(c[i]){ M++; pe[M]=i; pc[M]=c[i]; } } fac[0]=1; for(i=1;i<=s;i++)fac[i]=mul(fac[i-1],i); rfac[s]=pow(fac[s],mod-2); for(i=s;i>0;i--)rfac[i-1]=mul(rfac[i],i); inv[1]=1; for(i=2;i<=s;i++)inv[i]=-mul(mod/i,inv[mod%i]); h[1]=1; for(i=1;i<s;i++){ h[i+1]=h[i]; for(j=1;j<=M;j++)h[i+1]=mul(h[i+1],pow(mul(pe[j]+i,inv[i]),pc[j])); } pre(s<<1); for(i=0;i<=s;i++){ a[i]=(i&1?-1:1)*rfac[i]; b[i]=mul(h[i],rfac[i]); } ntt(a,1); ntt(b,1); for(i=0;i<N;i++)a[i]=mul(a[i],b[i]); ntt(a,-1); for(i=0;i<=s;i++)f[i]=mul(a[i],fac[i]); for(i=1;i<=n;i++){ t=0; for(j=1;j<=s;j++)inc(t,mul(pow(f[j+1],i-1),pow(f[j],n-i+1))); inc(t,mod); printf("%d\n",t); } }