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

· · 题解

本题小常数倍增NTTO(nlogn)解法!

虽然本蒟蒻人傻常数大跑了800+ms

但是作为一个求逆模板没有进入最优解前1000的人,这已经很快乐了

安利斯特林数小结,欢迎各位大佬前来吊打。

没有什么优美的公式帮忙了……我们来看看递推公式吧:

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

我们还是按套路,把答案设成多项式H_k=\sum\limits_{i=0}\begin{Bmatrix} i \\k \end{Bmatrix}x^i。

然后,我们发现H_k=\sum\limits_{i=0}\left(\begin{Bmatrix}i-1 \\ k-1\end{Bmatrix}+k \begin{Bmatrix}i-1 \\ k\end{Bmatrix}\right)x^i

=\sum\limits_{i=0}\begin{Bmatrix}i-1 \\ k-1\end{Bmatrix}x^i+k\sum\limits_{i=0}\begin{Bmatrix}i-1 \\ k\end{Bmatrix}x^i =x\sum\limits_{i=0}\begin{Bmatrix}i \\ k-1\end{Bmatrix}x^i+kx\sum\limits_{i=0}\begin{Bmatrix}i \\ k\end{Bmatrix}x^i=xH_{k-1}+kxH_k

总结一下就是H_k=xH_k+kxH_{k-1},那么H_k=\dfrac{x}{1-kx}H_{k-1}

根据递推公式能得到边界就是H_0=1,这里不再赘述。

那么,H_k=H_0\prod\limits_{i=1}^k\dfrac{x}{1-ix}=x^k\left(\prod\limits_{i=1}^k(1-ix)\right)^{-1}

重点就是怎么算\prod\limits_{i=1}^k(1-ix)。

\prod\limits_{i=1}^k(1-ix)=\prod\limits_{i=1}^kx(\dfrac{1}{x}-i)=x^k\prod\limits_{i=1}^k(x^{-1}-i)

使用一些小技巧,\prod\limits_{i=1}^k(x-i)翻转过来之后就是x^k\prod\limits_{i=1}^k(x^{-1}-i)

那么注意到\prod\limits_{i=1}^k(x-i)=\dfrac{x^{\underline{k+1}}}{x}我们只要求出x^{\underline{k+1}}就好。

咋整呢?我们还是考虑像这样来倍增(Orz rqy!):

我们知道x^{\underline{2n}}=x^{\underline n}*(x-n)^{\underline n}

我们知道一个多项式F(x)可以快速求出F(x-c),方法:

首先c=mod-c,那么问题变成F(x+c)

设F(x)的系数为f[0...n]

$=\sum\limits_{i=0}^nf[i]\sum\limits_{j=0}^i\dbinom{i}{j}x^jc^{i-j}$(二项式定理) $=\sum\limits_{i=0}^nf[i]\sum\limits_{j=0}^i\dfrac{i!}{j!(i-j)!}x^jc^{i-j} =\sum\limits_{i=0}^nf[i]i!\sum\limits_{j=0}^i\dfrac{x^j}{j!}\dfrac{c^{i-j}}{(i-j)!} 我们要求的就是$g[j]*j!=\sum\limits_{i=j}^nf[i]i!\dfrac{c^{i-j}}{(i-j)!}

我们设p[n]=f[n]n!;\ h[n]=\dfrac{c^{n}}{n!},则g[j]j!=\sum\limits_{i=j}^np[i]h[i-j]

但是我们发现i+i-j≠定值 , 我们把p翻转,得到g[j]j!=\sum\limits_{i=j}^np'[n-i]h[i-j]

那么n-i+i-j=n-j这就和i无关了。

卷起来之后翻转即可,记得最后乘上j!。

复杂度是T(n)=T(n/2)+O(nlogn)=O(nlogn)。

这东西跑得比EXP要快的多了,但是思考过程肯定没有EGF+EXP那样简洁。

于是乎现在可以考虑卡O(nlog^2n)的错解了(雾)

Code:

#include<algorithm>
#include<cstdio>
#define mod 167772161
#define G 3
#define Maxn 270000
using namespace std;
int n,k,r[Maxn<<2];
long long invn,invG;
long long fac[Maxn],inv[Maxn];
long long powM(long long a,long long t=mod-2)
{
  long long ans=1,buf=a;
  while(t){
    if(t&1)ans=(ans*buf)%mod;
    buf=(buf*buf)%mod;
    t>>=1;
  }return ans;
}
void NTT(long long *f,bool op,int n)
{
  for (int i=0;i<n;i++)
    if (r[i]<i)swap(f[r[i]],f[i]);
  for (int len=1;len<n;len<<=1){
    int w=powM(op==1 ? G:invG,(mod-1)/len/2);
    for (int p=0;p<n;p+=len+len){
      long long buf=1;
      for (int i=p;i<p+len;i++){
        int sav=f[i+len]*buf%mod;
        f[i+len]=f[i]-sav;
        if (f[i+len]<0)f[i+len]+=mod;
        f[i]=f[i]+sav;
        if (f[i]>=mod)f[i]-=mod;
        buf=buf*w%mod;
        }//F(x)=FL(x^2)+x*FR(x^2)
         //F(W^k)=FL(w^k)+W^k*FR(w^k)
         //F(W^{k+n/2})=FL(w^k)-W^k*FR(w^k)
    }
  }
}
long long g[Maxn<<2];
void rev(long long *f,int len)
{
  for (int i=0;i<len;i++)g[i]=f[i];
  for (int i=0;i<len;i++)f[len-i-1]=g[i];
}
//令f=f*g (mod x^lim) 
void times(long long *f,long long *gg,int len,int lim)
{
  int m=len+len,n;
  for(int i=0;i<len;i++)g[i]=gg[i];
  for(n=1;n<m;n<<=1);invn=powM(n);
  for(int i=len;i<n;i++)g[i]=0;
  for(int i=0;i<n;i++)
    r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
  NTT(f,1,n);NTT(g,1,n);
  for(int i=0;i<n;++i)f[i]=(f[i]*g[i])%mod;
  NTT(f,0,n);
  for(int i=0;i<lim;++i)f[i]=(f[i]*invn)%mod;
  for(int i=lim;i<n;++i)f[i]=0;
}
void Init(int lim)
{
  inv[1]=inv[0]=fac[0]=1;
  for (int i=1;i<=lim;i++)fac[i]=fac[i-1]*i%mod;
  for (int i=2;i<=lim;i++)
    inv[i]=inv[mod%i]*(mod-mod/i)%mod;
  for (int i=2;i<=lim;i++)inv[i]=inv[i-1]*inv[i]%mod;
  for (int i=1;i<=lim;i++)inv[i]=powM(fac[i]);
}
long long p[Maxn<<2];
//求出F(x-c) 
void fminus(long long *s,long long *f,int len,int c)
{
  c=mod-c;
  for (int i=0;i<len;i++)
    p[len-i-1]=f[i]*fac[i]%mod;
  long long buf=1;
  for (int i=0;i<len;i++,buf=buf*c%mod)
    s[i]=buf*inv[i]%mod;
  times(p,s,len,len);
  for (int i=0;i<len;i++)s[len-i-1]=p[i]*inv[len-i-1]%mod;
  for (int i=len;i<len+len;i++)s[i]=0;
}
long long f[Maxn<<2],s[Maxn<<2];
void solve(long long *f,int n)
{
  if (n==1){f[0]=0;f[1]=1;}
  else if (n&1){
    solve(f,n-1);f[n]=0;
    //再乘上(x-n+1)就好了
    for (int i=n;i>0;i--)
      f[i]=(f[i-1]+(mod-n+1)*f[i])%mod;
    f[0]=f[0]*(mod-n+1)%mod;
  }else {
    solve(f,n/2);
    //S(x)=F(x+n/2)
    fminus(s,f,n/2+1,n/2);
    times(f,s,n/2+1,n+1);
  }
}
void invp(long long *f,int len)
{
  for (int i=0;i<k+1;i++)s[i]=p[i]=0;
  //注意清空 
  long long *r=s,*rr=p;
  int n=1;for(;n<len;n<<=1);
  rr[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
    for (int i=0;i<len;i++)
      r[i]=rr[i]*2%mod;
    times(rr,rr,len/2,len);
    times(rr,f,len,len);
    for (int i=0;i<len;i++)
      rr[i]=(r[i]-rr[i]+mod)%mod;
  }for (int i=0;i<len;i++)
    f[i]=rr[i];
}
int main()
{
  scanf("%d%d",&n,&k);
  if (k>n){
    for (int i=0;i<=n;i++)printf("0 ");
    return 0;
  }invG=powM(G);
  Init(k);solve(f,k+1);
  for (int i=0;i<k+1;i++)f[i]=f[i+1];
  rev(f,k+1);
  for (int i=n-k+1;i<k+1;i++)f[i]=0;
  for (int i=k+1;i<n-k+1;i++)f[i]=0;
  invp(f,n-k+1);
  for (int i=0;i<k;i++)printf("0 ");
  for (int i=0;i<n-k+1;i++)printf("%lld ",f[i]);
  return 0;
}