劼爷上的课现在才去整理…
第一类Stirling数
s(n,m)表示n个元素组成m个圆排列
由以上定义我们可以得出递推公式:
s(n,m)=∑n−1i=0s(n−i,m−1)∗Ci−1n−1∗(i−1)!
以及
s(n,m)=s(n−1,m−1)+s(n−1,m)∗(n−1)
那么我们可以得到
∏i=ni=1(x+i−1)=∑ni=0s(n,i)xi
然后可以我们就可以在
O(nlog2n)
内求出一行Stirling数
假设我们已经有了
F(x,n)=∏i=ni=1(x+i−1)
则
F(x,2n)=∏i=ni=1(x+i−1)∗∏i=ni=1(x+n+i−1)
=F(x,n)∗G(x,n)
其中
[xi]G(x,n)=∑nj=iCij∗[xj]F[x,n]∗nj−i
然后就成了一个卷积形式
FNT即可
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstring>
using namespace std;
#define ll long long
const
int Mod=998244353,p=3,D=Mod-1;
int Add(int a,int b)
{
a+=b;
if(a>=Mod)a-=Mod;
if(a<0)a+=Mod;
return a;
}
int Mul(int a,int b)
{
return a*1ll*b%Mod;
}
int Pow(int a,int x)
{
int res=1;
for(;x;x>>=1,a=a*1ll*a%Mod)
if(x&1)res=res*1ll*a%Mod;
return res;
}
int Rev[400001];
void Rader(int n)
{
for(int i=0;i<(1<<n);i++)Rev[i]=Rev[i>>1]>>1|((i&1)<<n-1);
}
int FNT(int *a,int n,int d)
{
Rader(n);
int len=1<<n;
for(int i=0;i<len;i++)if(Rev[i]>i)swap(a[Rev[i]],a[i]);
for(int i=1;i<len;i*=2)
{
int W=(d==1?Pow(p,D/2/i):Pow(p,D-D/2/i));
for(int j=0;j<len;j+=2*i)
{
int W0=1;
for(int k=0;k<i;k++)
{
int X=a[j+k],Y=Mul(W0,a[j+k+i]);
a[j+k]=Add(X,Y);a[j+k+i]=Add(X,-Y);
W0=Mul(W0,W);
}
}
}
int In=Pow(len,D-1);
if(d==-1)for(int i=0;i<len;i++)a[i]=Mul(a[i],In);
}
void MUL(int*a,int alen,int *b,int blen,int *M,int &Ml)
{
int n=1;
while((1<<n)<=alen+blen)n++;
FNT(a,n,1);FNT(b,n,1);
for(int i=0;i<(1<<n);i++)
M[i]=Mul(a[i],b[i]);
FNT(M,n,-1);
Ml=1<<n;
while(Ml&&!M[Ml-1])Ml--;
}
int Last[100001],Now[100001],L2,Last2[100001];
int F[100001],flen,G[100001],glen;
int Fact[100001],Inv[100001];
void Solve(int n)
{
if(n==1){Last[0]=0,Last[1]=1;return;}
int t=n>>1;
Solve(t);
flen=t;
memset(F,0,sizeof(F));
memset(G,0,sizeof(G));
memset(Last2,0,sizeof(Last2));
for(int i=0;i<=t;i++)
F[t-i]=Mul(Last[i],Fact[i]),Last2[i]=Last[i];
G[0]=1;
for(int i=1;i<=t;i++)
G[i]=Mul(Mul(G[i-1],Fact[i-1]),Mul(Inv[i],t));
MUL(F,t+1,G,t+1,Now,L2);
L2--;
for(int i=0;i<=t;i++)
Last[t-i]=Mul(Now[i],Inv[t-i]);
if(n&1)
{
for(int i=t+1;i;i--)
Last[i]=Add(Last[i-1],Mul(n-1,Last[i]));
Last[0]=Mul(n-1,Last[0]);
}
MUL(Last,t+2,Last2,n-t+1,Last,L2);
}
int main()
{
int n;
scanf("%d",&n);
Fact[0]=1;
for(int i=1;i<=n;i++)Fact[i]=Mul(Fact[i-1],i);
Inv[n]=Pow(Fact[n],Mod-2);
for(int i=n-1;~i;i--)
Inv[i]=Mul(Inv[i+1],i+1);
Solve(n);
int l,r;
for(int i=0;i<=n;i++)
printf("%d ",Last[i]);
return 0;
}
第二类Stirling数:
S(n,m)表示n个元素分成m个集合的方案数
那么由定义可得递推式:
S(n,m)=S(n,m−1)+S(n−1,m)∗m
并且
kn=∑km=0Amk∗S(n,m)=∑km=0m!∗Cmk∗S(n,m)
可以用n个球放入k个篮子,允许存在空篮子的方案数
由上式 我们经过二项式反演之后可以得到
S(n,m)=∑k=0(−1)m−k∗Ckm∗(m−k)n
S(n,m)=∑k=0(−1)k∗1k!∗1m−k!∗(m−k)n
同样可以经过FFT
(nlog2n)
得出
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
#define ll long long
const
int Mod=998244353,p=3,D=Mod-1;
int Add(int a,int b)
{
a+=b;
if(a>=Mod)a-=Mod;
if(a<0)a+=Mod;
return a;
}
int Mul(int a,int b)
{
return a*1ll*b%Mod;
}
int Pow(int a,int x)
{
int res=1;
for(;x;x>>=1,a=a*1ll*a%Mod)
if(x&1)res=res*1ll*a%Mod;
return res;
}
int Rev[100001];
void Rader(int n)
{
for(int i=0;i<(1<<n);i++)Rev[i]=(Rev[i>>1]>>1)|((i&1)<<n-1);
}
void FNT(int *a,int n,int d)
{
Rader(n);
int len=1<<n,Inv=Pow(len,Mod-2);
for(int i=0;i<len;i++)if(Rev[i]>i)swap(a[Rev[i]],a[i]);
for(int i=1;i<len;i*=2)
{
int W0=(d==1?(Pow(p,D/2/i)):(Pow(p,D-D/2/i)));
for(int j=0;j<len;j+=2*i)
{
int W=1;
for(int k=0;k<i;k++)
{
int X=a[k+j],Y=Mul(a[k+j+i],W);
a[k+j]=Add(X,Y);
a[k+i+j]=Add(X,-Y);
W=Mul(W,W0);
}
}
}
if(d==-1)
for(int i=0;i<len;i++)a[i]=Mul(a[i],Inv);
}
int F[100001],G[100001],Fact[100001],Inv[100001];
int Ans[100001];
int main()
{
int n;
scanf("%d",&n);
int B=1;
while((1<<B)<=n)B++;B++;
Fact[0]=1;
for(int i=1;i<=n;i++)Fact[i]=Mul(Fact[i-1],i);
Inv[n]=Pow(Fact[n],Mod-2);
for(int i=n-1;~i;i--)Inv[i]=Mul(Inv[i+1],i+1);
for(int i=0;i<=n;i++)
{
F[i]=Add(0,Mul((i&1?-1:1),Inv[i]));
G[i]=Mul(Pow(i,n),Inv[i]);
}
FNT(F,B,1);
FNT(G,B,1);
for(int i=0;i<(1<<B);i++)
Ans[i]=Mul(F[i],G[i]);
FNT(Ans,B,-1);
for(int i=0;i<=n;i++)
printf("%d ",Ans[i]);
return 0;
}
第二类Stirling数还有一个好性质:
xn=∑ni=0Aix∗S(n,i)
依旧考虑x个盒子放n个不同的球,可以存在空盒子
则左式可以看做枚举非空盒子个数 然后按顺序取出i个盒子放入球
根据这个性质我们可以很方便的处理一类代价为n次方的问题
HDU4625 JZPTREE
给出一棵n个节点边权为1的树,对每一个节点i,输出
∑nj=1dist(i,j)m
n≤5×104,m≤500
注意到此题m较小,可以用上述方法求解
唯一的问题就是如何从
Aix
得出
Aix+1
对于
Aix+1
我们有递推公式
Aix+1=Aix+i∗Ai−1x
然后树形DP处理出每一个点
u
的
时间复杂度:
O(n∗m)
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
char c;
inline void read(int&a)
{
a=0;do c=getchar();while(c<'0'||c>'9');
while(c<='9'&&c>='0')a=(a<<3)+(a<<1)+c-'0',c=getchar();
}
const
int Mod=10007;
int S[501][603];
int Low[60001][603];
int High[60001][603];
int n,m;
struct Chain
{
Chain*next;
int u;
}*Head[50002];
Chain R[100001];
int A;
inline void Add(int a,int b)
{Chain*tp=R+A;A++;tp->next=Head[a];Head[a]=tp;tp->u=b;}
int F[60001];
inline int Up(int u,int i)
{
return ((i>=1?i*1ll*(Low[u][i-1]):0)+(i>=0?Low[u][i]:0))%Mod;
}
inline int add(int a,int b)
{
a+=b;
if(a>=Mod)a-=Mod;
if(a<0)a+=Mod;
return a;
}
inline int mul(int a,int b)
{
a=a*b;
a%=Mod;
return a;
}
void DFS(int u,int f)
{
F[u]=f;
for(Chain*tp=Head[u];tp;tp=tp->next)
if(tp->u!=f)
{
DFS(tp->u,u);
for(int i=0;i<=m;i++)
Low[u][i+1]=add(Low[u][i+1],Up(tp->u,i+1));
Low[u][0]=add(Low[u][0],Low[tp->u][0]);
}
Low[u][0]=add(1,Low[u][0]);
}
void DFS2(int u,int f)
{
for(Chain*tp=Head[u];tp;tp=tp->next)
if(tp->u!=f)
{
for(int i=0;i<=m;i++)
High[tp->u][i]
=add(add(add(High[u][i],-Up(tp->u,i)),mul(i,i>=1?add(High[u][i-1],-Up(tp->u,i-1)):0)),Low[tp->u][i]);
DFS2(tp->u,u);
}
}
int Calc(int u)
{
int res=0,f=F[u];
if(u==1)
{
for(int i=1;i<=m;i++)
{
res=add(res,mul(S[m][i],Low[u][i]));
}
return res;
}
for(int i=1;i<=m;i++)
{
res=
add(res,
mul(S[m][i],
High[u][i]
));
}
return res;
}
int main()
{
freopen("self.in","r",stdin);
freopen("self.out","w",stdout);
S[0][0]=1;
for(int i=1;i<=500;i++)
for(int j=1;j<=500;j++)
S[i][j]=(S[i-1][j-1]+S[i-1][j]*j)%Mod;
int T;
read(T);
while(T--)
{
A=0;
memset(F,0,sizeof(F));
memset(Low,0,sizeof(Low));
memset(High,0,sizeof(High));
memset(Head,0,sizeof(Head));
read(n),read(m);
for(int i=1;i<n;i++)
{
int a,b;
read(a),read(b);
Add(a,b),Add(b,a);
}
DFS(1,1);
for(int i=0;i<=m;i++)High[1][i]=Low[1][i];
DFS2(1,1);
for(int i=1;i<=n;i++)
printf("%d\n",Calc(i));
}
}
Hackerrank Costly Graphs
题目大意:
求所有节点为n的图的权值和
这里图
G
的权值和定义为
D(u)
为节点u的度
m
为给定常数
1≤m≤2⋅105
显然每个节点对于答案贡献独立且相等 我们考虑其中任意节点
u
的贡献
则有贡献
W(u)=2C2n−1∗∑n−1i=0Cin−1∗∑mt=0Ati∗S(m,t)
W(u)=2C2n−1∗∑n−1i=0(n−1)!(n−1−i)!∗∑min(m,i)t=01(i−t)!∗S(m,t)
W(u)=(n−1)!∗2C2n−1∗∑mt=0S(m,t)∗∑n−1i=t1(n−1−i)!∗(i−t)!
W(u)=(n−1)!∗2C2n−1∗∑mt=0S(m,t)∗∑n−1−ti=01(n−1−t−i)!∗(i)!
W(u)=(n−1)!∗2C2n−1∗∑mt=0S(m,t)∗1(n−1−t)!∗2n−1−t
W(u)=2C2n−1∗∑mt=0S(m,t)∗Atn−1∗2n−1−t
答案就是
n∗W(u)
注意幂次的取模需要费马小定理..
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdlib>
using namespace std;
const
int Mod=1005060097,p=7;
char c;
inline void read(int&a){a=0;do c=getchar();while(c<'0'||c>'9');while(c<='9'&&c>='0')a=(a<<3)+(a<<1)+c-'0',c=getchar();}
int Pow(int a,int x)
{
int res=1;
for(;x;x>>=1,a=a*1ll*a%Mod)
if(x&1)res=res*1ll*a%Mod;
return res;
}
int mul(int a,int b){return a*1ll*b%Mod;}
int mul(int a,int b,int Mod){return a*1ll*b%Mod;}
int add(int a,int b)
{
a=a+b;
if(a>=Mod)a-=Mod;
if(a<0)a+=Mod;
return a;
}
int Rev[1000001];
inline void Rader(int n)
{
for(int i=0;i<(1<<n);i++)Rev[i]=(Rev[i>>1]>>1)|((i&1)<<n-1);
}
void FNT(int *a,int n,int l)
{
int len=1<<n,I=Pow(len,Mod-2);
Rader(n);
for(int i=0;i<len;i++)
if(Rev[i]>i)swap(a[Rev[i]],a[i]);
for(int i=1;i<len;i*=2)
{
int W=Pow(p,l==1?(Mod-1)/2/i:(Mod-1-(Mod-1)/2/i));
for(int j=0;j<len;j+=2*i)
{
int W0=1;
for(int k=0;k<i;k++)
{
int x=a[j+k],y=mul(a[j+k+i],W0);
a[j+k]=add(x,y);
a[j+k+i]=add(x,-y);
W0=mul(W0,W);
}
}
}
if(l==-1)
for(int i=0;i<len;i++)a[i]=mul(a[i],I);
}
int Fact[1000001],Inv[1000001];
int F[1000001],G[1000001],S[1000001];
int C2(int n)
{
if(n&1)return mul(n,n-1>>1,Mod-1);
return mul(n>>1,n-1,Mod-1);
}
int main()
{
int T;
Fact[0]=1;
for(int i=1;i<=200000;i++)Fact[i]=mul(Fact[i-1],i);
Inv[200000]=Pow(Fact[200000],Mod-2);
for(int i=199999;~i;i--)Inv[i]=mul(Inv[i+1],i+1);
read(T);
while(T--)
{
int n,m;
read(n),read(m);
for(int i=0;i<=m;i++)
F[i]=add(0,(i&1?-1:1)*Inv[i]);
for(int i=0;i<=m;i++)
G[i]=mul(Pow(i,m),Inv[i]);
int B=1;
while((1<<B)<=m+m+1)B++;
FNT(F,B,1);FNT(G,B,1);
for(int i=0;i<(1<<B);i++)
S[i]=mul(F[i],G[i]);
FNT(S,B,-1);
int ans=0,A=1,To=Pow(2,n-1),I=(Mod>>1)+1;
for(int i=0;i<=m;To=mul(To,I),A=mul(A,(n-1-(i++))))
ans=add(ans,mul(A,mul(To,S[i])));
ans=mul(ans,mul(Pow(2,C2(n-1)),n));
printf("%d\n",ans);
for(int i=0;i<(1<<B);i++)F[i]=G[i]=S[i]=0;
}
return 0;
}
同时Stirling数可以应用于幂级数..
再来看看小火车在今年NOI十连测上的一道题:
一张无向图的权值为联通块个数的m次方,问所有n个点的带标号的无向图的权值和。
答案对998244353取模。
T≤1000,n≤30000,m≤15。
构造一个多项式
L=∑ni=1xi
其中
xi
表示连通块
i
是否存在
则
由扩展二项式定理可得
项
∏xtiai
的系数为
m!∏ti!
,其中
∑ti=m
考虑所有由
xa1
至
xak
组成的项的系数和
对于任意
∑ki=1ti=m
我们有
∑m!∏ti!=∑S(m,k)∗k!
则对于某
k
个联通块,如果它们在图中同时出现,那么贡献为
设fi为i个点的连通图数目
fi=2C2i−∑i−1j=1Cj−1i−1∗fj∗2C2i−j
我们可以使用多项式求逆的方法求出
f
有
gi,j=∑ik=1Cki∗fk∗gi−k,j−1
可以用
卷积算出
答案即为
∑mj=1j!∗S(m,j)∑ni=1Cin∗gi,j∗2C2n−i
可以把
∑ni=1Cin∗gi,j∗2C2n−i
看成
T(j,n)
m次卷积预处理之后O(m)累加即可
时间复杂度
O(Tm+mnlog2n)
多项式求逆的话具体看2015年集训队鏼爷论文
自己的代码常数巨大
#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
using namespace std;
const
int Mod=998244353,p=3,D=Mod-1;
int Pow(int a,int x)
{
int res=1;
for(;x;x>>=1,a=a*1ll*a%Mod)
if(x&1)res=res*1ll*a%Mod;
return res;
}
inline
int Mul(int a,int b)
{return a*1ll*b%Mod;}
inline int Add(int a,int b)
{
a+=b;
if(a>=Mod)a-=Mod;
if(a<0)a+=Mod;
return a;
}
int Rev[2000001];
int Rader(int n)
{
for(int i=0;i<1<<n;i++)Rev[i]=(Rev[i>>1]>>1)|((i&1)<<n-1);
}
int FNT(int *a,int n,int f)
{
int len=1<<n,I=Pow(len,D-1);
Rader(n);
for(int i=0;i<len;i++)if(Rev[i]>i)swap(a[Rev[i]],a[i]);
for(int i=1;i<len;i*=2)
{
int W=Pow(p,f==1?D/2/i:(D-D/2/i));
for(int j=0;j<len;j+=i*2)
{
int W0=1;
for(int k=0;k<i;k++)
{
int x=a[j+k],y=Mul(W0,a[i+j+k]);
a[j+k]=Add(x,y),a[i+j+k]=Add(x,-y);
W0=Mul(W0,W);
}
}
}
if(f==-1)
for(int i=0;i<len;i++)a[i]=Mul(a[i],I);
}
int T[200001],A[200001];
int Fact[100001],Inv[100001];
int C(int x)
{
return (x*1ll*(x-1)>>1)%D;
}
int B[200001],Ca[200001];
void Find(int Len)
{
if(Len==1)
{B[0]=1;return;}
int t=1;
Find(Len>>1);
memset(Ca,0,sizeof(Ca));
memset(A,0,sizeof(A));
while((1<<t)<=3*Len)t++;
for(int i=0;i<Len;i++)
A[i]=T[i];
FNT(A,t,1);
FNT(B,t,1);
for(int i=0;i<1<<t;i++)
Ca[i]=Add(Mul(2,B[i]),-Mul(B[i],Mul(B[i],A[i])));
FNT(A,t,-1);
FNT(Ca,t,-1);
FNT(B,t,-1);
memset(B,0,sizeof(B));
for(int i=0;i<Len;i++)
B[i]=Ca[i];
}
int S[101][101];
int G[2][100001];
int F[100001];
int ANS[101][100001];
int BIT[100001];
int main()
{
freopen("self.in","r",stdin);
freopen("self.out","w",stdout);
int N=16384*2-1;
Fact[0]=1;
int P=116195171;
for(int i=1;i<=N;i++)Fact[i]=Mul(i,Fact[i-1]);
Inv[N]=Pow(Fact[N],D-1);
for(int i=N-1;~i;i--)Inv[i]=Mul(i+1,Inv[i+1]);
for(int i=1;i<=N;i++)
T[i]=Mul(Pow(2,C(i)),Inv[i]);
T[0]=1;
Find(N+1);
int TTT=16;
FNT(B,TTT,1);
memset(T,0,sizeof(T));
for(int i=1;i<=N;i++)
T[i]=Mul(Inv[i-1],Pow(2,C(i)));
FNT(T,TTT,1);
for(int i=0;i<1<<TTT;i++)
F[i]=Mul(T[i],B[i]);
FNT(F,TTT,-1);
B[0]=0;
int cas;
S[1][1]=1;
for(int i=2;i<=100;i++)
for(int j=1;j<=100;j++)
S[i][j]=Add(S[i-1][j-1],Mul(j,S[i-1][j]));
N=1;
int n=16384*2,m=15,now=1,next=0;
G[now][1]=1;
while((1<<N)<n)N++;
N++;
for(int i=1;i<=n;i++)
G[now][i]=Mul(F[i],Fact[i-1]);
memset(T,0,sizeof(T));
for(int i=1;i<=n;i++)
T[i]=F[i];
T[0]=0;
FNT(T,N,1);
for(int i=1;i<n;i++)
BIT[i]=Mul(Inv[i],Pow(2,C(i)));
BIT[0]=1;
FNT(BIT,N,1);
for(int j=1;j<=m;j++,next^=1,now^=1)
{
for(int i=1;i<=n;i++)
G[now][i]=Mul(G[now][i],Inv[i]);
FNT(G[now],N,1);
for(int i=0;i<1<<N;i++)
ANS[j][i]=Mul(BIT[i],G[now][i]);
FNT(ANS[j],N,-1);
memset(G[next],0,sizeof(G[next]));
for(int i=0;i<1<<N;i++)
G[next][i]=Mul(G[now][i],T[i]);
FNT(G[next],N,-1);
for(int i=n+1;i<1<<N;i++)
G[next][i]=0;
for(int i=1;i<=n;i++)
G[next][i]=Mul(G[next][i],Fact[i-1]);
G[next][0]=0;
}
scanf("%d",&cas);
while(cas--)
{
int ans=0;
int n,m;
memset(G,0,sizeof(G));scanf("%d%d",&n,&m);
for(int j=1;j<=m;j++)
ans=Add(ans,Mul(S[m][j],Mul(Fact[n],Mul(Fact[j],ANS[j][n]))));
printf("%d\n",ans);
}
return 0;
}
还有利用反演的:
fi=∑ij=0S(i,j)gi
为
gi=∑ij=0(−1)i−js(i,j)fj
的充要条件
然后容斥的这里就埋个坑...什么时候会做了再说吧
关于容斥:
满足任意两行或者任意两列都不相同的 n × m 的数字矩阵有多少个,
每一个格子内的数必须是 [1, C] 内的整数。对一个大质数
(109+7)
取模。
n,m,C≤4000
首先考虑如何容斥..
fi表示列上最多有i个等价类每行不重复的方案数
gi表示列上刚好有i个等价类每行不重复的方案数
则有
fi=S(m,i)∗Anmi=∑ii=1gjS(m,j)∗S(m,i)
然后就可以
O(m2)
暴力啦
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
using namespace std;
const
int Mod=1000000007;
#define ll long long
inline int Mul(int a,int b)
{
return a*1ll*b%Mod;
}
inline int Add(int a,int b)
{
a+=b;
if(a>=Mod)a-=Mod;
if(a<0)a+=Mod;
return a;
}
inline int Pow(int a,int x)
{
int res=1;
for(;x;x>>=1,a=a*1ll*a%Mod)
if(x&1)res=res*1ll*a%Mod;
return res;
}
int F[10001],G[10001];
int S[4001][4001];
class CountTables{
public:
int howMany(int n,int m,int C)
{
S[1][1]=1;
for(int i=2;i<=4000;i++)
for(int j=1;j<=4000;j++)S[i][j]=Add(S[i-1][j-1],Mul(j,S[i-1][j]));
int INV=1;
INV=Pow(INV,Mod-2);
int D=C,P=1;
for(int i=1;i<=m;i++)
{
int V=1;
for(int j=0;j<n;j++)
V=Mul(V,Add(D,-j));
F[i]=V;
D=Mul(D,C);
P=Mul(P,i);
}
for(int i=1;i<=m;i++)
{
G[i]=F[i];
for(int j=1;j<i;j++)
G[i]=Add(G[i],-Mul(G[j],S[i][j]));
}
return Mul(G[m],1);
return Mul(G[m],Pow(INV,Mod-2));
}}Gh;
int main()
{
int n,m,C;
scanf("%d%d%d",&n,&m,&C);
printf("%d\n",Gh.howMany(n,m,C));
}
To Be continued