题解 P5396 【【模板】第二类斯特林数·列】

Great_Influence

2019-05-19 16:24:26

题解

我们根据第二类斯特林数的递推式来找方法。

\begin{Bmatrix}n\\ m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix}

我们设 S_n(x) 表示第二类斯特林数第 n 列的生成函数，即

S_m(x)=\sum_{i=0}^\infty \begin{Bmatrix}i\\ m\end{Bmatrix}x^i

那么，有:

S_m(x)=xS_{m-1}(x)+mxS_m

(1-mx)S_m(x)=xS_{m-1}(x)

S_m(x)=\frac{x}{1-mx}S_{m-1}(x)

S_m(x)=\frac{x^m}{\prod_{i=1}^m (1-ix)}

我们通过分治乘法算出下面，然后求逆后位移即可。

复杂度 O(n\log^2 n) 。

代码:

#include<cstdio>
#include<cstdlib>
#include<cctype>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cassert>
#include<iostream>
#define Rep(i,a,b) for(register int i=(a);i<=(b);++i)
#define Repe(i,a,b) for(register int i=(a);i>=(b);--i)
#define rep(i,a,b) for(register int i=(a);i<(b);++i)
#define pb push_back
#define mx(a,b) (a>b?a:b)
#define mn(a,b) (a<b?a:b)
typedef unsigned long long uint64;
typedef unsigned int uint32;
typedef long long ll;
using namespace std;

namespace IO
{
    const uint32 Buffsize=1<<15,Output=1<<24;
    static char Ch[Buffsize],*S=Ch,*T=Ch;
    inline char getc()
    {
        return((S==T)&&(T=(S=Ch)+fread(Ch,1,Buffsize,stdin),S==T)?0:*S++);
    }
    static char Out[Output],*nowps=Out;

    inline void flush(){fwrite(Out,1,nowps-Out,stdout);nowps=Out;}

    template<typename T>inline void read(T&x)
    {
        x=0;static char ch;T f=1;
        for(ch=getc();!isdigit(ch);ch=getc())if(ch=='-')f=-1;
        for(;isdigit(ch);ch=getc())x=x*10+(ch^48);
        x*=f;
    }

    template<typename T>inline void write(T x,char ch='\n')
    {
        if(!x)*nowps++='0';
        if(x<0)*nowps++='-',x=-x;
        static uint32 sta[111],tp;
        for(tp=0;x;x/=10)sta[++tp]=x%10;
        for(;tp;*nowps++=sta[tp--]^48);
        *nowps++=ch;
    }
}
using namespace IO;

void file()
{
#ifndef ONLINE_JUDGE
    FILE*DSD=freopen("water.in","r",stdin);
    FILE*CSC=freopen("water.out","w",stdout);
#endif
}

const int MAXN=1<<19;

namespace poly
{
    const int mod=167772161,gen=3;

    static int g[23][MAXN];

    inline int power(int u,int v)
    {
        register int sm=1;
        for(;v;v>>=1,u=(uint64)u*u%mod)if(v&1)
            sm=(uint64)sm*u%mod;
        return sm;
    }

    inline void predone()
    {
        Rep(i,1,21)
        {
            g[i][0]=1,g[i][1]=power(gen,(mod-1)>>i);
            Rep(j,2,4e5)g[i][j]=(ll)g[i][j-1]*g[i][1]%mod;
        }
    }

    static int Len,rev[MAXN];

    inline void calrev()
    {
        int II=log(Len)/log(2)-1;
        Rep(i,1,Len-1)rev[i]=rev[i>>1]>>1|(i&1)<<II;
    }

    inline int ad(int u,int v){return(u+=v)>=mod?u-mod:u;}

    inline void NTT(int*F,int typ)
    {
        Rep(i,1,Len-1)if(i<rev[i])swap(F[i],F[rev[i]]);
        for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
            for(register int j=0;j<Len;j+=i)rep(k,0,ii)
            {
                register int tt=(uint64)*(F+j+k+ii)*g[t][k]%mod;
                *(F+j+k+ii)=ad(*(F+j+k),mod-tt);
                *(F+j+k)=ad(*(F+j+k),tt);
            }
        if(typ==-1)
        {
            reverse(F+1,F+Len);
            register uint64 invn=power(Len,mod-2);
            rep(i,0,Len)*(F+i)=invn**(F+i)%mod;
        }
    }

