第二类斯特林数:列

第一、二类斯特林数算法推导与计算
博客围绕第一类和第二类斯特林数展开,指出代码与“第一类斯特林数:行”极为相似,通过一系列推导得出Sk(x)的表达式,还探讨了计算∏i=1k(x−i)的方法,其与“第一类斯特林数:列”相关,最后可通过递归多项式求逆的倍增算法计算,还给出了AC代码。

代码和 “第一类斯特林数:行” 极为相似(具体数学:第一类斯特林数的三角形往负数走就是第二类斯特林数)

Sk(x)=∑i=0{ik}xii!S_k(x) = \sum_{i=0} \begin{Bmatrix} i \\ k\end{Bmatrix}\frac {x^i}{i!}Sk(x)=i=0{ik}i!xi
易得:Sk(x)=kxSk(x)+Sk−1(x)xS_k(x) = kxS_k(x) + S_{k-1}(x)xSk(x)=kxSk(x)+Sk1(x)x
那么:Sk(x)=xk∏i=1k(1−ix)S_k(x) = \frac {x^k}{\prod_{i=1}^k (1-ix)}Sk(x)=i=1k(1ix)xk

Sk(x)=∏i=1k(1−ix)Sk(1/x)=∏i=1k(1−i/x)=∏i=1k(x−i)/xkSk(x)=∏i=1k(1/x−i)∗xkS_k(x) = \prod_{i=1}^k(1-ix)\\ S_k(1/x) = \prod_{i=1}^k(1-i/x) = \prod_{i=1}^k(x-i) / x^k\\ S_k(x) =\prod_{i=1}^k(1/x - i) * x^kSk(x)=i=1k(1ix)Sk(1/x)=i=1k(1i/x)=i=1k(xi)/xkSk(x)=i=1k(1/xi)xk

又因为∏i=1k(1/x−i)∗xk\prod_{i=1}^k(1/x - i) * x^ki=1k(1/xi)xk就是∏i=1k(x−i)\prod_{i=1}^k(x-i)i=1k(xi)按原点反转,再往正方向移k位,即对称翻转即可。
那么我们就需要算∏i=1k(x−i)\prod_{i=1}^k (x-i)i=1k(xi)
这个和"第一类斯特林数:列" 中∏i=0k(x+i)\prod_{i=0}^k (x+i)i=0k(x+i)
一个算法。
fk(x)=∏i=1k(x−i)=∑i=0kaixif_k(x) =\prod_{i=1}^k(x-i) = \sum_{i=0}^k a_ix^ifk(x)=i=1k(xi)=i=0kaixi
f2k(x)=fk(x)fk(x−k)f_2k(x) = f_k(x) f_k(x-k)f2k(x)=fk(x)fk(xk)
考虑由fk(x)f_k(x)fk(x)计算fk(x−k)f_k(x-k)fk(xk)
fk(x−k)=∑i=0kai(x−k)i=∑i=0kai∑j=0i(ij)xj(−k)i−j=∑j=0kxjj!∑i=jkaii!∗(−k)i−j(i−j)!f_k(x-k) = \sum_{i=0}^k a_i(x-k)^i = \sum_{i=0}^ka_i\sum_{j=0}^i\binom{i}{j}x^j(-k)^{i-j}=\sum_{j=0}^k\frac {x_j}{j!}\sum_{i=j}^ka_ii!*\frac{(-k)^{i-j}}{(i-j)!}fk(xk)=i=0kai(xk)i=i=0kaij=0i(ji)xj(k)ij=j=0kj!xji=jkaii!(ij)!(k)ij
后面把任意一个翻转一下就是个卷积。

然后就可以像递归多项式求逆一样倍增的算了。
(我这个老年选手在"第一类斯特林数:列"上加了几行居然惊奇的过了。。。。)

AC Code:

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define maxn 600005
#define mod 167772161
using namespace std;

int wlen,w[maxn]={1},lg[maxn],r[maxn];
int Pow(int base,int k){
    int ret=1;
    for(;k;k>>=1,base=1ll*base*base%mod)
        if(k&1)
            ret=1ll*ret*base%mod;
    return ret;
}
void Init(int n){
    for(wlen=1;n>=2*wlen;wlen<<=1);
    for(int i=1,pw=Pow(3,(mod-1)/(2*wlen));i<=2*wlen;i++)
        w[i] = w[i-1] * 1ll * pw % mod;
    for(int i=2;i<=2*wlen;i++) 
        lg[i] = lg[i>>1] + 1;
}

