第二类斯特林数：列

第一、二类斯特林数算法推导与计算

最新推荐文章于 2020-06-25 11:45:17 发布

Freopen

最新推荐文章于 2020-06-25 11:45:17 发布

阅读量572

点赞数

CC 4.0 BY-SA版权

分类专栏：模板数论 FFT

本文链接：https://blog.youkuaiyun.com/qq_35950004/article/details/95650995

模板同时被 3 个专栏收录

212 篇文章

订阅专栏

数论

157 篇文章

订阅专栏

FFT

37 篇文章

订阅专栏

博客围绕第一类和第二类斯特林数展开，指出代码与“第一类斯特林数：行”极为相似，通过一系列推导得出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!}$
易得： $S_k(x) = kxS_k(x) + S_{k-1}(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(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^k$

又因为 $∏i=1k(1/x−i)∗xk\prod_{i=1}^k(1/x - i) * x^k$ 就是 $∏i=1k(x−i)\prod_{i=1}^k(x-i)$ 按原点反转，再往正方向移k位，即对称翻转即可。
那么我们就需要算 $∏i=1k(x−i)\prod_{i=1}^k (x-i)$
这个和"第一类斯特林数：列" 中 $∏i=0k(x+i)\prod_{i=0}^k (x+i)$
一个算法。
设 $fk(x)=∏i=1k(x−i)=∑i=0kaixif_k(x) =\prod_{i=1}^k(x-i) = \sum_{i=0}^k a_ix^i$
则 $f_2k(x) = f_k(x) f_k(x-k)$
考虑由 $f_k(x)$ 计算 $f_k(x-k)$
$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)!}$
后面把任意一个翻转一下就是个卷积。

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

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]);
}