    static int X[MAXN],Y[MAXN],Iv[MAXN];

    inline void mul(int *a,int *b,int *c,int lenl,int lenr)
    {
        if((ll)lenl*lenr<=100)
        {
            Rep(i,0,lenl+lenr)X[i]=0;
            Rep(i,0,lenl)Rep(j,0,lenr)
                X[i+j]=(X[i+j]+(ll)a[i]*b[j])%mod;
            Rep(i,0,lenl+lenr)c[i]=X[i];
            return;
        }
        for(Len=2;Len<=lenl+lenr;Len<<=1);
        calrev();
        Rep(i,0,lenl)X[i]=a[i];
        Rep(i,0,lenr)Y[i]=b[i];
        rep(i,lenl+1,Len)X[i]=0;
        rep(i,lenr+1,Len)Y[i]=0;
        NTT(X,1),NTT(Y,1);
        rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
        NTT(X,-1);
        Rep(i,0,lenl+lenr)c[i]=X[i];
        rep(i,lenl+lenr+1,Len)c[i]=0;
    }

    inline void Inv(int*F,int*G,int ln)
    {
        Iv[0]=power(F[0],mod-2);
        for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
        {
            rep(i,ln+1,Ln)F[i]=0;
            rep(i,0,Ln)X[i]=F[i],Y[i]=0;
            rep(i,0,(Ln>>1))Y[i]=Iv[i];
            Len=Ln,calrev();
            NTT(X,1),NTT(Y,1);
            rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
            NTT(X,-1);
            rep(i,0,(Ln>>1))X[i]=0;
            NTT(X,1);
            rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
            NTT(X,-1);
            rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
        }
        Rep(i,0,ln)G[i]=Iv[i];
    }
}
using poly::mul;
using poly::power;
using poly::Len;
using poly::calrev;
using poly::NTT;
using poly::mod;
using poly::predone;
using poly::Inv;
using poly::ad;

static int n,k,F[MAXN];

static int st[27][2][MAXN];

void getmul(int l,int r,int lev,int dr)
{
    if(l==r){st[lev][dr][0]=1,st[lev][dr][1]=mod-l;return;}
    int md=(l+r)>>1;
    getmul(l,md,lev+1,0),getmul(md+1,r,lev+1,1);
    mul(st[lev+1][0],st[lev+1][1],st[lev][dr],md-l+1,r-md);
}

int main()
{
    file();
    predone();
    read(n),read(k);
    if(k>n)Rep(i,0,n)write(0,' ');
    else
    {
        Rep(i,0,k-1)write(0,' ');
        getmul(1,k,0,0);
        Inv(st[0][0],F,n-k);
        Rep(i,0,n-k)write(F[i],' ');
    }
    flush();
    return 0;
}

还有另一个式子:

\sum _{n=k}^{\infty }\left\{{n \atop k}\right\}{\frac {x^{n}}{n!}}={\frac {(e^{x}-1)^{k}}{k!}}

(从维基百科上扣的)

我们可以通过 \ln 和 \exp 来做到 O(n\log n) 。不过一般的 \exp 跑得都比分治乘法慢。

代码:

#include<cstdio>
#include<cstdlib>
#include<cctype>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cassert>
#include<iostream>
#define Rep(i,a,b) for(register int i=(a);i<=(b);++i)
#define Repe(i,a,b) for(register int i=(a);i>=(b);--i)
#define rep(i,a,b) for(register int i=(a);i<(b);++i)
#define pb push_back
#define mx(a,b) (a>b?a:b)
#define mn(a,b) (a<b?a:b)
typedef unsigned long long uint64;
typedef unsigned int uint32;
typedef long long ll;
using namespace std;

namespace IO
{
    const uint32 Buffsize=1<<15,Output=1<<24;
    static char Ch[Buffsize],*S=Ch,*T=Ch;
    inline char getc()
    {
        return((S==T)&&(T=(S=Ch)+fread(Ch,1,Buffsize,stdin),S==T)?0:*S++);
    }
    static char Out[Output],*nowps=Out;

    inline void flush(){fwrite(Out,1,nowps-Out,stdout);nowps=Out;}

    template<typename T>inline void read(T&x)
    {
        x=0;static char ch;T f=1;
        for(ch=getc();!isdigit(ch);ch=getc())if(ch=='-')f=-1;
        for(;isdigit(ch);ch=getc())x=x*10+(ch^48);
        x*=f;
    }

