斯特林数与斯特林反演
第一类斯特林数
定义:
$\Large {n \brack m} $ 表示 \(n\) 个元素分成 \(m\) 个环的方案数
显然:
理解:考虑从 \(n−1\) 个元素推过来,如果两个空环肯定是不符合的空一个环则单独成环,如果 \(n−1\) 的时候就没有空环就任意放在一个元素前
性质:
1:
理解: 其实本质就是置换,一个环则为一组轮换,每种排列都会对应着唯一 一种置换
2:
使用归纳法即可:
类似地,我们可以证明:
证明类上,不再赘述。
第二类斯特林数
定义:
\(\Large {n\brace m}\) 表示 \(n\) 个有区别的小球丢进 \(m\) 个无区别的盒子,无空盒子的方案数
显然:
理解: 考虑从 \(n−1\) 个小球推过来,如果两个空盒子肯定是不符合的空一个盒子则只能放到那个空盒子里面了,如果 \(n−1\) 的时候就没有空箱子就随便放。
性质:
当然也可以写成:
理解: \(m^n\) 为 \(n\) 个有区别的小球丢进 \(m\) 个有区别的盒子,允许空盒子。枚举有效盒子的个数,再从 \(m\) 个盒子选 \(i\) 个盒子,然后 \(n\) 个小球丢进 \(i\) 个盒子。
转换到组合表示:
第二类斯特林数显然是和排列组合有关系的,转换过来:
理解: 如果空箱子的情况我们也算进去,答案显然是 \(\dfrac {m^n}{m!}\) ,反过来求第二类斯特林数,又得减掉这种情况:选 \(k\) 个空盒子,然后小球放到其他的盒子里,但最后我们求出来的答案为有区别的盒子,转换过来要 \(×\dfrac 1 {m!}\)
第二类斯特林数与自然数幂的关系
求:
有:
其中 \(\sum_{i=0}^n C_i^j\) 是求杨辉三角一列的和。
通过归纳法可得: \(\sum_{i=0}^n C_i^j=C_{n+1}^{j+1}\) ,进而有:
例题:
[省选联考 2020 A 卷] 组合数问题
计算:
其中 \(n\), \(x\) \((1\le n, x\le 10^9)\) 为给定的整数,\(f(k)\) 为给定的一个 \(m\) 次多项式 \(f(k) = a_0 + a_1k + a_2k^2 + \cdots + a_mk^m\),\(0\le m \le \min(n,1000)\)
有:
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,x,md,m;
inline int pwr(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%md;
x=1ll*x*x%md;y>>=1;
}return res;
}
int s[1005][1005],a[1005];
inline void init(){
s[0][0]=1;
for(int i=1;i<=m;i++){
for(int j=1;j<=i;j++)s[i][j]=(s[i-1][j-1]+1ll*s[i-1][j]*j)%md;
}
}
int main(){
scanf("%d%d%d%d",&n,&x,&md,&m);init();
for(int i=0;i<=m;i++)scanf("%d",&a[i]);
long long ans=0;
for(int i=0;i<=m;i++){
long long res=0;
for(int j=i;j<=m;j++)res=(res+1ll*a[j]*s[j][i])%md;
res=res*pwr(x,i)%md*pwr(x+1,n-i)%md;
for(int j=1;j<=i;j++)res=res*(n-j+1)%md;ans=(ans+res)%md;
}
printf("%lld",ans);
return 0;
}
[清华集训2016] 如何优雅地求和
有一个多项式函数 \(f(x)\),最高次幂为 \(x^m\) ,定义变换 \(Q\) :
现在给定函数 \(f\) 和 \(n,x\) ,求 \(Q(f,n,x)\bmod998244353\) 。
出于某种原因,函数 \(f\) 由点值形式给出,即给定 \(a_0,a_1,⋯,a_m\) 共 \(m+1\) 个数,\(f(x)=a_x\) 。
考虑二项式反演:
由于 \(\Large{x\choose i}\) 是一个 \(i\) 次多项式,\(f(x)\) 可以写成 \(f(x)=\sum_{i=0}^m\limits s_i{x\choose i}\) 的形式。
有:
卷积即可求出 \(s\) ,将 \(f(x)=\sum_{i=0}^m\limits s_i{x\choose i}\) ,代入原式中,有:
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,m,x;
const int md=998244353,G=3,Gi=332748118;
int r[100005];
inline int pwr(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%md;
x=1ll*x*x%md;y>>=1;
}return res;
}
inline void NTT(int *dp,int N,int W){
int lim=0;while((1<<lim)<=N)lim++;
for(int i=0;i<(1<<lim);i++)r[i]=(r[i>>1]>>1)+((i&1)<<(lim-1));
for(int i=0;i<(1<<lim);i++)if(i<r[i])swap(dp[i],dp[r[i]]);
for(int i=0;i<lim;i++){
int w=pwr(W,(md-1)/(1<<(i+1)));
for(int j=0;j<(1<<lim);j+=(1<<(i+1))){
int Pw=1;
for(int t=0;t<(1<<i);t++,Pw=1ll*Pw*w%md){
int x=dp[j+t],y=1ll*Pw*dp[j+(1<<i)+t]%md;
dp[j+t]=(x+y)%md;dp[j+(1<<i)+t]=(x-y+md)%md;
}
}
}
}
int fac[100005],inv[100005],a[100005],s[100005];
inline void init(){
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=2;i<=2*m;i++)fac[i]=1ll*fac[i-1]*i%md;
for(int i=2;i<=2*m;i++)inv[i]=1ll*(md-md/i)*inv[md%i]%md;
for(int i=2;i<=2*m;i++)inv[i]=1ll*inv[i]*inv[i-1]%md;
}
int main(){
scanf("%d%d%d",&n,&m,&x);init();
for(int i=0;i<=m;i++)scanf("%d",&a[i]);
for(int i=0;i<=m;i++)a[i]=1ll*a[i]*inv[i]%md;
for(int i=0;i<=m;i++)s[i]=(i&1?md-inv[i]:inv[i]);
NTT(a,2*m+2,G);NTT(s,2*m+2,G);
for(int i=0;i<=4*m+4;i++)s[i]=1ll*a[i]*s[i]%md;
NTT(s,2*m+2,Gi);
int lim=0;while((1<<lim)<=2*m+2)lim++;lim=pwr(1<<lim,md-2);
for(int i=0;i<=m;i++)s[i]=1ll*s[i]*fac[i]%md*lim%md;
long long ans=0,Pw=1;
for(int i=0;i<=m;i++){
ans=(ans+1ll*s[i]*Pw%md*inv[i])%md;
Pw=Pw*(n-i)%md*x%md;
}
printf("%lld",ans);
return 0;
}
CF961G Partitions
给出 \(n\) 个物品,每个物品有一个权值 \(w_i\) ,定义一个集合 \(S\) 的权值 \(W(S)=|S|\sum\limits_{x\in S}w_x\) 。
定义一个划分的权值为 \(W'(R)=\sum\limits_{S\in R}W(S)\) ,求将 \(n\) 个物品划分成 \(k\) 个集合的所有方案的权值和。
考虑每个物品的贡献,对于每种方案,它都会和自己产生 \(w_i\times\Large{n\brace k}\) 的贡献,另外,对于剩下 \(n-1\) 个物品的每种划分方案,它都可以和每个物品产生贡献,共计 \(w_i\times(n-1)\times \Large{n-1\brace k}\)
我们要求的就是:
有容斥公式:
\(O(n+k)\) 求解即可。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,k;
int w[200005],fac[200005],inv[200005];
const int md=1e9+7;
inline void init(){
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=2;i<=n;i++)fac[i]=1ll*fac[i-1]*i%md;
for(int i=2;i<=n;i++)inv[i]=1ll*(md-md/i)*inv[md%i]%md;
for(int i=2;i<=n;i++)inv[i]=1ll*inv[i]*inv[i-1]%md;
}
inline int C(int x,int y){
if(x<y)return 0;
return 1ll*fac[x]*inv[y]%md*inv[x-y]%md;
}
inline int pwr(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%md;
x=1ll*x*x%md;y>>=1;
}return res;
}
inline int S(int x,int y){
int res=0;
for(int i=0;i<=y;i++)res=(res+(i&1?-1ll:1ll)*C(y,i)*pwr(y-i,x)%md)%md;
return 1ll*res*inv[y]%md;
}
int main(){
scanf("%d%d",&n,&k);init();
for(int i=1;i<=n;i++)scanf("%d",&w[i]);
int sum=0;
for(int i=1;i<=n;i++)sum=(sum+w[i])%md;
printf("%lld",(1ll*sum*(S(n,k)+(n-1ll)*S(n-1,k)%md)%md+md)%md);
return 0;
}
[国家集训队] Crash 的文明世界
给定一棵树,对于每个点 \(x\),求 \(\sum_{i=1}^n\limits dist(i,x)^k\) 。
做法1.
直接维护 \(0\) ~ \(k\) 次方的和单次转移是 \(k^2\) 的吗,考虑维护 \(0\) ~ \(k\) 次下降幂。
转移时考虑到:
因此,考虑转移 \(u\to x\),有:
最终答案为:
做法2.
对每个 \(x\) ,对于 \(j\in [0,k]\) 维护:
考虑到:
有转移:
最终答案为:
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,k;
int ver[100005],ne[100005],head[100005],cnt;
inline void link(int x,int y){
ver[++cnt]=y;
ne[cnt]=head[x];
head[x]=cnt;
}
const int md=10007;
int s[155][155],dp[100005][155],fac[100005];
void dfs1(int x,int fi){
dp[x][0]=1;
for(int i=head[x];i;i=ne[i]){
int u=ver[i];
if(u==fi)continue;
dfs1(u,x);
for(int j=1;j<=k;j++)dp[x][j]=(dp[x][j]+dp[u][j]+dp[u][j-1])%md;
dp[x][0]=(dp[x][0]+dp[u][0])%md;
}
}
int ans[100005];
void dfs2(int x,int fi){
for(int j=0;j<=k;j++)ans[x]=(ans[x]+dp[x][j]*fac[j]%md*s[k][j])%md;
for(int i=head[x];i;i=ne[i]){
int u=ver[i];
if(u==fi)continue;
for(int j=1;j<=k;j++)dp[x][j]=(dp[x][j]-dp[u][j]-dp[u][j-1])%md;
dp[x][0]=(dp[x][0]-dp[u][0])%md;
for(int j=1;j<=k;j++)dp[u][j]=(dp[u][j]+dp[x][j]+dp[x][j-1])%md;
dp[u][0]=(dp[u][0]+dp[x][0])%md;
dfs2(u,x);
for(int j=1;j<=k;j++)dp[u][j]=(dp[u][j]-dp[x][j]-dp[x][j-1])%md;
dp[u][0]=(dp[u][0]-dp[x][0])%md;
for(int j=1;j<=k;j++)dp[x][j]=(dp[x][j]+dp[u][j]+dp[u][j-1])%md;
dp[x][0]=(dp[x][0]+dp[u][0])%md;
}
}
int main(){
scanf("%d%d",&n,&k);
for(int i=1;i<n;i++){
int x,y;scanf("%d%d",&x,&y);
link(x,y);link(y,x);
}
fac[0]=fac[1]=s[0][0]=s[1][1]=1;
for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%md;
for(int i=2;i<=k;i++){
for(int j=1;j<=i;j++)s[i][j]=(1ll*j*s[i-1][j]+s[i-1][j-1])%md;
}
dfs1(1,1);
dfs2(1,1);
for(int i=1;i<=n;i++)printf("%d\n",(ans[i]+md)%md);
return 0;
}
[HEOI2016/TJOI2016]求和
计算:
有:
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n;
int inv[400005],fac[400005];
const int md=998244353,G=3,Gi=332748118;
inline void init(){
inv[0]=inv[1]=fac[0]=fac[1]=1;
for(int i=2;i<=n;i++)fac[i]=1ll*fac[i-1]*i%md;
for(int i=2;i<=n;i++)inv[i]=1ll*(md-md/i)*inv[md%i]%md;
for(int i=2;i<=n;i++)inv[i]=1ll*inv[i]*inv[i-1]%md;
}
int r[400005];
inline int pwr(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%md;
x=1ll*x*x%md;y>>=1;
}return res;
}
inline void NTT(int *dp,int N,int W){
int lim=0;while((1<<lim)<=N)lim++;
for(int i=0;i<(1<<lim);i++)r[i]=(r[i>>1]>>1)+((i&1)<<(lim-1));
for(int i=0;i<(1<<lim);i++)if(i<r[i])swap(dp[i],dp[r[i]]);
for(int i=0;i<lim;i++){
int w=pwr(W,(md-1)/(1<<(i+1)));
for(int j=0;j<(1<<lim);j+=(1<<(i+1))){
int Pw=1;
for(int t=0;t<(1<<i);t++){
int x=dp[j+t],y=1ll*Pw*dp[j+(1<<i)+t]%md;
dp[j+t]=(x+y)%md;dp[j+(1<<i)+t]=(x-y+md)%md;Pw=1ll*Pw*w%md;
}
}
}
}
int a[400005],b[400005];
int main(){
scanf("%d",&n);init();
for(int i=0;i<=n;i++)a[i]=(i&1?md-1ll:1ll)*inv[i]%md;
for(int i=0;i<=n;i++)b[i]=1ll*(i!=1?1ll*(pwr(i,n+1)-1)*pwr(i-1,md-2)%md:n+1)%md*inv[i]%md;b[0]=1;
NTT(a,2*n+2,G);NTT(b,2*n+2,G);
for(int i=0;i<=4*n+4;i++)a[i]=1ll*a[i]*b[i]%md;
NTT(a,2*n+2,Gi);
int lim=0;while((1<<lim)<=2*n+2)lim++;lim=pwr(1<<lim,md-2);
for(int i=0;i<=n;i++)a[i]=1ll*a[i]*lim%md;
long long ans=0;
for(int i=0;i<=n;i++)ans=(ans+1ll*fac[i]*pwr(2,i)%md*a[i]%md)%md;
printf("%lld",ans);
return 0;
}
CF932E Team Work
求:
有:
点击查看代码
#include<bits/stdc++.h>
using namespace std;
const long long md=1e9+7;
int n,k;
long long S[5005][5005],ans;
inline long long pwr(long long x,long long y){
long long res=1;
while(y){
if(y&1)res=res*x%md;
x=x*x%md;y>>=1;
}return res;
}
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]*j+S[i-1][j-1])%md;
long long fac=1,cho=1;
for(int i=1;i<=k&&i<=n;i++){
cho=cho*(n-i+1)%md*pwr(i,md-2)%md;
fac=fac*i%md;
ans=(ans+fac*cho%md*pwr(2,n-i)%md*S[k][i]%md)%md;
}
printf("%lld",ans);
return 0;
}
斯特林反演
定义:
证明:
我们先证这个反转公式
反转公式1:
反转公式2:
由第一反演公式可得,对于任意列向量 \(f\),有:
令 \(A_{i,j}=∑_{t=i}^j\limits (−1)^{j−t}{j\brack t}{t\brace i}\)、\(B_{i,j}=∑_{t=i}^j\limits (−1)^{i−t}{j\brace t}{t\brack i}\) ,将上式写成矩乘形式,有:
显然 \(A=B=I\) ,因此两个反转公式得证。
进一步,令 \(X_{i,j}=(-1)^{i-j}{i\brack j}\)、\(Y_{i,j}={i\brace j}\),有:
因此 \(X\) 与 \(Y\) 互逆。
同样,令 \(X_{i,j}'=(-1)^{i-j}{i\brack j}\)、\(Y_{i,j}'={i\brace j}\),有:
可得 \(X'\) 与 \(Y'\) 互逆,进而可知斯特林反演公式成立。
例题:
#4671. 异或图
直接处理显然很难,考虑容斥。
\(f[i]\) 表示至少有 \(i\) 个联通块的方案,形如设立 \(i\) 个联通块轮廓,联通块内连边随意,联通块与联通块之间无连边。
\(g[i]\) 表示恰好有 \(i\) 个联通块的方案,形如设立 \(i\) 个联通块轮廓,在保证内部联通的情况下,外部块与块间无连边。
显然:
根据斯特林反演:
因此:
而 \({i\brack 1}\) 是阶乘形式:\({i\brack 1}=(i−1)!\)
化简答案为:\(g[1]=∑_{i=1}^n\limits (−1)^{i−1}(i−1)!f[i]\)
考虑 \(f[i]\) 如何求出:枚举点所属联通块状态,则我们使块与块之间无边,状压每个图的S表示点与点之间的连边(不属同一联通块),压到线性基里去,\(cnt\) 表示线性基元素,这些元素是不能选择的(相异),故答案为 \(2^{S-cnt}\)
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int S;
int n,vis[15];
char s[105];
bool ver[65][15][15];
long long bas[65],fac[15],ans;
void dfs(int pos=1,int lim=0){
if(pos>=n+1){
int cnt=0;
memset(bas,0,sizeof(bas));
for(int i=1;i<=S;i++){
int tot=0;long long s0=0;
for(int j=1;j<=n;j++){
for(int t=j+1;t<=n;t++){
if(vis[j]!=vis[t])s0|=(1ll*ver[i][j][t]<<tot),tot++;
}
}
for(int j=0;j<tot;j++){
if((s0>>j)&1){
if(!bas[j]){bas[j]=s0;cnt++;break;}
else s0^=bas[j];
}
}
}
ans+=(lim&1?1:-1)*fac[lim-1]*(1ll<<(S-cnt));
return ;
}
for(int i=1;i<=lim+1;i++){
vis[pos]=i;
dfs(pos+1,max(lim,i));
}
}
int main(){
scanf("%d",&S);
for(int i=1;i<=S;i++){
scanf("%s",s+1);
int m=strlen(s+1);while(n*(n-1)/2<m)n++;
int k=0;
for(int j=1;j<=n;j++){
for(int t=j+1;t<=n;t++)ver[i][j][t]=s[++k]-'0';
}
}
fac[0]=1;for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i;
dfs();
printf("%lld",ans);
return 0;
}