[BZOJ 4555][Tjoi2016&Heoi2016]求和
题意
给定 $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要除以阶乘...出答案的时候还得乘回来...
1 #include <bits/stdc++.h> 2 3 const int G=3; 4 const int DFT=1; 5 const int IDFT=-1; 6 const int MAXN=550000; 7 const int MOD=998244353; 8 const int INV2=(MOD+1)>>1; 9 const int PHI=MOD-1; 10 11 typedef std::vector<int> Poly; 12 13 Poly Sqrt(Poly); 14 void Read(Poly&); 15 Poly Inverse(Poly); 16 Poly Ln(const Poly&); 17 Poly Exp(const Poly&); 18 void Print(const Poly&); 19 void NTT(Poly&,int,int); 20 Poly Pow(const Poly&,int); 21 Poly Integral(const Poly&); 22 Poly Derivative(const Poly&); 23 Poly operator*(Poly,Poly); 24 Poly operator/(Poly,Poly); 25 Poly operator%(Poly,Poly); 26 Poly operator+(const Poly&,const Poly&); 27 Poly operator-(const Poly&,const Poly&); 28 29 int rev[MAXN]; 30 31 int NTTPre(int); 32 int Pow(int,int,int); 33 34 int main(){ 35 int n; 36 scanf("%d",&n); 37 Poly h(n+1),one(1,1); 38 h[1]=2; 39 for(int i=2;i<=n;i++) 40 h[i]=1ll*h[i-1]*Pow(i,MOD-2,MOD)%MOD; 41 Poly g(Inverse(one-h)); 42 int ans=0; 43 int fact=1; 44 for(int i=0;i<=n;i++){ 45 (ans+=1ll*g[i]*fact%MOD)%=MOD; 46 fact=1ll*fact*(i+1)%MOD; 47 } 48 printf("%d\n",ans); 49 return 0; 50 } 51 52 void Read(Poly& a){ 53 for(auto& i:a) 54 scanf("%d",&i); 55 } 56 57 void Print(const Poly& a){ 58 for(auto i:a) 59 printf("%d ",i); 60 puts(""); 61 } 62 63 Poly Pow(const Poly& a,int k){ 64 Poly log=Ln(a); 65 for(auto& i:log) 66 i=1ll*i*k%MOD; 67 return Exp(log); 68 } 69 70 Poly Sqrt(Poly a){ 71 int len=a.size(); 72 if(len==1) 73 return Poly(1,int(sqrt(a[0]))); 74 else{ 75 Poly b=a; 76 b.resize((len+1)>>1); 77 b=Sqrt(b); 78 b.resize(len); 79 Poly inv=Inverse(b); 80 int bln=NTTPre(inv.size()+a.size()); 81 NTT(a,bln,DFT); 82 NTT(inv,bln,DFT); 83 for(int i=0;i<bln;i++) 84 a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD; 85 NTT(a,bln,IDFT); 86 for(int i=0;i<len;i++) 87 b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD; 88 return b; 89 } 90 } 91 92 Poly Exp(const Poly& a){ 93 size_t len=1; 94 Poly ans(1,1),one(1,1); 95 while(len<(a.size()<<1)){ 96 len<<=1; 97 Poly b=a; 98 b.resize(len); 99 ans=ans*(one-Ln(ans)+b); 100 ans.resize(len); 101 } 102 ans.resize(a.size()); 103 return ans; 104 } 105 106 Poly Ln(const Poly& a){ 107 Poly ans=Integral(Derivative(a)*Inverse(a)); 108 ans.resize(a.size()); 109 return ans; 110 } 111 112 Poly Integral(const Poly& a){ 113 int len=a.size(); 114 Poly ans(len+1); 115 for(int i=1;i<len;i++) 116 ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD; 117 return ans; 118 } 119 120 Poly Derivative(const Poly& a){ 121 int len=a.size(); 122 Poly ans(len-1); 123 for(int i=1;i<len;i++) 124 ans[i-1]=1ll*a[i]*i%MOD; 125 return ans; 126 } 127 128 Poly operator/(Poly a,Poly b){ 129 int n=a.size()-1,m=b.size()-1; 130 Poly ans(1); 131 if(n>=m){ 132 std::reverse(a.begin(),a.end()); 133 std::reverse(b.begin(),b.end()); 134 b.resize(n-m+1); 135 ans=Inverse(b)*a; 136 ans.resize(n-m+1); 137 std::reverse(ans.begin(),ans.end()); 138 } 139 return ans; 140 } 141 142 Poly operator%(Poly a,Poly b){ 143 int n=a.size()-1,m=b.size()-1; 144 Poly ans; 145 if(n<m) 146 ans=a; 147 else 148 ans=a-(a/b)*b; 149 ans.resize(m); 150 return ans; 151 } 152 153 Poly operator*(Poly a,Poly b){ 154 int len=a.size()+b.size()-1; 155 int bln=NTTPre(len); 156 NTT(a,bln,DFT); 157 NTT(b,bln,DFT); 158 for(int i=0;i<bln;i++) 159 a[i]=1ll*a[i]*b[i]%MOD; 160 NTT(a,bln,IDFT); 161 a.resize(len); 162 return a; 163 } 164 165 Poly operator+(const Poly& a,const Poly& b){ 166 Poly ans(std::max(a.size(),b.size())); 167 std::copy(a.begin(),a.end(),ans.begin()); 168 for(size_t i=0;i<b.size();i++) 169 ans[i]=(ans[i]+b[i])%MOD; 170 return ans; 171 } 172 173 Poly operator-(const Poly& a,const Poly& b){ 174 Poly ans(std::max(a.size(),b.size())); 175 std::copy(a.begin(),a.end(),ans.begin()); 176 for(size_t i=0;i<b.size();i++) 177 ans[i]=(ans[i]+MOD-b[i])%MOD; 178 return ans; 179 } 180 181 Poly Inverse(Poly a){ 182 int len=a.size(); 183 if(len==1) 184 return Poly(1,Pow(a[0],MOD-2,MOD)); 185 else{ 186 Poly b(a); 187 b.resize((len+1)>>1); 188 b=Inverse(b); 189 int bln=NTTPre(b.size()*2+a.size()); 190 NTT(a,bln,DFT); 191 NTT(b,bln,DFT); 192 for(int i=0;i<bln;i++) 193 b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD; 194 NTT(b,bln,IDFT); 195 b.resize(len); 196 return b; 197 } 198 } 199 200 void NTT(Poly& a,int len,int opt){ 201 a.resize(len); 202 for(int i=0;i<len;i++) 203 if(rev[i]>i) 204 std::swap(a[i],a[rev[i]]); 205 for(int i=1;i<len;i<<=1){ 206 int step=i<<1; 207 int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD); 208 for(int j=0;j<len;j+=step){ 209 int w=1; 210 for(int k=0;k<i;k++,w=1ll*w*wn%MOD){ 211 int x=a[j+k]; 212 int y=1ll*a[j+k+i]*w%MOD; 213 a[j+k]=(x+y)%MOD; 214 a[j+k+i]=(x-y+MOD)%MOD; 215 } 216 } 217 } 218 if(opt==IDFT){ 219 int inv=Pow(len,MOD-2,MOD); 220 for(int i=0;i<len;i++) 221 a[i]=1ll*a[i]*inv%MOD; 222 } 223 } 224 225 int NTTPre(int n){ 226 int bln=1,bct=0; 227 while(bln<n){ 228 bln<<=1; 229 ++bct; 230 } 231 for(int i=0;i<bln;i++) 232 rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1)); 233 return bln; 234 } 235 236 inline int Pow(int a,int n,int p){ 237 int ans=1; 238 while(n>0){ 239 if(n&1) 240 ans=1ll*a*ans%p; 241 a=1ll*a*a%p; 242 n>>=1; 243 } 244 return ans; 245 }
本博客已弃用, 新个人主页: https://rvalue.moe, 新博客: https://blog.rvalue.moe