    template<typename T>inline void write(T x,char ch='\n')
    {
        if(!x)*nowps++='0';
        if(x<0)*nowps++='-',x=-x;
        static uint32 sta[111],tp;
        for(tp=0;x;x/=10)sta[++tp]=x%10;
        for(;tp;*nowps++=sta[tp--]^48);
        *nowps++=ch;
    }
}
using namespace IO;

void file()
{
#ifndef ONLINE_JUDGE
    FILE*DSD=freopen("water.in","r",stdin);
    FILE*CSC=freopen("water.out","w",stdout);
#endif
}

const int MAXN=1<<19;

namespace poly
{
    const int mod=167772161,gen=3;

    static int g[23][MAXN],iv[MAXN];

    inline int power(int u,int v)
    {
        register int sm=1;
        for(;v;v>>=1,u=(uint64)u*u%mod)if(v&1)
            sm=(uint64)sm*u%mod;
        return sm;
    }

    inline void predone()
    {
        Rep(i,1,21)
        {
            g[i][0]=1,g[i][1]=power(gen,(mod-1)>>i);
            Rep(j,2,4e5)g[i][j]=(ll)g[i][j-1]*g[i][1]%mod;
        }
        iv[1]=1;
        Rep(i,2,4e5)iv[i]=mod-(uint64)mod/i*iv[mod%i]%mod;
    }

    static int Len,rev[MAXN];

    inline void calrev()
    {
        int II=log(Len)/log(2)-1;
        Rep(i,1,Len-1)rev[i]=rev[i>>1]>>1|(i&1)<<II;
    }

    inline int ad(int u,int v){return(u+=v)>=mod?u-mod:u;}

    inline void NTT(int*F,int typ)
    {
        Rep(i,1,Len-1)if(i<rev[i])swap(F[i],F[rev[i]]);
        for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
            for(register int j=0;j<Len;j+=i)rep(k,0,ii)
            {
                register int tt=(uint64)*(F+j+k+ii)*g[t][k]%mod;
                *(F+j+k+ii)=ad(*(F+j+k),mod-tt);
                *(F+j+k)=ad(*(F+j+k),tt);
            }
        if(typ==-1)
        {
            reverse(F+1,F+Len);
            register uint64 invn=power(Len,mod-2);
            rep(i,0,Len)*(F+i)=invn**(F+i)%mod;
        }
    }

    static int X[MAXN],Y[MAXN],Iv[MAXN];

    inline void mul(int *a,int *b,int *c,int lenl,int lenr)
    {
        if((ll)lenl*lenr<=100)
        {
            Rep(i,0,lenl+lenr)X[i]=0;
            Rep(i,0,lenl)Rep(j,0,lenr)
                X[i+j]=(X[i+j]+(ll)a[i]*b[j])%mod;
            Rep(i,0,lenl+lenr)c[i]=X[i];
            return;
        }
        for(Len=2;Len<=lenl+lenr;Len<<=1);
        calrev();
        Rep(i,0,lenl)X[i]=a[i];
        Rep(i,0,lenr)Y[i]=b[i];
        rep(i,lenl+1,Len)X[i]=0;
        rep(i,lenr+1,Len)Y[i]=0;
        NTT(X,1),NTT(Y,1);
        rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
        NTT(X,-1);
        Rep(i,0,lenl+lenr)c[i]=X[i];
        rep(i,lenl+lenr+1,Len)c[i]=0;
    }

    inline void Inv(int*F,int*G,int ln)
    {
        Iv[0]=power(F[0],mod-2);
        for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
        {
            rep(i,ln+1,Ln)F[i]=0;
            rep(i,0,Ln)X[i]=F[i],Y[i]=0;
            rep(i,0,(Ln>>1))Y[i]=Iv[i];
            Len=Ln,calrev();
            NTT(X,1),NTT(Y,1);
            rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
            NTT(X,-1);
            rep(i,0,(Ln>>1))X[i]=0;
            NTT(X,1);
            rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
            NTT(X,-1);
            rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
        }
        Rep(i,0,ln)G[i]=Iv[i];
    }

    static int ExX[MAXN],ExY[MAXN],Op[MAXN];

    inline void Deriv(int*F,int*G,int ln)
    {Rep(i,1,ln)G[i-1]=(uint64)F[i]*i%mod;G[ln]=0;}

