斯特林数与贝尔数
https://blog.csdn.net/qq_33229466/article/details/75042895
https://www.cnblogs.com/gzy-cjoier/p/8426987.html
http://www.cnblogs.com/zhouzhendong/p/Stirling-Number.html
题解:
这个东西挺容易忘的。。
第一类斯特林数
定义:n个不同小球排成m个循环排列(不同排列交换顺序是相同的)的方案数(每个都不为空)
暴力算的代码:
for (int i=1;i<=n;++i) for (int j=1;j<=m;++j) for (register int k=1;k<=i;++k) f[i][j]=(f[i][j]+f[i-k][j-1]*inv[k])%mod;
答案是$\frac{f[n][m]*n!}{m!}$
原理就是将求出来的方案数除以每个块的元素个数(圆排列的性质)
考虑用组合意义来计算
用$S(n,m)$表示第一类斯特林数,我们去考虑最后一个放进去的元素,要么新开一个排列,要么加入之前的排列,所以
$$S(n,m)=(n-1)*S(n-1,m)+S(n-1,m-1)$$
#include <bits/stdc++.h> using namespace std; #define rint register int #define IL inline #define rep(i,h,t) for(int i=h;i<=t;i++) #define dep(i,t,h) for(int i=t;i>=h;i--) #define ll long long #define me(x) memset(x,0,sizeof(x)) #define mep(x,y) memcpy(x,y,sizeof(y)) #define mid ((h+t)>>1) namespace IO{ char ss[1<<24],*A=ss,*B=ss; IL char gc() { return A==B&&(B=(A=ss)+fread(ss,1,1<<24,stdin),A==B)?EOF:*A++; } template<class T> void read(T &x) { rint f=1,c; while (c=gc(),c<48||c>57) if (c=='-') f=-1; x=(c^48); while (c=gc(),c>47&&c<58) x=(x<<3)+(x<<1)+(c^48); x*=f; } char sr[1<<24],z[20]; ll Z,C1=-1; template<class T>void wer(T x) { if (x<0) sr[++C1]='-',x=-x; while (z[++Z]=x%10+48,x/=10); while (sr[++C1]=z[Z],--Z); } IL void wer1() { sr[++C1]=' '; } IL void wer2() { sr[++C1]='\n'; } template<class T>IL void maxa(T &x,T y) {if (x<y) x=y;} template<class T>IL void mina(T &x,T y) {if (x>y) x=y;} template<class T>IL T MAX(T x,T y){return x>y?x:y;} template<class T>IL T MIN(T x,T y){return x<y?x:y;} }; using namespace IO; int f[100][100],n,m; const int mo=998244353; int main() { freopen("1.in","r",stdin); freopen("2.out","w",stdout); read(n); read(m); f[0][0]=1; rep(i,1,90) rep(j,1,90) f[i][j]=(f[i-1][j-1]+1ll*(i-1)*f[i-1][j])%mo; cout<<f[n][m]<<endl; return 0; }
然后呢这个递推式依旧只能做成上面那个复杂度的
另一种形式
这个是用来化简一些式子使得能上fft之类用的
$A(n,m)=n*(n-1)*...*(n-m+1)$ (也就是所说的下降幂)
首先有结论$$A(n,m)=\sum \limits_{i=0}^{n} {(-1)}^{m-i)}*n^i*S(m,i)$$
考虑证明一下这个东西
它前面的系数等价于要在$0..(m-1)$之中选择$(m-i)$个的乘积和
我们看$S(n,m)$的递推式子,m每增加1说明少取了1个数,并且所取的数都不同,于是证毕
按照这个套路我们可以知道
x的上升幂$$x(x+1)...(x+n)=\sum\limits_{i=0}^{n} x^i*S(n,i)$$
更快的计算方法
我们总不能用大常数的线性递推来解决这个问题
我们使用它的上升幂
那么只需要求出$x(x+1)...(x+n-1)$的乘积就可以了
比较一般的方法是分治$fft$,$nlog^2{n}$
有一种更快的方法
我们考虑如果知道了$x$的$n$次上升幂,$x+n$的$n$次上升幂
这两个乘起来得到的就是$x$的$2n$次上升幂
所以我们要的就是从$f(x)=x$的$n$次上升幂得到$g(x)=x+n$的$n$次上升幂
$$f(x)=\sum\limits_{i=0}^{n} {a[i]*x^i}$$
$$g(x)=\sum\limits_{i=0}^{n} {a[i]*(x+n)^i}$$
将g的组合数展开,得到
那么$$g(x)=\sum\limits_{j=0}^{n} {\frac{x^j}{j!} \sum\limits_{i=j}^{n} {a[i]*i!*\frac{n^{i-j}}{(i-j)!}}}$$
我们令$f1(i)=a[n-i]*(n-i)!,f2(i)=\frac{n^i}{i!}$即可 那么$ans[n-i]=g[i]$
$$T(n)=T(n/2)+nlogn$$
总复杂度$nlogn$
#include <bits/stdc++.h> using namespace std; #define rint register int #define IL inline #define rep(i,h,t) for(int i=h;i<=t;i++) #define dep(i,t,h) for(int i=t;i>=h;i--) #define ll long long #define me(x) memset(x,0,sizeof(x)) #define mep(x,y) memcpy(x,y,sizeof(y)) #define mid ((h+t)>>1) namespace IO{ char ss[1<<24],*A=ss,*B=ss; IL char gc() { return A==B&&(B=(A=ss)+fread(ss,1,1<<24,stdin),A==B)?EOF:*A++; } template<class T> void read(T &x) { rint f=1,c; while (c=gc(),c<48||c>57) if (c=='-') f=-1; x=(c^48); while (c=gc(),c>47&&c<58) x=(x<<3)+(x<<1)+(c^48); x*=f; } char sr[1<<24],z[20]; ll Z,C1=-1; template<class T>void wer(T x) { if (x<0) sr[++C1]='-',x=-x; while (z[++Z]=x%10+48,x/=10); while (sr[++C1]=z[Z],--Z); } IL void wer1() { sr[++C1]=' '; } IL void wer2() { sr[++C1]='\n'; } template<class T>IL void maxa(T &x,T y) {if (x<y) x=y;} template<class T>IL void mina(T &x,T y) {if (x>y) x=y;} template<class T>IL T MAX(T x,T y){return x>y?x:y;} template<class T>IL T MIN(T x,T y){return x<y?x:y;} }; using namespace IO; const int N=2.1e5; const int mo=998244353; const int G=3; int n,m,l,r[N],w[N],a[N],b[N],k1[N],k2[N],ans[N],jc[N],inv[N]; IL int fsp(int x,int y) { ll now=1; while (y) { if (y&1) now=now*x%mo; x=(1ll*x*x)%mo; y>>=1; } return now; } IL void ntt_init() { l=0; for (n=1;n<=m;n<<=1) l++; for (int i=0;i<n;i++) r[i]=(r[i/2]/2)|((i&1)<<(l-1)); } IL void clear() { rep(i,0,n) a[i]=b[i]=0; } IL void ntt(int *a,int o) { rep(i,0,n-1) if (i>r[i]) swap(a[i],a[r[i]]); for (int i=1;i<n;i*=2) { int wn=fsp(G,(mo-1)/(i*2)); w[0]=1; rep(j,1,i-1) w[j]=1ll*w[j-1]*wn%mo; for (int j=0;j<n;j+=(i*2)) for (int k=0;k<i;k++) { int x=a[j+k],y=1ll*a[i+j+k]*w[k]%mo; a[j+k]=(x+y)%mo,a[i+j+k]=(x-y)%mo; } } if (o==-1) { reverse(&a[1],&a[n]); int inv=fsp(n,mo-2); for (int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mo; } } IL void get_cj(int *A,int *B,int len) { m=len*2; ntt_init(); rep(i,0,len-1) a[i]=A[i],b[i]=B[i]; ntt(a,1);ntt(b,1); rep(i,0,n-1) a[i]=1ll*a[i]*b[i]%mo; ntt(a,-1); rep(i,0,m) B[i]=a[i]; clear(); } void js(int *A,int k) { if (k==1) { A[1]=1; return; } int now=k/2; js(A,now); rep(i,0,now) k1[i]=1ll*jc[now-i]*A[now-i]%mo,k2[i]=1ll*fsp(now,i)*inv[i]%mo; get_cj(k1,k2,now+1); rep(i,0,now) k1[i]=1ll*k2[now-i]*inv[i]%mo; get_cj(A,k1,now+1); if (k%2==1) { rep(i,1,k) A[i]=k1[i-1]; rep(i,0,k-1) (A[i]+=1ll*k1[i]*(k-1)%mo)%=mo; } else { rep(i,0,k) A[i]=k1[i]; } } int main() { freopen("1.in","r",stdin); freopen("1.out","w",stdout); int m1; read(n); read(m1); jc[0]=1; rep(i,1,n) jc[i]=(1ll*jc[i-1]*i)%mo; inv[n]=fsp(jc[n],mo-2); inv[0]=1; dep(i,n-1,1) inv[i]=1ll*inv[i+1]*(i+1)%mo; js(ans,n); cout<<(ans[m1]+mo)%mo<<endl; return 0; }
第二类斯特林数
定义:n个不同小球放进m个相同盒子的方案数(每个都不为空)
暴力算的代码:
for (int i=1;i<=n;++i) for (int j=1;j<=m;++j) for (register int k=1;k<=i;++k) f[i][j]=(f[i][j]+f[i-k][j-1]*inv[jc[k]]%mod;
答案是$\frac{f[n][m]*n!}{m!}$
原理就是将求出来的方案数除以每个块内元素的排序
同样考虑通过其组合意义来计算
$$S(n,m)=S(n-1,m-1)+m*S(n-1,m) \tag1$$
一个性质
$$x^k=\sum\limits_{i=0}^{k} {S(k,i)*i!*C(x,i)}$$
其中$i!*C(x,i)$也是x的i次下降幂
那么如果证明这个呢,考虑它的意义
左边是k个不同小球放入x个不同盒子(可空)
右边枚举不空的盒子的数量
这个可以用来化简一些式子
快速计算第二类斯特林数,这个比较显然用容斥来计算,我们把盒子先当成不同的来处理
$$S(n,m)=\frac{\sum\limits_{i=0}^{m} {{(-1)}^{i}*C(m,i)*(m-i)^k}}{m!}$$
这个要计算一项是可以暴力算的,计算n相同时的多项可以fft
贝尔数
贝尔数的定义为大小为k的集合的划分方式