【生成函数+容斥原理+NTT】HDU6036 Division Game

【题目】
$k$堆石子，每堆石子初始数量均为$n$，编号$0$~$k-1$，第$i$次操作对第$(i-1)\%k$堆石子操作，可以从该堆石子中拿走若干石子(至少要拿走一颗)，要求拿走后这堆石子的个数是拿走前这堆石子个数的一个约数。当某堆石子被取走若干石子后变成$1$时结束操作。问最终操作结束于第$i$堆的方案数 。答案对$985661441$取模。
其中$n={\prod }_{i=1}^m{p}_i^{e_i}，m,k\le 10,2\le p_i\le 10^9,\sum e_i\le 10^5,1\le e_i$
原题地址

【解题思路】
这个题很难啊qwq，反正我是没有想到的。那么下面引用官方题解

显然每个石子堆最多做$w={\sum }_{i=1}^m{e}_i$次操作。
现在我们定义$f(x)$为某堆石子做$x$次操作后恰好变为$1$个的生成函数。那么我们可以得到每个数字做少于$x$次操作不变为$1$的生成函数也是$f(x)$。（这两个东西显然是一一对应的）

现在统计结束于第$i$堆石子的情况，我们枚举它是第几次操作时结束的，不妨设为$x$，那么对应的方案数就是$f(x+1{)}^{i-1}f(x{)}^{k-i+1}$，这是因为编号$\le i$的石子堆在$x$次操作后仍然没有变为$1$，编号$\ge i$的石子堆在$x-1$次操作后没有变为$1$，只有编号为$i$的石子堆恰好变为$1$。
那么如果我们知道$f(x)$，则可以用$O(wk)$的时间解决问题。

考虑$f(x)$的组成，不妨设某种方案第$j$次操作$e_i$减少了$d(i,j)$，我们知道对于每个$e_i$有${\sum }_{j=1}^{x}d(i,j)=e_i$，且对于每个$j$有${\sum }_{i=1}^{m}d(i,j)>0$。
那么$f(x)$就是分配$d(i,j)(1\le i\le m,1\le j\le x)$使得它满足上面两个条件的方案数。

如果只考虑第一个条件对应的方案数$g(x)$，我们可以规约到一个经典的组合问题，求一个方程非负整数解的个数。
$$g(x)=\prod _{i=1}^m\left(\genfrac{}{}{0ex}{}{e_i+x-1}{x-1}\right)$$

这个公式可以通过对$x-1$个隔板，$e_i$个物品的排列数量得到，隔板把区间分成了$x$块，每块区间中物品的数量对应一个变量，排列的数量等同于隔板选位置的方法数。

观察到如果某些$j$与第二个条件产生了矛盾，与之相关的$d(i,j)$都会是零。利用容斥原理可以得到
$$f(x)=\sum _{y=0}^x(-1{)}^{x-y}(\genfrac{}{}{0ex}{}{x}{y})g(y)$$
可以发现：
$$\frac{f(x)}{x!}=\sum _{y=0}^x\frac{(-1{)}^{x-y}}{(x-y)!}\cdot \frac{g(y)}{y!}$$
用$拆系数FFT$或$NTT$即可。$(985661441=235\times 2^{22}+1)$

总的时间复杂度是$O(w\log n+wk)$。

【参考代码】

#include<bits/stdc++.h>
using namespace std;

typedef long long LL;
const int N=2e5+10;
const int mod=985661441,g=985661438;
int n,m,k,L,pris;
LL fac[N],ifac[N],f[N],e[N],rev[N<<1],a[N<<1],b[N<<1];

LL qpow(LL x,LL y){LL ret=1;for(;y;y>>=1,x=x*x%mod) if(y&1) ret=ret*x%mod;return ret;}
LL C(int n,int m){return fac[n]*ifac[m]%mod*ifac[n-m]%mod;}
void up(LL &x,LL y){x+=y;if(x>=mod)x-=mod;}

void totinit()
{
	fac[0]=fac[1]=1;
	for(int i=2;i<N;++i) fac[i]=fac[i-1]*i%mod;
	ifac[N-1]=qpow(fac[N-1],mod-2);
	for(int i=N-1;i;--i) ifac[i-1]=ifac[i]*i%mod;
}

void ntt(LL *a,int n,int f)
{
	for(int i=0;i<n;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=1;i<n;i<<=1)
	{
		LL wn=qpow(g,f==1?(mod-1)/i/2:mod-1-(mod-1)/i/2);
		for(int j=0;j<n;j+=(i<<1))
		{
			LL w=1;
			for(int k=0;k<i;++k,w=w*wn%mod)
			{
				LL x=a[j+k],y=w*a[i+j+k]%mod;
				a[j+k]=(x+y)%mod;a[i+j+k]=(x-y+mod)%mod;
			}
		}
	}
	LL inv=qpow(n,mod-2);
	if(f==-1) for(int i=0;i<n;++i) a[i]=a[i]*inv%mod;
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("HDU6036.in","r",stdin);
	freopen("HDU6036.out","w",stdout);
#endif
	totinit();int cas=0;
	while(scanf("%d%d",&pris,&k)!=EOF)
	{
		++cas;printf("Case #%d: ",cas);n=0;
		for(int i=1,x;i<=pris;++i) scanf("%d%d",&x,&e[i]),n+=e[i];
		for(L=0,m=1;m<=n<<1;m<<=1,++L);
		for(int i=0;i<m;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)),a[i]=b[i]=0;
		for(int i=0;i<=n;++i) a[i]=(i&1)?mod-ifac[i]:ifac[i];
		for(int i=1;i<=n;++i)
		{
			b[i]=ifac[i];
			for(int j=1;j<=pris;++j) b[i]=b[i]*C(e[j]+i-1,i-1)%mod;
		}
		ntt(a,m,1);ntt(b,m,1);
		for(int i=0;i<m;++i) a[i]=a[i]*b[i]%mod;
		ntt(a,m,-1);
		for(int i=0;i<=n;++i) f[i]=a[i]*fac[i]%mod; f[n+1]=0;
		for(int i=1;i<=k;++i)
		{
			LL ans=0;
			for(int j=1;j<=n;++j) up(ans,qpow(f[j],k-i+1)*qpow(f[j+1],i-1)%mod);
			if(i<k) printf("%lld ",ans); else printf("%lld\n",ans);
		}
	}

	return 0;
}

【总结】
论如何构造生成函数计算。