    inline void Inter(int*F,int*G,int ln)
    {Repe(i,ln,1)G[i]=(uint64)F[i-1]*iv[i]%mod;G[0]=0;}

    static int LnX[MAXN];

    inline void Ln(int*F,int*G,int ln)
    {
        Deriv(F,LnX,ln),Inv(F,G,ln);
        for(Len=2;Len<=ln<<1;Len<<=1);
        rep(i,ln+1,Len)LnX[i]=G[i]=0;
        calrev();
        NTT(LnX,1),NTT(G,1);
        rep(i,0,Len)G[i]=(uint64)G[i]*LnX[i]%mod;
        NTT(G,-1);
        Inter(G,G,ln);
    }

    void cdq_Exp(int*a,int*F,int l,int r)
    {
        if(l==r)
        {
            if(!l)F[l]=1;
            else F[l]=(ll)F[l]*iv[l]%mod;
            return;
        }
        int md=(l+r)>>1;cdq_Exp(a,F,l,md);
        for(Len=2;Len<=r-l;Len<<=1);
        calrev();
        memset(X,0,sizeof(int)*Len);
        memset(Y,0,sizeof(int)*Len);
        memcpy(X,F+l,sizeof(int)*(md-l+1));
        memcpy(Y,a,sizeof(int)*(r-l));
        NTT(X,1),NTT(Y,1);
        rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
        NTT(X,-1);
        Rep(i,md+1,r)F[i]=ad(F[i],X[i-l-1]);
        cdq_Exp(a,F,md+1,r);
    }

    inline void Exp(int *F,int *G,int ln)
    {
        Rep(i,1,ln)Op[i-1]=(ll)F[i]*i%mod,Op[i]=0;
        memset(G,0,sizeof(int)*(ln+1));
        cdq_Exp(Op,G,0,ln);
        /*Op[0]=1;
        for(register int Lx=2;Lx>>1<=ln;Lx<<=1)
        {
            rep(i,Lx>>1,Lx)Op[i]=0;
            Ln(Op,ExX,Lx);
            rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
            rep(i,0,Lx>>1)ExY[i]=Op[i];
            rep(i,Lx>>1,Lx)ExX[i]=ExY[i]=0;
            Len=Lx;
            calrev();
            NTT(ExY,1),NTT(ExX,1);
            rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
            NTT(ExX,-1);
            rep(i,0,Lx>>1)Op[i+(Len>>1)]=ExX[i];
        }
        Rep(i,0,ln)G[i]=Op[i];*/
    }

    inline void Pow(int*F,int*G,int ln,ll k)
    {
        int lst=ln+1;
        Rep(i,0,ln)if(F[i]){lst=i;break;}
        if(ln&&(__int128)lst*k>ln)
        {memset(G,0,sizeof(int)*(ln+1));return;}
        int md=ln-k*lst,bs=F[lst],iv=power(bs,mod-2);
        Rep(i,lst,ln)G[i]=(ll)G[i]*iv%mod;
        Ln(G+lst,G,md);
        k%=mod;
        Rep(i,0,md)G[i]=(ll)G[i]*k%mod;
        Exp(G,G,md);
        bs=power(bs,k);
        Repe(i,md,0)G[i+lst*k]=(ll)G[i]*bs%mod;
        memset(G,0,sizeof(int)*(lst*k));
    }
}
using poly::mul;
using poly::power;
using poly::Len;
using poly::calrev;
using poly::NTT;
using poly::mod;
using poly::predone;
using poly::Inv;
using poly::Inter;
using poly::Deriv;
using poly::Ln;
using poly::Exp;
using poly::Pow;
using poly::ad;

static int n,k,X[MAXN],fac[MAXN],inv[MAXN];

int main()
{
    file();
    read(n),read(k);
    if(k>n)Rep(i,0,n)write(0,' ');
    else
    {
        predone();
        fac[0]=1;
        Rep(i,1,n)fac[i]=(ll)fac[i-1]*i%mod;
        inv[n]=power(fac[n],mod-2);
        Repe(i,n,1)inv[i-1]=(ll)inv[i]*i%mod;
        Rep(i,1,n)X[i]=inv[i];
        Pow(X,X,n,k);
        Rep(i,0,n)X[i]=(ll)X[i]*inv[k]%mod;
        Rep(i,0,n)write((ll)X[i]*fac[i]%mod,' ');
    }
    flush();
    return 0;
}