一个计数题

我也不知道在哪里见的题 qwq

description

给定 \(n,k\)，定义一个满二叉树（每个非叶子节点都有两个儿子的二叉树）权值为其每条从根出发的链经过的向左的边的数量的最大值。对于每个 \(i\in[1,n]\) 求出 \(i\) 个恰有叶子的权值不超过 \(k\) 的满二叉树的数量。

• \(n,k\leq 5000\)

solution

多项式的题解

设 \(f_{i,j}\) 表示 \(i\) 个叶子结点，最长的左偏链长度是 \(j\)，设 \(s_{i,j}=\sum\limits_{k=0}^j f_{i,k}\)，有转移：

• \(f_{1,0}=1\)
• \(f_{i,j}=\sum\limits_{x=0}^{i} s_{i-x,j-1}f_{x,j-1}+s_{x,j-1}f_{i-x,j}\)

设 \(F_j(x)\) 是 \(\{f_{i,j}\}_{i=0}^{+\infty}\) 的生成函数，\(S_j(x)\) 是 \(\{s_{i,j}\}_{i=0}^{+\infty}\) 的生成函数。我们有：

• \(S_j(x)=S_{j-1}(x)+F_{j}(x)\)
• \(F_j(x)=\dfrac{F_{j-1}(x)S_{j-1}(x)}{1-S_{j-1}(x)}\)

这时候直接多项式乘法加多项式求逆就可以做了。答案就是 \(S_{k-1}(x)\) 的各项系数。时间复杂度 \(O(n^2\log n)\)。但是常数巨大而且复杂度比正解劣，只能跑 65 分 awa。

注意到 \(S_j(x)=S_{j-1}(x)+F_{j}(x)\)，第二个转移式可以写成 \(S_{j}(x)=\dfrac{1-S_{j-2}(x)}{1-S_{j-1}(x)} S_{j-1}(x)\)。（设 \(S_{-1}(x)=0\)）

把右边的 \(S_{j-1}(x)\) 按照转移式展开，可以得到，\(S_{k-1}(x)=\dfrac{x}{1-\dfrac{x}{1-\dfrac{x}{1-\dots}}}\)。连分数有 \(k-1\) 层，最下面是个 \(S_0(x)=x\)。

然后按照 noi2021 密码箱的经典做法，把分数的分子分母塞到向量里，然后这么做的一次转移就可以用矩阵刻画了。

具体来说，\(\dfrac{a}{b}\) 用向量 \([a,b]\) 表示，乘上矩阵 \(A=\begin{bmatrix}1&1 \\ -x & 0 \end{bmatrix}\) 就转移好了。

然后多项式矩阵快速幂。

\(O(n\log^2 n)\)

常数巨大，没有 \(O(n^2)\) 正解跑得快

code

#include<bits/stdc++.h>

using namespace std;

using E=long long;
using ui=uint32_t;
using lint=__int128;
constexpr E inf=1e16,mod=998244353;
const int N=10010;

namespace Comb{
  vector<E> fac(1,1ll),ifac(1,1ll);

  inline E ksm(E a,E b,const E MOD=mod){
    E ret=1;
    assert(b>=0);
    while(b){
      if(b&1) ret=ret*a%MOD;
      a=a*a%MOD;
      b>>=1;
    }
    return ret;
  }

  void CombPrework(int n){
    fac.resize(n+1),ifac.resize(n+1);
    fac[0]=ifac[0]=1;
    for(int i=1; i<=n; i++) fac[i]=fac[i-1]*i%mod;
    ifac[n]=ksm(fac[n],mod-2);
    for(int i=n-1; i; i--) ifac[i]=ifac[i+1]*(i+1)%mod;
  }

  inline E C(E n,E m){
    if(n<0||m<0||n<m) return 0;
    while(fac.size()<=n){
      fac.emplace_back(fac.back()*(fac.size())%mod);
      ifac.emplace_back(ifac.back()*ksm(ifac.size(),mod-2)%mod);
    }
    return fac[n]*ifac[m]%mod*ifac[n-m]%mod;
  }

  inline E A(E n,E m){
    if(m<0) return 0;
    return C(n,m)*fac[m]%mod;
  }
};
using namespace Comb;

namespace Poly{
  typedef vector<E> poly;
  using namespace Comb;
  const unsigned __int128 brt=((unsigned __int128)1<<64)/mod;
  inline int rdc(E a){ return a-mod*(brt*a>>64); }
  inline E add(E a,E b){ return a+b>=mod?a+b-mod:a+b; }
  inline E del(E a,E b){ return a-b<0?a-b+mod:a-b; }
  inline E mul(E a,E b){ return a*b%mod; }

  void cpy(poly &x,poly &y,int n){
    x.resize(n);
    assert(y.size()>=n);
    for(int i=0; i<n; i++){
      x[i]=y[i];
    }
  }

  poly Add(const poly &x,const E k){
    auto ret=x;
    if(ret.size()==0) ret.emplace_back(k%mod);
    else ret[0]=add(ret[0],k);
    return ret;
  }

  poly Add(const poly &x,const poly &y){
    poly ret=x;
    for(int i=0; i<x.size(); i++){
      if(i<y.size()) ret[i]=add(x[i],y[i]);
      else return ret;
    }
    for(int i=x.size(); i<y.size(); i++) ret.emplace_back(y[i]);
    return ret;
  }

  poly Minus(const poly &x){
    auto ret=x;
    for(int i=0; i<ret.size(); i++) ret[i]=mod-x[i];
    return ret;
  }

  poly px(const poly &x,const poly &y,int t){
    poly ret=x;
    assert(t<=x.size()&&t<=y.size());
    for(int i=0; i<t; i++) ret[i]=mul(ret[i],y[i]);
    return ret;
  }

