LibreOJ #3120.「CTS2019 | CTSC2019」珍珠 生成函数

题意

有$n$个范围在$[1,D]$之间的整数均匀随机变量，问能够选出$m$对变量，满足每个变量至多只出现在一对变量中，且同一对变量的值相同的概率。
$n,m\le 10^9,D\le 10^5$

分析

显然一个数值出现偶数次的$EGF$为$\frac{e^x+e^{-x}}{2}$，奇数次则为$\frac{e^x-e^{-x}}{2}$
那么答案就是$$n!\sum _{k=0}^{n-2m}(\frac{e^x+e^{-x}}{2}+y\frac{e^x-e^{-x}}{2}{)}^D[x^n][y^k]$$
$$=\frac{n!}{2^D}\sum _{k=0}^{n-2m}[e^x(1+y)+e^{-x}(1-y){]}^D[x^n][y^k]$$
$$=\frac{n!}{2^D}\sum _{k=0}^{n-2m}\sum _{i=0}^D{C}_D^i{e}^{(2i-D)x}(1+y{)}^i(1-y{)}^{D-i}[x^n][y^k]$$
$$=\frac{1}{2^D}\sum _{i=0}^D{C}_D^i(2i-D{)}^n\sum _{k=0}^{n-2m}(1+y{)}^i(1-y{)}^{D-i}[y^k]$$对于后面那部分有$$\sum _{k=0}^{n-2m}(1+y{)}^i(1-y{)}^{D-i}[y^k]$$
$$=(1+y{)}^i(1-y{)}^{D-i}(1+y+y^2+...)[y^{n-2m}]$$
$$=(1+y{)}^i(1-y{)}^{D-i}\frac{1}{1-y}[y^{n-2m}]$$
当$D=i$时可以暴力求，当$D≠i$时上式则等于$$(1+y{)}^i(1-y{)}^{D-i-1}[y^{n-2m}]$$
二项式展开之后就变成了一个卷积形式，用$FFT$直接做就好。
时间复杂度$O(DlogD)$。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>

typedef long long LL;

const int N=400005;
const int MOD=998244353;

int n,m,D,jc[N],ny[N],a[N],b[N],rev[N],L;

int ksm(int x,int y)
{
	int ans=1;
	while (y)
	{
		if (y&1) ans=(LL)ans*x%MOD;
		x=(LL)x*x%MOD;y>>=1;
	}
	return ans;
}

int C(int n,int m)
{
	if (n<m) return 0;
	else return (LL)jc[n]*ny[m]%MOD*ny[n-m]%MOD;
}

void NTT(int *a,int f)
{
	for (int i=0;i<L;i++) if (i<rev[i]) std::swap(a[i],a[rev[i]]);
	for (int i=1;i<L;i<<=1)
	{
		int wn=ksm(3,f==1?(MOD-1)/i/2:MOD-1-(MOD-1)/i/2);
		for (int j=0;j<L;j+=(i<<1))
		{
			int w=1;
			for (int k=0;k<i;k++)
			{
				int u=a[j+k],v=(LL)a[j+k+i]*w%MOD;
				a[j+k]=(u+v)%MOD;a[j+k+i]=(u+MOD-v)%MOD;
				w=(LL)w*wn%MOD;
			}
		}
	}
	if (f==-1) for (int i=0,ny=ksm(L,MOD-2);i<L;i++) a[i]=(LL)a[i]*ny%MOD;
}

int main()
{
	scanf("%d%d%d",&D,&n,&m);
	jc[0]=jc[1]=ny[0]=ny[1]=1;
	for (int i=2;i<=D;i++) jc[i]=(LL)jc[i-1]*i%MOD,ny[i]=(LL)(MOD-MOD/i)*ny[MOD%i]%MOD;
	for (int i=2;i<=D;i++) ny[i]=(LL)ny[i-1]*ny[i]%MOD;
	int k=n-m*2,d=D-1;
	for (int i=0;i<=D;i++)
	{
		if (i<=k&&k-i<=D)
		{
			a[i]=(LL)ny[i]*ny[k-i]%MOD;
			if ((k-i)&1) a[i]=MOD-a[i];
		}
		if (d-k-i>=0) b[i]=(LL)ny[i]*ny[d-k-i]%MOD;
 	}
 	int lg=0;
 	for (L=1;L<=D*2;L<<=1,lg++);
 	for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
 	NTT(a,1);NTT(b,1);
 	for (int i=0;i<L;i++) a[i]=(LL)a[i]*b[i]%MOD;
 	NTT(a,-1);
 	int ans=0;
 	for (int i=0;i<D;i++) (ans+=(LL)C(D,i)*ksm(i*2-D,n)%MOD*jc[i]%MOD*jc[d-i]%MOD*a[i]%MOD)%=MOD;
 	for (int i=0;i<=std::min(D,k);i++) (ans+=(LL)C(D,i)*ksm(D,n)%MOD)%=MOD;
 	ans=(LL)ans*ksm(ksm(2,D),MOD-2)%MOD;
 	ans+=ans<0?MOD:0;
 	printf("%d\n",ans);
 	return 0;
 }