斯特林数和贝尔数的求法
前置知识：
第一篇是基础多项式工业。
第二篇是关于各种数性质的推导，如果笔记里有式子这里就不会再证明了。
目录：
第二类斯特林数·行
[egin{Bmatrix}n\mend{Bmatrix}=dfrac{1}{m!}sum_{i=0}^{m}inom{m}{i}(-1)^{m-i}i^n\
=sum_{i=0}^{m}dfrac{i^n}{i!}dfrac{(-1)^{m-i}}{(m-i)!}
]
卷积即可。
注意一下这个毒瘤模数（原根是 (3) ），数组长度是 (2 imes 10^5) 不要开小。
void stirling2_row(int*g,int n){//注意！这个函数左闭右闭
static int A[M],B[M];
for(int i=0;i<=n;++i)A[i]=1ll*qpow(i,n)*math::ifc[i]%mod;
for(int i=0;i<=n;++i)B[i]=i&1?mod-math::ifc[i]:math::ifc[i];
poly_mul(A,B,A,n+1,n+1);
for(int i=0;i<=n;++i)g[i]=A[i];
}
第一类斯特林数·行
[x^{overline{n}} =sum_{i=0}^{n}egin{bmatrix}n\iend{bmatrix}x^i
]
发现等式右边就是第一类斯特林数的生成函数。
所以求出 (x^{overline{n}}) 即可。
不知道是谁想出来的，可以倍增搞。
[x^{overline{n+1}}=x^{overline{n}}(x+n)
]
这个可以 (O(n)) ，暴力扫一遍即可。
[x^{overline{2n}}=x^{overline{n}}(x+n)^{overline{n}}
]
有个技巧叫做多项式平移，不会可以看开头的笔记一。
所以我们可以直接由 (x^{overline{n}}) 得到 ((x+n)^{overline{n}}) ，直接平移即可。
注意一下这个毒瘤模数（原根是 (3) ），数组长度是 (262144) 不要开小。
void poly_shift(int*g,int*f,int n,int c){
static int A[M],B[M];
for(int i=0;i<n;++i)A[i]=1ll*math::fac[i]*f[i]%mod;
for(int i=0,j=1;i<n;++i,j=1ll*j*c%mod)B[i]=1ll*j*math::ifc[i]%mod;
reverse(B,B+n),poly_mul(A,B,A,n,n);
for(int i=0;i<n;++i)g[i]=1ll*A[n+i-1]*math::ifc[i]%mod;
}
void up_pow(int*g,int n){//注意！这个函数左闭右闭
static int A[M];
if(n==1)return g[0]=0,g[1]=1,void();
if(n&1){
up_pow(g,n-1);
for(int i=n;i>0;--i)g[i]=(g[i-1]+1ll*(n-1)*g[i]%mod)%mod;
}else{
up_pow(g,n/2);
clr(A,n/2+1),poly_shift(A,g,n/2+1,n/2),poly_mul(A,g,g,n/2+1,n/2+1);
}
}
void stirling1_row(int*g,int n){up_pow(g,n);}//注意！这个函数左闭右闭
第一类斯特林数·列
第一类斯特林数列的 ( m EGF) 是 (dfrac{(-ln(1-x))^k}{k!}) ，多项式快速幂即可。
注意 (-ln(1-x)=sumlimits_{i=1}dfrac{x^i}{i}) ，所以系数要左移一位不然没法快速幂。
注意一下这个毒瘤模数（原根是 (3) ），数组长度是 (131072) 不要开小。
void stirling1_column(int*g,int n,int k){//注意！这个函数左闭右闭
if(n<k){
for(int i=0;i<=n;++i)g[i]=0;
return;
}
static int A[M];
for(int i=0;i<=n;++i)A[i]=math::inv[i+1];
clr(g,n+1),poly_qpow(g,A,k,n+1);
for(int i=n;i>=k;--i)g[i]=g[i-k];
for(int i=0;i<k;++i)g[i]=0;
for(int i=k;i<=n;++i)g[i]=1ll*g[i]*math::fac[i]%mod*math::ifc[k]%mod;
}
第二类斯特林数·列
[egin{Bmatrix}n\mend{Bmatrix}=megin{Bmatrix}n-1\mend{Bmatrix}+egin{Bmatrix}n-1\m-1end{Bmatrix}\
]
设答案为 (F_{m}(x)=sumlimits_{i}egin{Bmatrix}i\mend{Bmatrix}x^i)
那么
[F_{m}(x)=mxF_{m}(x)+xF_{m-1}(x)\
F_{m}(x)=dfrac{xF_{m-1}(x)}{1-mx}
]
所以现在要求的就是
[dfrac{x^m}{prod_{i=1}^{m}(1-ix)}=dfrac{1}{prod_{i=1}^{m}(frac{1}{x}-i)}
]
把分母算出来求逆再系数平移就好了。
令 (G(x)=prodlimits_{i=1}^{m}(x-i)) ，那么 (G(dfrac{1}{x})) 就是把 (G(x)) 的系数翻转（reverse）的结果，而 (G(dfrac{1}{x})=prodlimits_{i=1}^{m}(dfrac{1}{x}-i)) ，(G(dfrac{1}{x})) 就是分母，也就是说我们求出 (G(x)) 就行了，显然 (G(x)=(x-1)^{underline{m}}=dfrac{x^{underline{m+1}}}{x}) 可以和上面一样倍增。
[x^{underline{2n}}=x^{underline{n}}(x-n)^{underline{n}}\
x^{underline{n+1}}=x^{underline{n}}(x-n)
]
void down_pow(int*g,int n){//注意！这个函数左闭右闭
static int A[M];
if(n==1)return g[0]=0,g[1]=1,void();
if(n&1){
down_pow(g,n-1);
for(int i=n;i>0;--i)g[i]=(g[i-1]+1ll*(mod-n+1)*g[i]%mod)%mod;
}else{
down_pow(g,n/2);
clr(A,n/2+1),poly_shift(A,g,n/2+1,mod-n/2),poly_mul(A,g,g,n/2+1,n/2+1);
}
}
void stirling2_column(int*g,int n,int k){//注意！这个函数左闭右闭
static int A[M];
if(n<k){
for(int i=0;i<=n;++i)g[i]=0;
return;
}
clr(A,k+2),down_pow(A,k+1);
for(int i=0;i<k+1;++i)A[i]=A[i+1];A[k+1]=0;
reverse(A,A+k+1);
clr(g,k+1),poly_inv(g,A,n-k+1);
for(int i=n;i>=k;--i)g[i]=g[i-k];
for(int i=0;i<k;++i)g[i]=0;
}
集合划分计数
贝尔数板子，就是让你求出贝尔数 (B_{1,cdots,100000})。
- 定义：(n) 个元素划分成任意个非空集合的方案数，就是一行第二类斯特林数的和。
- 递推：
[B_{n+1}=sum_{i=0}^{n}inom{n}{i}B_{i}
]
就是枚举有 (n-i) 个数和第 (n) 个数在同一个集合，有 (dbinom{n}{i}) 种方法选出这么多数。剩下数的划分方案就是 (B_i) 。
- 多项式科技
单个集合的 ( m EGF) 是 ({0,1,1,1,cdots}=e^x-1)
把它当作元素去生成集合，贝尔数的 ( m EGF) 就是 (exp(e^x-1)) 。
void bell(int*g,int n){
static int A[M];A[0]=0;
for(int i=1;i<n;++i)A[i]=math::ifc[i];
clr(g,n),poly_exp(g,A,n);
for(int i=0;i<n;++i)g[i]=1ll*g[i]*math::fac[i]%mod;
}
伯努利数
直接用( m EGF) ：(dfrac{x}{e^x-1}) ，求逆即可。
void bernoulli(int*g,int n){
static int A[M];
for(int i=0;i<n;++i)A[i]=ifc[i+1];
poly_inv(g,A,n);
for(int i=0;i<n;++i)g[i]=1ll*g[i]*fac[i]%mod;
}
欧拉数 ·行
用
[left<egin{matrix}n\mend{matrix}
ight>=sum_{k=0}^{m}inom{n+1}{k}(m+1-k)^n(-1)^{k}
]
设 (A_i=dbinom{n+1}{i}(-1)^{i},B_i=(i+1)^n) ，卷积即可。
void eula(int*g,int n){
static int A[M];
for(int i=0;i<n;++i)A[i]=qpow(i+1,n);
for(int i=0;i<n;++i)g[i]=i&1?mod-comb(n+1,i):comb(n+1,i);
poly_mul(g,A,g,n,n);
for(int i=n;i<n<<1;++i)g[i]=0;
}
完整代码
附送这部分的完整代码。poly的其他部分在开头《多项式笔记（一）》里面有，加到poly这个namespace最后就可以用了。
void stirling2_row(int*g,int n){//注意！这个函数左闭右闭
static int A[M],B[M];
for(int i=0;i<=n;++i)A[i]=1ll*qpow(i,n)*math::ifc[i]%mod;
for(int i=0;i<=n;++i)B[i]=i&1?mod-math::ifc[i]:math::ifc[i];
poly_mul(A,B,A,n+1,n+1);
for(int i=0;i<=n;++i)g[i]=A[i];
}
void poly_shift(int*g,int*f,int n,int c){
static int A[M],B[M];
for(int i=0;i<n;++i)A[i]=1ll*math::fac[i]*f[i]%mod;
for(int i=0,j=1;i<n;++i,j=1ll*j*c%mod)B[i]=1ll*j*math::ifc[i]%mod;
reverse(B,B+n),poly_mul(A,B,A,n,n);
for(int i=0;i<n;++i)g[i]=1ll*A[n+i-1]*math::ifc[i]%mod;
}
void up_pow(int*g,int n){//注意！这个函数左闭右闭
static int A[M];
if(n==1)return g[0]=0,g[1]=1,void();
if(n&1){
up_pow(g,n-1);
for(int i=n;i>0;--i)g[i]=(g[i-1]+1ll*(n-1)*g[i]%mod)%mod;
}else{
up_pow(g,n/2);
clr(A,n/2+1),poly_shift(A,g,n/2+1,n/2),poly_mul(A,g,g,n/2+1,n/2+1);
}
}
void stirling1_row(int*g,int n){up_pow(g,n);}//注意！这个函数左闭右闭
void stirling1_column(int*g,int n,int k){//注意！这个函数左闭右闭
if(n<k){
for(int i=0;i<=n;++i)g[i]=0;
return;
}
static int A[M];
for(int i=0;i<=n;++i)A[i]=math::inv[i+1];
clr(g,n+1),poly_qpow(g,A,k,n+1);
for(int i=n;i>=k;--i)g[i]=g[i-k];
for(int i=0;i<k;++i)g[i]=0;
for(int i=k;i<=n;++i)g[i]=1ll*g[i]*math::fac[i]%mod*math::ifc[k]%mod;
}
void down_pow(int*g,int n){//注意！这个函数左闭右闭
static int A[M];
if(n==1)return g[0]=0,g[1]=1,void();
if(n&1){
down_pow(g,n-1);
for(int i=n;i>0;--i)g[i]=(g[i-1]+1ll*(mod-n+1)*g[i]%mod)%mod;
}else{
down_pow(g,n/2);
clr(A,n/2+1),poly_shift(A,g,n/2+1,mod-n/2),poly_mul(A,g,g,n/2+1,n/2+1);
}
}
void stirling2_column(int*g,int n,int k){//注意！这个函数左闭右闭
static int A[M];
if(n<k){
for(int i=0;i<=n;++i)g[i]=0;
return;
}
clr(A,k+2),down_pow(A,k+1);
for(int i=0;i<k+1;++i)A[i]=A[i+1];A[k+1]=0;
reverse(A,A+k+1);
clr(g,k+1),poly_inv(g,A,n-k+1);
for(int i=n;i>=k;--i)g[i]=g[i-k];
for(int i=0;i<k;++i)g[i]=0;
}
void bell(int*g,int n){
static int A[M];A[0]=0;
for(int i=1;i<n;++i)A[i]=math::ifc[i];
clr(g,n),poly_exp(g,A,n);
for(int i=0;i<n;++i)g[i]=1ll*g[i]*math::fac[i]%mod;
}
void bernoulli(int*g,int n){
static int A[M];
for(int i=0;i<n;++i)A[i]=ifc[i+1];
poly_inv(g,A,n);
for(int i=0;i<n;++i)g[i]=1ll*g[i]*fac[i]%mod;
}
void eula(int*g,int n){
static int A[M];
for(int i=0;i<n;++i)A[i]=qpow(i+1,n);
for(int i=0;i<n;++i)g[i]=i&1?mod-comb(n+1,i):comb(n+1,i);
poly_mul(g,A,g,n,n);
for(int i=n;i<n<<1;++i)g[i]=0;
}