void NTT(int *A,int n,int tp){
    int lgn = lg[n];
    for(int i=1;i<n;i++) r[i] = (r[i>>1]>>1)|((i&1)<<(lgn-1));
    for(int i=1;i<n;i++) if(i<r[i]) swap(A[i],A[r[i]]);
    for(int L=2;L<=n;L<<=1)
        for(int st=0,l=L>>1,inc=wlen/l;st<n;st+=L)
            for(int k=st,x=0;k<st+l;k++,x+=inc){
                int t=  1ll * (tp == 1 ? w[x] : w[2*wlen-x]) * A[k+l] % mod;
                A[k+l] = (A[k] - t) % mod , A[k] = (A[k] + t) % mod;
            }
    if(tp == -1){
        for(int i=0,inv=Pow(n,mod-2);i<n;i++)
            A[i] = 1ll * A[i] * inv % mod;
    }
}

void Inv(int *A,int *B,int n){
    B[0] = Pow(A[0],mod-2),B[1]=0;
    static int tmp[maxn];
    for(int k=2;k<(n<<1);k<<=1){
        for(int i=0;i<(k<<1);i++) tmp[i]=(i<k?A[i]:0),B[i]=i<k?B[i]:0;
        NTT(tmp,k<<1,1),NTT(B,k<<1,1);
        for(int i=0;i<(k<<1);i++) B[i] = 1ll * B[i] * (2 - 1ll * B[i] * tmp[i] % mod) % mod;
        NTT(B,k<<1,-1);
        for(int i=min(n,k);i<=(k<<1);i++) B[i]=0;
    }
}

int n,k,fac[maxn]={1,1},invf[maxn]={1,1},A[maxn],B[maxn];

void Solve(int *A,int n){
    if(!n){
        A[0] = 1;
        return;
    }
    if(n&1){
        Solve(A,n-1);
        for(int i=n;i>=0;i--)
            A[i] = ((i?A[i-1]:0)-1ll*A[i]*n) % mod;
        return;
    }
    static int tmp[3][maxn];
    int len = (1<<lg[n]+1);
    Solve(A,n/2);
    for(int i=0;i<len;i++) 
        tmp[0][i] = (i<=n/2 ? 1ll * A[i] * fac[i] % mod : 0),
        tmp[1][i] = (i<=n/2 ? (i?1ll * tmp[1][i-1] * (-n/2) : 1) % mod: 0),
        tmp[2][i] = 0;
    for(int i=0;i<=n/2;i++)
        tmp[2][n/2-i] = 1ll * tmp[1][i] * invf[i] % mod;
    NTT(tmp[0],len,1),NTT(tmp[2],len,1);
    for(int i=0;i<len;i++)
        tmp[0][i] = 1ll * tmp[0][i] * tmp[2][i] % mod;
    NTT(tmp[0],len,-1);
    for(int i=0;i<len;i++)
        tmp[0][i] = (i<=n/2 ? 1ll * tmp[0][i+n/2] * invf[i] % mod : 0),
        A[i] = (i<=n/2 ? A[i] : 0);
    NTT(tmp[0],len,1),NTT(A,len,1);
    for(int i=0;i<len;i++)
        A[i] = 1ll * A[i] * tmp[0][i] % mod;
    NTT(A,len,-1);
    for(int i=n+1;i<=len;i++) A[i] = 0; 
}

int main(){
    scanf("%d%d",&n,&k); 
    
    for(int i=2;i<=k;i++) fac[i] = 1ll * fac[i-1] * i % mod;
    invf[k] = Pow(fac[k] , mod-2);
    for(int i=k-1;i>=2;i--) invf[i] = 1ll * invf[i+1] * (i+1) % mod;
    
    Init(n*2);
    Solve(A,k);
    for(int i=0;i<k-i;i++) swap(A[i],A[k-i]);
    
    Inv(A,B,n+1);
    
    for(int i=0;i<=n;i++)
        printf("%d%c",i>=k?(B[i-k]+mod)%mod:0," \n"[i==n]);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值