「学习笔记」斯特林数

「学习笔记」斯特林数

最近做题感觉两类斯特林数挺有用的，特地总结一下。

第一类斯特林数：

\(S(n,m)\) 表示将一个有 \(n\) 个数的序列划分成 \(m\) 个圆排列的方案数。

\(S(n,m)=S(n-1,m-1)+(n-1)\times S(n-1,m)\)

若有 \(m-1\) 个圆排列，那么直接单独组成一个圆排列就行了。

若有 \(m\) 个圆排列，那么有 \(n-1\) 个位置可放。

若 \(n\leq 10^5\) 怎么办呢？

考虑它的生成函数

\[\prod_{i=0}^{n-1}(x+i) \]

可以用分治 \(FFT\)，也可以倍增。

所以其实“在 \(n\) 个选择 \(m\) 个数相乘的和”这个问题可以用类似第一类斯特林数的生成函数解决。

放两道例题：

「FJOI2016」建筑师

我们从最大值剪开，变成第一类斯特林数的问题。

\[ans=C(A+B-2,A-1)\times S(n-1,A+B-2) \]

\(Code\ Below:\)

#include <bits/stdc++.h>
using namespace std;
const int maxn=50000+10;
const int mod=1e9+7;
int n,A,B,fac[maxn],inv[maxn],S[maxn][210];

inline int C(int n,int m){
	if(n<m) return 0;
	return 1ll*fac[n]*inv[m]%mod*inv[n-m]%mod;
}

int main()
{
	fac[0]=fac[1]=inv[0]=inv[1]=1;
	for(int i=2;i<=200;i++) fac[i]=1ll*fac[i-1]*i%mod;
	for(int i=2;i<=200;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
	for(int i=2;i<=200;i++) inv[i]=1ll*inv[i]*inv[i-1]%mod;
	S[0][0]=1;
	for(int i=1;i<=50000;i++)
		for(int j=1;j<=200;j++) S[i][j]=(S[i-1][j-1]+1ll*(i-1)*S[i-1][j]%mod)%mod;
	int T;
	scanf("%d",&T);
	while(T--){
		scanf("%d%d%d",&n,&A,&B);
		printf("%d\n",1ll*C(A+B-2,A-1)*S[n-1][A+B-2]%mod);
	}
	return 0;
}

CF960G Bandit Blues

可以用分治 \(FFT\)，时间复杂度 \(O(n\log^2 n)\)

\(Code\ Below:\)

#include <bits/stdc++.h>
using namespace std;
const int maxn=400000+10;
const int mod=998244353;
int n,A,B,fac[maxn],inv[maxn],r[maxn],G[40][2];
vector<int> P[maxn];

inline int C(int n,int m){
    if(n<m) return 0;
    return 1ll*fac[n]*inv[m]%mod*inv[n-m]%mod;
}

inline int fpow(int a,int b){
	int ret=1;
	for(;b;b>>=1,a=1ll*a*a%mod)
		if(b&1) ret=1ll*ret*a%mod;
	return ret;
}

inline void Gpre(){
	for(int len=1,l=1;len<=mod;len<<=1,l++){
		G[l][0]=fpow(3,(mod-1)/(len<<1));
		G[l][1]=fpow(G[l][0],mod-2);
	}
}

inline void NTT(int *f,int n,int op){
	for(int i=0;i<n;i++)
		if(i<r[i]) swap(f[i],f[r[i]]);
	int buf,tmp,x,y;
	for(int len=1,l=1;len<n;len<<=1,l++){
		tmp=(op==1)?G[l][0]:G[l][1];
		for(int i=0;i<n;i+=len<<1){
			buf=1;
			for(int j=0;j<len;j++){
				x=f[i+j];y=1ll*buf*f[i+j+len]%mod;
				f[i+j]=(x+y)%mod;f[i+j+len]=(x-y+mod)%mod;
				buf=1ll*buf*tmp%mod;
			}
		}
	}
	if(op==1) return ;
	int inv=fpow(n,mod-2);
	for(int i=0;i<n;i++) f[i]=1ll*f[i]*inv%mod;
}

inline void Mul(int *A,int *B,int *C,int n,int m){
	int lim;
	for(lim=1;lim<(n+m);lim<<=1);
	for(int i=n;i<lim;i++) A[i]=0;
	for(int i=m;i<lim;i++) B[i]=0;
	for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)?(lim>>1):0);
	NTT(A,lim,1);NTT(B,lim,1);
	for(int i=0;i<lim;i++) C[i]=1ll*A[i]*B[i]%mod;
	NTT(C,lim,-1);
}

