题意
给定 $n$ , 求下式的值:
$$ f(n)= \sum_{i=0}^n\sum_{j=0}^i\begin{Bmatrix}i\\ j\end{Bmatrix}\times 2^j\times j!$$
题解
这题比较神仙...
那么我们可以思考如何来求一个比较简单的转移式.
首先我们发现, $f(n)$ 表达式中的第一重和式包含了 $f(n-1)$, 那么我们对 $f$ 的值做差分, 于是我们有 $f(n)-f(n-1)=\sum\limits_{i=0}^n\begin{Bmatrix}i\\ j\end{Bmatrix}\times 2^j\times j!$ .设 $g(n)=f(n)-f(n-1)$.
而我们知道第二类斯特林数的组合意义是在 $i$ 个物品中选 $j$ 个子集的方案数. 那么后面的 $j!$ 的组合意义相当于让划分的子集带标号, $2^j$ 的组合意义相当于每个子集贡献两种状态.
知道了表达式的组合意义, 我们可以尝试用朴素DP来解它. 我们枚举一下最后一个子集中的元素个数为 $i$ (因为子集带标号, 我们相当于枚举最后一个), 那么我们就有递推式:
$$g(n)=[n=0]+\sum_{i=1}^n {n\choose i} g(n-i)\times 2$$
最后乘 $2$ 是因为每个集合会贡献两种状态, $[n=0]$ 是边界条件.
观察这个式子, 发现是二项卷积的形式, 我们果断上EGF来解决. 设 $G(x)$ 是 $g(n)$ 的EGF, $H(x)$ 是数列 $\langle 0,2,2,2,\dots \rangle$ 的EGF, 那么我们有:
$$ G(x)=G(x)H(x)+1 $$
注意到 $H(x)$ 代表的数列的第 $0$ 项是 $0$, 因为 $g(n)$ 的递推式中求和下指标是 $1$.
那么我们移项之后就可以开心地得到:
$$ G(x)=\frac 1 {1-H(x)} $$
NTT爆算一发就可以了
代码实现
记得EGF要除以阶乘...出答案的时候还得乘回来...
#include <bits/stdc++.h> const int G=;
const int DFT=;
const int IDFT=-;
const int MAXN=;
const int MOD=;
const int INV2=(MOD+)>>;
const int PHI=MOD-; typedef std::vector<int> Poly; Poly Sqrt(Poly);
void Read(Poly&);
Poly Inverse(Poly);
Poly Ln(const Poly&);
Poly Exp(const Poly&);
void Print(const Poly&);
void NTT(Poly&,int,int);
Poly Pow(const Poly&,int);
Poly Integral(const Poly&);
Poly Derivative(const Poly&);
Poly operator*(Poly,Poly);
Poly operator/(Poly,Poly);
Poly operator%(Poly,Poly);
Poly operator+(const Poly&,const Poly&);
Poly operator-(const Poly&,const Poly&); int rev[MAXN]; int NTTPre(int);
int Pow(int,int,int); int main(){
int n;
scanf("%d",&n);
Poly h(n+),one(,);
h[]=;
for(int i=;i<=n;i++)
h[i]=1ll*h[i-]*Pow(i,MOD-,MOD)%MOD;
Poly g(Inverse(one-h));
int ans=;
int fact=;
for(int i=;i<=n;i++){
(ans+=1ll*g[i]*fact%MOD)%=MOD;
fact=1ll*fact*(i+)%MOD;
}
printf("%d\n",ans);
return ;
} void Read(Poly& a){
for(auto& i:a)
scanf("%d",&i);
} void Print(const Poly& a){
for(auto i:a)
printf("%d ",i);
puts("");
} Poly Pow(const Poly& a,int k){
Poly log=Ln(a);
for(auto& i:log)
i=1ll*i*k%MOD;
return Exp(log);
} Poly Sqrt(Poly a){
int len=a.size();
if(len==)
return Poly(,int(sqrt(a[])));
else{
Poly b=a;
b.resize((len+)>>);
b=Sqrt(b);
b.resize(len);
Poly inv=Inverse(b);
int bln=NTTPre(inv.size()+a.size());
NTT(a,bln,DFT);
NTT(inv,bln,DFT);
for(int i=;i<bln;i++)
a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
NTT(a,bln,IDFT);
for(int i=;i<len;i++)
b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
return b;
}
} Poly Exp(const Poly& a){
size_t len=;
Poly ans(,),one(,);
while(len<(a.size()<<)){
len<<=;
Poly b=a;
b.resize(len);
ans=ans*(one-Ln(ans)+b);
ans.resize(len);
}
ans.resize(a.size());
return ans;
} Poly Ln(const Poly& a){
Poly ans=Integral(Derivative(a)*Inverse(a));
ans.resize(a.size());
return ans;
} Poly Integral(const Poly& a){
int len=a.size();
Poly ans(len+);
for(int i=;i<len;i++)
ans[i]=1ll*a[i-]*Pow(i,MOD-,MOD)%MOD;
return ans;
} Poly Derivative(const Poly& a){
int len=a.size();
Poly ans(len-);
for(int i=;i<len;i++)
ans[i-]=1ll*a[i]*i%MOD;
return ans;
} Poly operator/(Poly a,Poly b){
int n=a.size()-,m=b.size()-;
Poly ans();
if(n>=m){
std::reverse(a.begin(),a.end());
std::reverse(b.begin(),b.end());
b.resize(n-m+);
ans=Inverse(b)*a;
ans.resize(n-m+);
std::reverse(ans.begin(),ans.end());
}
return ans;
} Poly operator%(Poly a,Poly b){
int n=a.size()-,m=b.size()-;
Poly ans;
if(n<m)
ans=a;
else
ans=a-(a/b)*b;
ans.resize(m);
return ans;
} Poly operator*(Poly a,Poly b){
int len=a.size()+b.size()-;
int bln=NTTPre(len);
NTT(a,bln,DFT);
NTT(b,bln,DFT);
for(int i=;i<bln;i++)
a[i]=1ll*a[i]*b[i]%MOD;
NTT(a,bln,IDFT);
a.resize(len);
return a;
} Poly operator+(const Poly& a,const Poly& b){
Poly ans(std::max(a.size(),b.size()));
std::copy(a.begin(),a.end(),ans.begin());
for(size_t i=;i<b.size();i++)
ans[i]=(ans[i]+b[i])%MOD;
return ans;
} Poly operator-(const Poly& a,const Poly& b){
Poly ans(std::max(a.size(),b.size()));
std::copy(a.begin(),a.end(),ans.begin());
for(size_t i=;i<b.size();i++)
ans[i]=(ans[i]+MOD-b[i])%MOD;
return ans;
} Poly Inverse(Poly a){
int len=a.size();
if(len==)
return Poly(,Pow(a[],MOD-,MOD));
else{
Poly b(a);
b.resize((len+)>>);
b=Inverse(b);
int bln=NTTPre(b.size()*+a.size());
NTT(a,bln,DFT);
NTT(b,bln,DFT);
for(int i=;i<bln;i++)
b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
NTT(b,bln,IDFT);
b.resize(len);
return b;
}
} void NTT(Poly& a,int len,int opt){
a.resize(len);
for(int i=;i<len;i++)
if(rev[i]>i)
std::swap(a[i],a[rev[i]]);
for(int i=;i<len;i<<=){
int step=i<<;
int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
for(int j=;j<len;j+=step){
int w=;
for(int k=;k<i;k++,w=1ll*w*wn%MOD){
int x=a[j+k];
int y=1ll*a[j+k+i]*w%MOD;
a[j+k]=(x+y)%MOD;
a[j+k+i]=(x-y+MOD)%MOD;
}
}
}
if(opt==IDFT){
int inv=Pow(len,MOD-,MOD);
for(int i=;i<len;i++)
a[i]=1ll*a[i]*inv%MOD;
}
} int NTTPre(int n){
int bln=,bct=;
while(bln<n){
bln<<=;
++bct;
}
for(int i=;i<bln;i++)
rev[i]=(rev[i>>]>>)|((i&)<<(bct-));
return bln;
} inline int Pow(int a,int n,int p){
int ans=;
while(n>){
if(n&)
ans=1ll*a*ans%p;
a=1ll*a*a%p;
n>>=;
}
return ans;
}
BZOJ 4555