(资料图片仅供参考)
题目链接
令\(c_i\)表示第i种颜色的珍珠的数量，显然我们最多能装的瓶数是\(\sum \lfloor \frac {c_i}2 \rfloor\)。也就是说，\(c_i\)为奇数的\(i\)的数量不能太多，这个数量要在某个范围内才行，且这个范围是很好求的。考虑对每个n求出\(f_n\)表示刚好有n个\(c_i\)为奇数的方案数。
如何指定某些颜色一定选奇数个，并把所有不同的颜色排列组合成一行n个？考虑使用指数型生成函数(EGF)。复习一下EGF相乘的组合意义：\(A(x)=\sum_{i\ge 0}a_i\frac {x^i}{i!},B(x)=\sum_{i\ge 0}b_i\frac {x^i}{i!}\)，其中\(a_i\)表示i个元素的集合做A操作的方案数，\(b_i\)表示i个元素的集合做B操作的方案数。若\(C(x)=A(x)\cdot B(x)=\sum_{i\ge 0}c_i\frac {x^i}{i!}\)，则\(c_i\)表示i个元素的集合分成两个子集，其中第一个做A操作，第二个做B操作的总方案数。那么可以把i个元素的集合"有奇数个"也看成一种操作，\(0,1,2,3\cdots\)个元素的集合对应操作数是\(0,1,0,1\cdots\)。由于存在如下恒等式\(e^x=\sum_{i\ge 0}\frac{x^i}{i!}\)，所以\(e^x\)的系数序列是\(1,1,1,1\cdots\)，\(e^{-x}\)的系数序列是\(1,-1,1,-1\cdots\)，那么\(0,1,0,1\cdots\)其实就是\(\frac{e^x-e^{-x}}2\)对应的系数序列。i个元素有偶数个的方案数"系数序列"也一样，可以用\(\frac{e^x+e^{-x}}2\)表示。
那么\(f_i\)就能表示成\(([n](\frac{e^x-e^{-x}}2)^i(\frac{e^x+e^{-x}}2)^{D-i})\cdot n!\cdot \binom di\)。其中\([n]\)表示这个函数的\(n\)次项系数，乘\(n!\)是因为这个函数是EGF。但是这样似乎不太能继续推了，考虑重新定义\(f\)。令\(g_i\)表示原来\(f_i\)的定义。对于\(f_i\)：钦定任意i个颜色必须选奇数个，其他颜色随便，并把方案数贡献到\(f_i\)。这样，对于\(i现在的\(f_i=([n](\frac{e^x-e^{-x}}2)^i(e^x)^{D-i})\cdot n!\cdot \binom di\)。直接二项式定理暴力展开：
\[\begin{align}f_i&=(\frac 12)^in!\binom di \cdot\sum_{j=0}^i [n]e^{D+2j-2i}(-1)^{i-j}\binom ij\\&=(\frac 12)^in!\binom di \cdot\sum_{j=0}^i \frac{(D+2j-2i)^n}{n!} (-1)^{i-j}\binom ij\\&=(\frac 12)^i\binom di \cdot\sum_{j=0}^i (D+2j-2i)^n(-1)^{i-j}\binom ij\\\end{align}\]显然，这玩意儿可以写成一个以\(j\)为下标的序列和一个\(i-j\)为下标的序列卷积得到\(f_i\)的形式，所以可以\(O(DlogD)\)求出\(f\)。
总复杂度\(O(DlogD)\)。
点击查看代码
#include #define rep(i,n) for(int i=0;i#define fi first#define se second#define mpr make_pair#define pb push_backvoid fileio(){ #ifdef LGS freopen("in.txt","r",stdin); freopen("out.txt","w",stdout); #endif}void termin(){ #ifdef LGS std::cout<<"\n\nEXECUTION TERMINATED"; #endif exit(0);}using namespace std;const LL MOD=998244353;LL qpow(LL x,LL a){LL res=x,ret=1;while(a>0){if((a&1)==1) ret=ret*res%MOD;a>>=1;res=res*res%MOD;}return ret;}namespace NTT{ vector rev; void ntt(vector &a,LL G) { LL nn=a.size(),gn,g,x,y;vector tmp=a; rep(i,nn) a[i]=tmp[rev[i]]; for(int len=1;len convolution(vector a,vector b,LL G) { if(a.size()<=23&&b.size()<=23) { vector ret(a.size()+b.size()-1,0); rep(i,a.size()) rep(j,b.size()) (ret[i+j]+=a[i]*b[j])%=MOD; return ret; } LL nn=1,bt=0,sv=a.size()+b.size()-1;while(nn>1]>>1)|((i&1)<<(bt-1)); } ntt(a,G);ntt(b,G); rep(i,nn) (a[i]*=b[i])%=MOD; ntt(a,qpow(G,MOD-2)); while(a.size()>sv) a.pop_back(); LL inv=qpow(nn,MOD-2); rep(i,a.size()) (a[i]*=inv)%=MOD; return a; }}LL d,n,m,fac[100010],inv[100010],f[100010],g[100010];LL CC(LL nn,LL mm){return fac[nn]*inv[mm]%MOD*inv[nn-mm]%MOD;}int main(){ fileio(); fac[0]=1;repn(i,100005) fac[i]=fac[i-1]*i%MOD; rep(i,100005) inv[i]=qpow(fac[i],MOD-2); cin>>d>>n>>m; if(n==1) { puts("0"); termin(); } vector A,B,C; rep(i,d+2) A.pb(inv[i]); rep(i,d+2) { LL val=inv[i]; if(i%2==1) val=(MOD-val)%MOD; (val*=qpow((MOD+d-i-i)%MOD,n))%=MOD; B.pb(val); } C=NTT::convolution(A,B,3); rep(i,d+1) f[i]=C[i]*CC(d,i)%MOD*qpow(qpow(2,i),MOD-2)%MOD*fac[i]%MOD; //cout<=m+m) (ans+=g[i])%=MOD; cout<