inline void solve(int l,int r,int x){
	if(l==r){
		P[x].push_back(l);
		P[x].push_back(1);
		return ;
	}
	int mid=(l+r)>>1;
	solve(l,mid,x<<1);solve(mid+1,r,x<<1|1);
	int n=P[x<<1].size(),m=P[x<<1|1].size();
	static int A[maxn],B[maxn],C[maxn];
	for(int i=0;i<n;i++) A[i]=P[x<<1][i];
	for(int i=0;i<m;i++) B[i]=P[x<<1|1][i];
	Mul(A,B,C,n,m);
	for(int i=0;i<n+m-1;i++) P[x].push_back(C[i]);
	P[x<<1].clear();P[x<<1|1].clear();
}

int main()
{
	Gpre();
    fac[0]=fac[1]=inv[0]=inv[1]=1;
    for(int i=2;i<=maxn-10;i++) fac[i]=1ll*fac[i-1]*i%mod;
    for(int i=2;i<=maxn-10;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    for(int i=2;i<=maxn-10;i++) inv[i]=1ll*inv[i]*inv[i-1]%mod;
    scanf("%d%d%d",&n,&A,&B);
    if(A==0||B==0||A+B-1>n){
    	printf("0\n");
    	return 0;
	}
    if(n==1){
    	printf("1\n");
    	return 0;
	}
    solve(0,n-2,1);
    printf("%d\n",1ll*C(A+B-2,A-1)*P[1][A+B-2]%mod);
    return 0;
}

第二类斯特林数

\(S(n,m)\) 表示将 \(n\) 个不同的小球放到 \(m\) 个相同的箱子内的方案数。

\[S(n,m)=S(n-1,m-1)+m\times S(n-1,m) \]

若有 \(m-1\) 个箱子，那么直接单独占有一个箱子就行了。

若有 \(m\) 个箱子，那么有 \(m\) 个位置可放。

若 \(n\leq 10^5\) 怎么办？

考虑一个神奇的式子

\[S(i,j)=\frac {1}{j!}\sum_{k=0}^{j}(-1)^k{j\choose k}(j-k)^{i} \]

容斥原理。虽然我不会证，不过会用就行了。上述式子用 \(NTT\) 优化即可。

还有一个神奇的式子

\[n^k=\sum_{i=0}^{k}S(k,i)\times i!\times {n\choose i} \]

左边就是把 \(k\) 个球任意放到 \(n\) 个箱子。

右边就是枚举放了几个非空箱子，哪几个非空箱子，再乘上第二类斯特林数。不过这里的箱子是不同的，所以要乘上 \(i!\)

「国家集训队」Crash 的文明世界

树形 \(dp\) 好题！

这题还有用到

\[{i\choose j}={i-1\choose j-1}+{i-1\choose j} \]

好巧啊！

\(Code\ Below:\)

// luogu-judger-enable-o2
#include <bits/stdc++.h>
using namespace std;
const int maxn=50000+10;
const int maxm=150+10;
const int mod=10007;
int n,k,dp[maxn][maxm],now[maxm],fac[maxm],S[maxm][maxm];
int head[maxn],to[maxn<<1],nxt[maxn<<1],tot;

inline int read(){
	register int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return (f==1)?x:-x;
}

inline void addedge(int x,int y){
	to[++tot]=y;
	nxt[tot]=head[x];
	head[x]=tot;
}

void dfs1(int x,int f){
	dp[x][0]=1;
	for(int i=head[x],y;i;i=nxt[i]){
		y=to[i];
		if(y==f) continue;
		dfs1(y,x);dp[x][0]=(dp[x][0]+dp[y][0])%mod;
		for(int j=1;j<=k;j++) dp[x][j]=(dp[x][j]+dp[y][j]+dp[y][j-1])%mod;
	}
}

void dfs2(int x,int f){
	for(int i=head[x],y;i;i=nxt[i]){
		y=to[i];
		if(y==f) continue;
		now[0]=(dp[x][0]-dp[y][0])%mod;
		for(int j=1;j<=k;j++) now[j]=(dp[x][j]-dp[y][j]-dp[y][j-1])%mod;
		dp[y][0]=(dp[y][0]+now[0])%mod;
		for(int j=1;j<=k;j++) dp[y][j]=(dp[y][j]+now[j]+now[j-1])%mod;
		dfs2(y,x);
	}
}

int main()
{
	n=read(),k=read();
	int x,y;
	for(int i=1;i<n;i++){
		x=read(),y=read();
		addedge(x,y);addedge(y,x);
	}
	fac[0]=S[0][0]=1;
	for(int i=1;i<=k;i++) fac[i]=fac[i-1]*i%mod;
	for(int i=1;i<=k;i++)
		for(int j=1;j<=i;j++) S[i][j]=(S[i-1][j-1]+S[i-1][j]*j)%mod;
	dfs1(1,0);dfs2(1,0);
	int ans;
	for(int i=1;i<=n;i++){
		ans=0;
		for(int j=0;j<=k;j++) ans=(ans+S[k][j]*fac[j]%mod*dp[i][j])%mod;
		ans=(ans+mod)%mod;
		printf("%d\n",ans);
	}
	return 0;
}

CF932E Team Work

咕咕咕。

\(Code\ Below:\)

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=5000+10;
const int mod=1e9+7;
int n,k,S[maxn][maxn];

int fast_pow(int a,int b){
	int ret=1;
	for(;b;b>>=1,a=1ll*a*a%mod)
		if(b&1) ret=1ll*ret*a%mod;
	return ret;
}

int main()
{
	scanf("%d%d",&n,&k);
	S[0][0]=1;
	for(int i=1;i<=k;i++)
		for(int j=1;j<=i;j++) S[i][j]=(S[i-1][j-1]+1ll*S[i-1][j]*j%mod)%mod;
	int now=fast_pow(2,n),inv2=fast_pow(2,mod-2),tmp=1,ans=0;
	for(int i=0;i<=k;i++){
		ans=(ans+1ll*S[k][i]*now%mod*tmp%mod)%mod;
		now=1ll*now*inv2%mod;tmp=1ll*tmp*(n-i)%mod;
	}
	printf("%d\n",ans);
	return 0;
}

CF1097G Vladislav and a Great Legend

树形 \(dp\) 好题！

\[f(X)^k=\sum_{i=0}^{k}S(k,i)\times i!\times {f(X)\choose i} \]

考虑 \(\Large{f(X)\choose i}\) 的组合意义，可以树形 \(dp\)

\(dp[i][j]\) 表示 \(i\) 结点选了 \(j\) 条边的答案。我们在深度最小的节点上计算贡献，每一次去合并两个子树的信息。

\[f[i+j]=\sum_dp[x][i]\times dp[y][j] \]

我本来以为可以用 \(NTT\) 优化的，结果模数没原根，\(k\leq 200\)。。。

不过有一些细节需要注意的，比如开始 \(dp[x][0]=2\)（选/不选），最后 \(dp[x][1]\) 要 \(-1\)（子树空了就不行）

\(Code\ Below:\)

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=100000+10;
const int mod=1e9+7;
int n,k;ll dp[maxn][210],f[210],ans[210],S[210][210],fac[maxn];
int siz[maxn],head[maxn],to[maxn<<1],nxt[maxn<<1],tot;

inline int read(){
	register int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
	return (f==1)?x:-x;
}

inline void addedge(int x,int y){
	to[++tot]=y;
	nxt[tot]=head[x];
	head[x]=tot;
}

void dfs(int x,int fa){
	siz[x]=1;dp[x][0]=2;
	for(int i=head[x],y;i;i=nxt[i]){
		y=to[i];
		if(y==fa) continue;
		dfs(y,x);
		for(int i=0;i<=k;i++) f[i]=0;
		for(int i=0;i<=min(k,siz[x]);i++)
			for(int j=0;j<=min(k-i,siz[y]);j++) f[i+j]=(f[i+j]+dp[x][i]*dp[y][j])%mod;
		siz[x]+=siz[y];
		for(int i=0;i<=k;i++) dp[x][i]=f[i];
		for(int i=0;i<=k;i++) ans[i]=(ans[i]-dp[y][i]+mod)%mod;
	}
	for(int i=0;i<=k;i++) ans[i]=(ans[i]+dp[x][i])%mod;
	for(int i=k;i>=1;i--) dp[x][i]=(dp[x][i]+dp[x][i-1])%mod;
	dp[x][1]=(dp[x][1]-1+mod)%mod;
}

int main()
{
	n=read(),k=read();
	int x,y;
	for(int i=1;i<n;i++){
		x=read(),y=read();
		addedge(x,y);addedge(y,x);
	}
	fac[0]=1;S[0][0]=1;
	for(ll i=1;i<=k;i++) fac[i]=fac[i-1]*i%mod;
	for(ll i=1;i<=k;i++)
		for(ll j=1;j<=i;j++) S[i][j]=(S[i-1][j-1]+S[i-1][j]*j)%mod;
	dfs(1,0);
	ll Ans=0;
	for(int i=0;i<=k;i++) Ans=(Ans+S[k][i]*fac[i]%mod*ans[i])%mod;
	printf("%I64d\n",Ans);
	return 0;
}