斯特林数和贝尔数的求法
前置知识:
第一篇是基础多项式工业。
第二篇是关于各种数性质的推导,如果笔记里有式子这里就不会再证明了。
目录:
第二类斯特林数·行
卷积即可。
注意一下这个毒瘤模数(原根是 (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}}) 即可。
不知道是谁想出来的,可以倍增搞。
这个可以 (O(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;
}
第二类斯特林数·列
设答案为 (F_{m}(x)=sumlimits_{i}egin{Bmatrix}i\mend{Bmatrix}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}) 可以和上面一样倍增。
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) 个元素划分成任意个非空集合的方案数,就是一行第二类斯特林数的和。
- 递推:
就是枚举有 (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;
}
欧拉数 ·行
用
设 (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;
}