最主要的两个式子:
套路1:
$$\begin{array}{l}x^k=\sum_{i=0}^k\begin{pmatrix}x\\i\end{pmatrix}\begin{Bmatrix}k\\i\end{Bmatrix}i!\\\end{array}$$
左边的式子可以看成将k个球放到x个有标号的盒子中,总共的方案数。
右边的式子可以看成先选出i个非空的盒子,将k个球放到这i个盒子中且不空的方案数。但由于第二类斯特林数是无序的,还要乘上i阶乘。
套路2:
$$\begin{array}{l}\begin{Bmatrix}n\\m\end{Bmatrix}=\frac1{m!}\;\sum_{i=0}^m\;(-1)^i\;\;\begin{pmatrix}m\\i\end{pmatrix}(m-i)^n\\\end{array}$$
这是第二类斯特林数的容斥形式。
右边相当于枚举有多少个空的盒子,然后剩下的盒子随便放。由于这是有标号的结果,最后要除以m阶乘。
有了套路2,就能NTT在O(nlogn)时间内算出$$\begin{array}{l}\\\begin{Bmatrix}n\\0\end{Bmatrix}\sim\begin{Bmatrix}n\\m\end{Bmatrix}\end{array}$$>。
1.求第二类斯特林数S(n,0)~S(n,m),n,m<=5E5,模998244353。
朴素NTT即可。
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long int ll; 4 const ll mod=998244353; 5 const ll G=3; 6 const ll Gi=332748118; 7 const ll maxn=1E6+5; 8 ll a[maxn],b[maxn],fac[maxn],n,m,len,limit,sum[maxn],ans; 9 int r[maxn]; 10 int re(int x) 11 { 12 int ans=0; 13 for(int i=0;i<len;++i)ans=ans*2+((x&(1<<i))>0); 14 return ans; 15 } 16 ll qpow(ll x,ll y) 17 { 18 ll ans=1,base=x; 19 while(y) 20 { 21 if(y&1)ans=ans*base%mod; 22 base=base*base%mod; 23 y>>=1; 24 } 25 return ans; 26 } 27 void init() 28 { 29 fac[0]=1; 30 for(int i=1;i<=n;++i)fac[i]=fac[i-1]*i%mod; 31 for(int i=0;i<=n;++i)b[i]=qpow(fac[i],mod-2)*qpow(i,n)%mod; 32 for(int i=0;i<=n;++i)a[i]=(qpow(-1,i)*qpow(fac[i],mod-2)%mod+mod)%mod; 33 } 34 void NTT(ll*A,int g) 35 { 36 for(int i=0;i<limit;++i) 37 if(r[i]<i)swap(A[i],A[r[i]]); 38 for(int i=2;i<=limit;i*=2) 39 { 40 ll w; 41 if(g==1)w=qpow(G,(mod-1)/i); 42 else w=qpow(Gi,(mod-1)/i); 43 for(int j=0;j<limit/i;++j) 44 { 45 ll d=1; 46 for(int k=0;k<i/2;++k) 47 { 48 ll a=A[i*j+k],b=A[i*j+i/2+k]*d%mod; 49 A[i*j+k]=(a+b)%mod; 50 A[i*j+i/2+k]=(a-b+mod)%mod; 51 d=d*w%mod; 52 } 53 } 54 } 55 } 56 int main() 57 { 58 ios::sync_with_stdio(false); 59 cin>>n; 60 init(); 61 len=0,limit=1; 62 while(limit<=n*2) 63 { 64 ++len; 65 limit<<=1; 66 } 67 for(int i=n+1;i<=limit;++i)a[i]=0; 68 for(int i=n+1;i<=limit;++i)b[i]=0; 69 for(int i=0;i<limit;++i)r[i]=re(i); 70 NTT(a,1); 71 NTT(b,1); 72 for(int i=0;i<limit;++i)a[i]=a[i]*b[i]%mod; 73 NTT(a,-1); 74 ll g=qpow(limit,mod-2); 75 for(int i=0;i<limit;++i)a[i]=a[i]*g%mod; 76 for(int i=0;i<=n;++i)cout<<a[i]<<" ";cout<<endl; 77 return 0; 78 }
2.n<=1E5,求$$\sum_{i=0}^n\;\sum_{j=0}^i\;\begin{Bmatrix}i\\j\end{Bmatrix}\times2^j\times(j!)\;$$
改变求和顺序,将斯特林数展开即可。
$$\begin{array}{l}\sum_{i=0}^n\;\sum_{j=0}^i\;\begin{Bmatrix}i\\j\end{Bmatrix}\times2^j\times(j!)\;\\=\sum_{j=0}^n\;2^j\times(j!)\sum_{i=j}^n\;\begin{Bmatrix}i\\j\end{Bmatrix}\\=\sum_{j=0}^n\;2^j\times(j!)\sum_{i=0}^n\;\begin{Bmatrix}i\\j\end{Bmatrix}\end{array}$$
令$$\begin{array}{l}\\f(m)=\sum_{i=0}^n\;\begin{Bmatrix}i\\m\end{Bmatrix}\\f(m)=\sum_{i=0}^n\;\frac{\displaystyle1}{\displaystyle m!}\;\sum_{j=0}^n\;(-1)^j\begin{pmatrix}m\\j\end{pmatrix}(m-j)^i\\=\sum_{i=0}^n\sum_{j=0}^n\;\frac{\displaystyle1}{\displaystyle m!}\frac{\displaystyle m!}{\displaystyle j!(m-j)!}(m-j)^i(-1)^j\\=\sum_{i=0}^n\sum_{j=0}^n\;\frac{(m-j)^i}{(m-j)!}\times\frac{(-1)^j}{j!}\\=\sum_{j=0}^n\;\frac{\displaystyle(-1)^j}{\displaystyle j!}\sum_{i=0}^n\;\frac{\displaystyle(m-j)^i}{\displaystyle(m-j)!}\\\end{array}$$
右边只是一个等比数列求和。求出所有的值后,显然可以卷积。
最后答案为$$\sum_{i=0}^n\;2^i\times(i!)\times f(i)$$
注意等比数列求和时若x=1,需特判。
1 // luogu-judger-enable-o2 2 #include<bits/stdc++.h> 3 using namespace std; 4 typedef long long int ll; 5 const ll mod=998244353; 6 const ll G=3; 7 const ll Gi=332748118; 8 const ll maxn=1E6+5; 9 ll a[maxn],b[maxn],fac[maxn],n,m,len,limit,sum[maxn],ans; 10 int r[maxn]; 11 int re(int x) 12 { 13 int ans=0; 14 for(int i=0;i<len;++i)ans=ans*2+((x&(1<<i))>0); 15 return ans; 16 } 17 ll qpow(ll x,ll y) 18 { 19 ll ans=1,base=x; 20 while(y) 21 { 22 if(y&1)ans=ans*base%mod; 23 base=base*base%mod; 24 y>>=1; 25 } 26 return ans; 27 } 28 ll get(ll x){return (qpow(x,n+1)-1)*qpow(x-1,mod-2)%mod;} 29 void init() 30 { 31 fac[0]=1; 32 for(int i=1;i<=n;++i)fac[i]=fac[i-1]*i%mod; 33 for(int i=0;i<=n;++i)a[i]=(qpow(-1,i)*qpow(fac[i],mod-2)%mod+mod)%mod; 34 for(int i=0;i<=n;++i)b[i]=get(i)*qpow(fac[i],mod-2)%mod; 35 } 36 void NTT(ll*A,int g) 37 { 38 for(int i=0;i<limit;++i) 39 if(r[i]<i)swap(A[i],A[r[i]]); 40 for(int i=2;i<=limit;i*=2) 41 { 42 ll w; 43 if(g==1)w=qpow(G,(mod-1)/i); 44 else w=qpow(Gi,(mod-1)/i); 45 for(int j=0;j<limit/i;++j) 46 { 47 ll d=1; 48 for(int k=0;k<i/2;++k) 49 { 50 ll a=A[i*j+k],b=A[i*j+i/2+k]*d%mod; 51 A[i*j+k]=(a+b)%mod; 52 A[i*j+i/2+k]=(a-b+mod)%mod; 53 d=d*w%mod; 54 } 55 } 56 } 57 } 58 int main() 59 { 60 ios::sync_with_stdio(false); 61 cin>>n; 62 init(); 63 len=0,limit=1; 64 while(limit<=n*2) 65 { 66 ++len; 67 limit<<=1; 68 } 69 for(int i=0;i<limit;++i)r[i]=re(i); 70 b[1]=n+1; 71 NTT(a,1); 72 NTT(b,1); 73 for(int i=0;i<limit;++i)a[i]=a[i]*b[i]%mod; 74 NTT(a,-1); 75 ll g=qpow(limit,mod-2); 76 for(int i=0;i<limit;++i)a[i]=a[i]*g%mod; 77 for(int i=0;i<=n;++i)ans=(ans+qpow(2,i)*fac[i]%mod*a[i]%mod)%mod; 78 cout<<ans<<endl; 79 return 0; 80 }
https://www.luogu.org/problemnew/show/P4091
3.n个点,求((所有有编号无向图(每个点度数的k次方的和))的和)。
考虑每个点的贡献,答案为:$$n\;\sum_{i=0}^{n-1}\;i^k\times\begin{pmatrix}n-1\\i\end{pmatrix}\times2^{\begin{pmatrix}n-1\\2\end{pmatrix}}$$
即枚举每个点的度数,n-1个点中要有i个点连向他,剩下的点随便连边。
方便起见,只需要求
$$\begin{array}{l}\sum_{i=0}^n\;i^k\times\begin{pmatrix}n\\i\end{pmatrix}\;\\=\sum_{i=0}^n\;\begin{pmatrix}n\\i\end{pmatrix}\;\sum_{j=0}^k\;\begin{Bmatrix}k\\j\end{Bmatrix}\begin{pmatrix}i\\j\end{pmatrix}j!\\=\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}j!\;\;\sum_{i=0}^n\;\begin{pmatrix}n\\i\end{pmatrix}\begin{pmatrix}i\\j\end{pmatrix}\\=\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}j!\;\;\sum_{i=0}^n\;\begin{pmatrix}n\\j\end{pmatrix}\begin{pmatrix}n-j\\i-j\end{pmatrix}\\=\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}j!\begin{pmatrix}n\\j\end{pmatrix}\;\;\sum_{i=0}^n\;\begin{pmatrix}n-j\\i-j\end{pmatrix}\;\\=\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}j!\begin{pmatrix}n\\j\end{pmatrix}\;\;2^{n-j}\end{array}$$
最后几步看不懂翻翻具体数学。P143
1 #pragma GCC optimize(2) 2 #include<bits/stdc++.h> 3 using namespace std; 4 typedef long long int ll; 5 const ll mod=998244353; 6 const ll G=3; 7 const ll Gi=332748118; 8 const ll maxn=2E6+5; 9 ll min(ll x,ll y){return x<y?x:y;} 10 ll max(ll x,ll y){return x>y?x:y;} 11 int r[maxn],len,limit; 12 ll fac[maxn],a[maxn],b[maxn],k,n,ans,inv[maxn]; 13 int re(int x) 14 { 15 int ans=0; 16 for(int i=0;i<len;++i)ans=ans*2+((x&(1<<i))>0); 17 return ans; 18 } 19 ll qpow(ll x,ll y) 20 { 21 ll ans=1,base=x; 22 while(y) 23 { 24 if(y&1)ans=ans*base%mod; 25 base=base*base%mod; 26 y>>=1; 27 } 28 return ans; 29 } 30 void init() 31 { 32 fac[0]=1; 33 for(int i=1;i<=k;++i)fac[i]=fac[i-1]*i%mod; 34 for(int i=0;i<=k;++i)inv[i]=qpow(fac[i],mod-2); 35 for(int i=0;i<=k;++i)a[i]=(qpow(-1,i)*inv[i]%mod+mod)%mod; 36 for(int i=0;i<=k;++i)b[i]=qpow(i,k)*inv[i]%mod; 37 limit=1; 38 while(limit<=2*k) 39 { 40 ++len; 41 limit<<=1; 42 } 43 for(int i=0;i<limit;++i)r[i]=re(i); 44 } 45 void NTT(ll*A,int g) 46 { 47 for(int i=0;i<limit;++i) 48 if(i<r[i])swap(A[i],A[r[i]]); 49 for(int i=2;i<=limit;i*=2) 50 { 51 ll w; 52 if(g==1)w=qpow(G,(mod-1)/i); 53 else w=qpow(Gi,(mod-1)/i); 54 for(int j=0;j<limit/i;++j) 55 { 56 ll d=1; 57 for(int k=0;k<i/2;++k) 58 { 59 ll a=A[i*j+k],b=d*A[i*j+i/2+k]%mod; 60 A[i*j+k]=(a+b)%mod; 61 A[i*j+i/2+k]=(a-b+mod)%mod; 62 d=d*w%mod; 63 } 64 } 65 } 66 } 67 int main() 68 { 69 ios::sync_with_stdio(false); 70 cin>>n>>k; 71 --n; 72 init(); 73 NTT(a,1); 74 NTT(b,1); 75 for(int i=0;i<=limit;++i)a[i]=a[i]*b[i]%mod; 76 NTT(a,-1); 77 ll g=qpow(limit,mod-2); 78 for(int i=0;i<=limit;++i)a[i]=a[i]*g%mod; 79 ll w=qpow(2,n),GG=qpow(2,mod-2); 80 for(int i=0;i<=min(k,n);++i) 81 { 82 ans=(ans+a[i]*w)%mod; 83 w=w*(n-i)%mod*GG%mod; 84 } 85 ans=ans*(n+1)%mod*qpow(2,n*(n-1)/2)%mod; 86 cout<<ans<<endl; 87 return 0; 88 }
https://www.lydsy.com/JudgeOnline/problem.php?id=5093
4.求树上所有点对距离的k次方和,点对中两个点不一样。n<=5E4,k<=150。
还是展开。
1 // luogu-judger-enable-o2 2 #include<bits/stdc++.h> 3 using namespace std; 4 typedef long long int ll; 5 const ll mod=10007; 6 ll n,f1[50005][255],f2[50005][255],head[50005*2],size,x,y,s[555][555],fac[555],k,tmp[555]; 7 struct edge{int to,next;}E[50005*2]; 8 void add(int u,int v) 9 { 10 E[++size].to=v; 11 E[size].next=head[u]; 12 head[u]=size; 13 } 14 void init() 15 { 16 s[0][0]=1; 17 for(int i=1;i<=200;++i) 18 for(int j=1;j<=200;++j) 19 s[i][j]=(j*s[i-1][j]%mod+s[i-1][j-1])%mod; 20 fac[0]=1; 21 for(int i=1;i<=200;++i)fac[i]=fac[i-1]*i%mod; 22 } 23 void dfs1(int u,int F) 24 { 25 f1[u][0]=1; 26 for(int e=head[u];e;e=E[e].next) 27 { 28 int v=E[e].to; 29 if(v==F)continue; 30 dfs1(v,u); 31 for(int i=k;i>=1;--i)f1[u][i]=(f1[u][i]+f1[v][i]+f1[v][i-1])%mod; 32 f1[u][0]=(f1[u][0]+f1[v][0])%mod; 33 } 34 } 35 void dfs2(int u,int F) 36 { 37 for(int i=0;i<=k;++i)f2[u][i]=f1[u][i]; 38 if(F!=0) 39 { 40 for(int i=1;i<=k;++i) 41 tmp[i]=(f2[F][i]-f1[u][i]-f1[u][i-1]+mod*2)%mod; 42 tmp[0]=(f2[F][0]-f1[u][0]+mod)%mod; 43 for(int i=1;i<=k;++i) 44 f2[u][i]=(f2[u][i]+tmp[i]+tmp[i-1])%mod; 45 f2[u][0]=(f2[u][0]+tmp[0])%mod; 46 } 47 for(int e=head[u];e;e=E[e].next) 48 { 49 int v=E[e].to; 50 if(v==F)continue; 51 dfs2(v,u); 52 } 53 } 54 int main() 55 { 56 // freopen("a.in","r",stdin); 57 ios::sync_with_stdio(false); 58 cin>>n>>k; 59 for(int i=1;i<=n-1;++i) 60 { 61 cin>>x>>y; 62 add(x,y); 63 add(y,x); 64 } 65 init(); 66 dfs1(1,0); 67 dfs2(1,0); 68 for(int u=1;u<=n;++u) 69 { 70 int ans=0; 71 for(int i=0;i<=k;++i)ans=(ans+s[k][i]*fac[i]%mod*f2[u][i]%mod)%mod; 72 cout<<ans<<endl; 73 } 74 return 0; 75 }
https://www.luogu.org/problemnew/show/P4827