  vector<ui> rev;
  vector<vector<E>> gg,invgg;
  int lstn;
  void prework_for_ntt(int n){
    if(n==lstn) return ;
    rev.resize(n+1);
    gg.clear(),invgg.clear();
    for(int i=1; i<n; i++) rev[i]=(rev[i^(i&-i)]|(1<<(__lg(n)-1-__lg(i&-i))));
    for(int i=2; i<=n; i<<=1){
      E w=ksm(3,(mod-1)/i),iw=ksm(w,mod-2);
      vector<E> now1,now2;
      now1.emplace_back(1),now2.emplace_back(1);
      int len=i>>1;
      for(int i=1; i<len; i++){
        now1.emplace_back(mul(now1.back(),w));
        now2.emplace_back(mul(now2.back(),iw));
      }
      gg.emplace_back(now1),invgg.emplace_back(now2);
    }
    return lstn=n,void();
  }

  void ntt(poly &p,int n,int op){
    if(op==1){
      while(p.size()<n){
        p.emplace_back(0);
      }
    }
    prework_for_ntt(n);
    for(int i=1; i<n; i++) if(i<rev[i]) swap(p[i],p[rev[i]]);

    for(int i=2,t=0; i<=n; t++,i<<=1){
      int len=i>>1;
      for(int pos=0; pos<n; pos+=i){
        for(int j=pos; j<pos+len; j++){
          E u=p[j],v=mul(p[j+len],op==-1?invgg[t][j-pos]:gg[t][j-pos]);
          p[j]=add(u,v),p[j+len]=del(u,v);
        }
      }
    }

    if(op==-1){
      E tmp=ksm(n,mod-2);
      for(int i=0; i<n; i++) p[i]=mul(p[i],tmp);
    }
  }

  poly mul(const poly &x,const poly &y,int n){
    poly tmp1=x,tmp2=y;
    int t;
    for(t=1;t<=(x.size()+y.size());t<<=1);
    ntt(tmp1,t,1),ntt(tmp2,t,1);
    tmp1=px(tmp1,tmp2,t);
    ntt(tmp1,t,-1);
    while(tmp1.size()>n) tmp1.pop_back();
    while(tmp1.size()&&tmp1.back()==0) tmp1.pop_back();
    while(tmp1.size()<n) tmp1.emplace_back(0);
    return tmp1;
  }

  poly square(const poly &x,int n){
    poly tmp1=x;
    int t;
    for(t=1;t<=(x.size()*2);t<<=1) ;
    ntt(tmp1,t,1);
    tmp1=px(tmp1,tmp1,1);
    ntt(tmp1,t,-1);
    while(tmp1.size()>n) tmp1.pop_back();
    while(tmp1.size()&&tmp1.back()==0) tmp1.pop_back();
    while(tmp1.size()<n) tmp1.emplace_back(0);
    return tmp1;
  }

  poly polyinv(const poly &A,int n){
    auto x=A;
    while(x.size()<n) x.emplace_back(0);
    poly tmp,now,ret;
    ret.resize(n+1);
    ret[0]=ksm(x[0],mod-2);
    for(int i=2; i<=n; i<<=1){
      cpy(now,x,i);
      tmp.resize(i);
      for(int j=0; j<(i>>1); j++) tmp[j]=mul(ret[j],2);
      ntt(ret,i<<1,1); ret=px(ret,ret,i<<1);
      ntt(now,i<<1,1); ret=px(ret,now,i<<1);
      ntt(ret,i<<1,-1);
      while(ret.size()>i) ret.pop_back();
      for(int j=0; j<i; j++){
        ret[j]=del(tmp[j],ret[j]);
      }
    }
    return ret;
  }

};
using namespace Poly;

E n,k;

/*

1000 1000

5000 5000

*/

struct Matrix{
  poly s[2][2];
  int t;

  Matrix (int sz){
    t=sz;
  }

  friend Matrix operator *(const Matrix &x,const Matrix &y){
    Matrix z(x.t);
    for(int i=0; i<2; i++){
      for(int j=0; j<2; j++){
        for(int k=0; k<2; k++){
          z.s[i][j]=Add(z.s[i][j],mul(x.s[i][k],y.s[k][j],z.t));
        }
      }
    }
    return z;
  }
};

int main(){

   freopen("tree.in","r",stdin);
   freopen("out.out","w",stdout);

  cin>>n>>k;
  auto stttttt=clock();
  int t;
  for(t=1;t<=n;t<<=1);
  Matrix trans(t);

  trans.s[0][0].emplace_back(1),trans.s[0][1].emplace_back(1);
  trans.s[1][0].emplace_back(0),trans.s[1][0].emplace_back(mod-1);

  Matrix ret(t);
  ret.s[0][0].emplace_back(1),ret.s[1][1].emplace_back(1);
  k--;
  while(k){
    if(k&1) ret=ret*trans;
    trans=trans*trans;
    k>>=1;
  }

  Matrix init(t);
  init.s[0][0].emplace_back(1),init.s[0][0].emplace_back(mod-1);
  init.s[0][1].emplace_back(1);

  init=init*ret;

  auto s=mul(init.s[0][0],polyinv(init.s[0][1],t),t);

  /*poly f,s;s.emplace_back(0); 这是暴力做法
  s.emplace_back(1);
  f.emplace_back(0),f.emplace_back(1);

  int t;
  for(t=1;t<=n;t<<=1);
  for(int j=1; j<k; j++){
    poly tmp2=polyinv(Add(Minus(s),1ll),t);
    s=mul(f,tmp2,t);
  }
  */
  for(int i=1; i<=n; i++) cout<<(mod-s[i]%mod)%mod<<'\n';

  cerr<<(clock()-stttttt)<<endl;

  return 0;
}