题解 CF1349F2 【Slime and Sequences (Hard Version)】

command_block · 2020-05-15 19:45:17 · 题解

官方题解比较简略,这里我来补充一个详细一点的。

此题 : 国际计数水平.jpg

题意 : 定义好序列为 : 都是正整数 , 如果某个数 k 出现过 , 那么一定有 k-1 出现在其前面。

现在对 k=1...n ,分别统计其在长度为 n 的所有好序列中,出现次数的和。

Easy \ n\leq 5000 , 时限\texttt{2s}
Hard \ n\leq 10^5 , 时限\texttt{3s}

中国特色

计数工业总产值跃居世界第一!

这位老哥您先别走,能不能告诉我 GF 是什么意思, 拉格朗日反演又是啥?

0.前置芝士

分式域
(扩展)拉格朗日反演
欧拉数
二元 GF

1.转化

这个好序列的限制十分刻意,首先考虑转化一下。

事实上,可以建立起长度为 n 的排列和长度为 n 的好序列之间的对应关系。

在排列 p 相邻数之间填写不等号,序列中 a_{p_i} 的值就是 p_{1...i} 之间的小于号个数+1。

证明其满足好序列的条件 :

考虑一个 < ,以此为分界线,前面为 p_{i-1} ,后面为 p_{i} 。有 a_{p_{i-1}}+1=a_{p_i} ,因为中间隔了一个小于号。

同时又有 p_{i-1}<p_i(位置) ,也是因为小于号。实际上小于号就是把两个限制结合在了一起。

通过这个构造,也容易从好序列还原排列。具体证明略。

~~学长说这可以入选“十大炫酷映射”~~

以本人的憨憨水平,考场上根本想不出来这种精妙的构造。甚至还根本没去想等价转化,直接对着原问题分析,连正解的边都没摸着,直接被拒之门外。这可能也就是我正式比赛中没切过任何一道计数工业的原因吧……

目前大概想到了几点。

首先如果原问题复杂(条件信息量太大),感觉不可做,那么大概率是出题人做了神仙转化,要去寻找构造或者充要条件。

其次,可以统计需要计数的东西的总数,比如这个题,如果发现总数是 n! ,就更容易想到构造成排列。

2.初步分析

现在就是要对排列统计 < 的个数。不难发现其实就是 P5825 排列计数中的欧拉数。

设 \left\langle\begin{matrix}n \\ k\end{matrix}\right\rangle 为 : 长度为 n 的排列,小于号的个数为 k 的方案数。

我们分别考虑每个位置的贡献。

{\rm Ans}[k]=\sum\limits_{i=\max(k,1)}^n\left\langle\begin{matrix}i \\ k\end{matrix}\right\rangle\dbinom{n}{i}(n-i)!

意思是 : 枚举在好序列中产生 k 的位置。此时整个排列的方案数,就是对局部排列进行标号分配。后面的部分没有顺序限制,还能随意排列。

其实只有 {\rm Ans}[0] 会被 \max(k,1) 影响,我们可以特判,下界就会变得更漂亮。

存在一个O(n^2)递推欧拉数的做法,然后暴力计算上式,即可通过Easy版。

#include<algorithm>
#include<cstdio>
#define MaxN 5050
#define ll long long
using namespace std;
const int mod=998244353;
ll powM(ll a,int t=mod-2){
  ll ret=1;
  while(t){
    if (t&1)ret=ret*a%mod;
    a=a*a%mod;t>>=1;
  }return ret;
}
ll fac[MaxN],ifac[MaxN],s[MaxN],s2[MaxN],ans[MaxN];
int n; 
int main()
{
  scanf("%d",&n);
  fac[0]=1;
  for (int i=1;i<=n;i++)
    fac[i]=fac[i-1]*i%mod;
  ifac[n]=powM(fac[n]);
  for (int i=n;i;i--)
    ifac[i-1]=ifac[i]*i%mod;
  s2[0]=1;
  for (int i=1;i<=n;i++){
    s[0]=1;ans[0]=(ans[0]+ifac[i])%mod;
    for (int j=1;j<=i;j++){
      s[j]=(s2[j]*(j+1)+s2[j-1]*(i-j))%mod;
      ans[j]=(ans[j]+s[j]*ifac[i])%mod;
    }for (int j=0;j<=i;j++)s2[j]=s[j];
  }for (int i=0;i<n;i++)
    printf("%lld ",ans[i]*fac[n]%mod);
  return 0;
}

3.加速计算

这个式子涉及到了O(n^2)个的大量欧拉数,显然分离考虑是不行的。我们需要把 “黑盒” 欧拉数换成更易于批量分析的形式。

先设 \left|\begin{matrix}n \\ k\end{matrix}\right| 为 : 长度为 n 的排列,钦定 k 个小于号,其余随意的方案数(可能重复计数)。

那么,根据二项式反演,非常经典地有

\left|\begin{matrix}n \\ k\end{matrix}\right|=\sum\limits_{i=k}^n\dbinom{i}{k}\left\langle\begin{matrix}n \\ i\end{matrix}\right\rangle\quad\Longrightarrow\quad \left\langle\begin{matrix}n \\ k\end{matrix}\right\rangle=\sum\limits_{i=k}^n\dbinom{i}{k}(-1)^{i-k}\left|\begin{matrix}n \\ i\end{matrix}\right|

我们考虑 \left|\begin{matrix}n \\ k\end{matrix}\right| 如何计算,注意到一串连续的 < 可以视作一块,然后块中放置数的顺序就是固定的。

考虑单个块的EGF,为 \{0,1,1,1,\dots\}\longrightarrow e^x-1 (不考虑空集)

如果恰好需要 k 个 < ,则有 n-k 个块。

这是个有序组合并分配标号,能得到 \left|\begin{matrix}n \\ k\end{matrix}\right|=n![x^n](e^x-1)^{n-k}

则 \left\langle\begin{matrix}n \\ k\end{matrix}\right\rangle=n!\sum\limits_{i=k}^n\dbinom{i}{k}(-1)^{i-k}[x^n](e^x-1)^{n-i}

代入Easy的式子中,得到

{\rm Ans}[k]=\sum\limits_{i=k}^n\dbinom{n}{i}(n-i)!i!\sum\limits_{j=k}^i\dbinom{j}{k}(-1)^{j-k}[x^i](e^x-1)^{i-j}

=n!\sum\limits_{i=k}^n\sum\limits_{j=k}^i\dbinom{j}{k}(-1)^{j-k}[x^i](e^x-1)^{i-j}

=\dfrac{n!}{k!}\sum\limits_{j=k}^n\dfrac{(-1)^{j-k}j!}{(j-k)!}\sum\limits_{i=j}^n[x^i](e^x-1)^{i-j}

设S[k]=\sum\limits_{i=k}^n[x^i](e^x-1)^{i-k}

=\dfrac{n!}{k!}\sum\limits_{j=k}^n\dfrac{(-1)^{j-k}}{(j-k)!}S[j]j!

得知 S[0...k] 之后显然可以差卷积计算。

观察 S[k]=\sum\limits_{i=k}^n[x^i](e^x-1)^{i-k}
=\sum\limits_{i=k}^n[x^k]\Big(\dfrac{e^x-1}{x}\Big)^{i-k} =\sum\limits_{i=0}^{n-k}[x^k]\Big(\dfrac{e^x-1}{x}\Big)^{i}
形式统一了许多。

设 F(x)=\dfrac{e^x-1}{x} 。

我们的目标就是对 k=0...(n-1) 求出 S[k]=\sum\limits_{i=0}^{n-k}[x^k]F(x)^{i}

=[x^k]\bigg(\dfrac{1-F(x)^{n-k+1}}{1-F(x)}\bigg)

=[x^k]\dfrac{1}{1-F(x)}-[x^k]\dfrac{F(x)^{n-k+1}}{1-F(x)}

第一部分可以多项式求逆,但是注意到分母的常数项为 0 ,不能直接上。

考虑修改求逆的定义 : G(x)=x^kF(x), 使得 F(x) 有常数项，那么 \dfrac{1}{G(x)} = x^{-k}\dfrac{1}{F(x)}

附 : 这需要引入分式域(负数次数),拉格朗日反演在这里仍然成立,所以不用担心干扰推导。

通俗地讲,就是先把整个式子乘以 x^k ,相当于对分母除 x^k ,然后可以经典求逆。提取系数的时候要考虑多余的 x^k ,提取高k位,所以多项式的界要大一点。

这里实际上最低只有 [x^{-1}] 项,所以也不是很麻烦。

接下来考虑第二部分。(分母的逆同样拥有负指数,不过没关系)

把次数补齐 : [x^k]\dfrac{F(x)^{n-k+1}}{1-F(x)}=[x^{n+1}]\dfrac{\big(xF(x)\big)^{n-k+1}}{1-F(x)}

考虑加入一个元以区分信息。

=[x^{n+1}y^{n-k+1}]\sum\limits_{i=0}\dfrac{\big(xF(x)y\big)^i}{1-F(x)}

=[x^{n+1}y^{n-k+1}]\dfrac{1}{1-F(x)}\dfrac{1}{1-xF(x)y}

注意到右边是与 k 无关的,我们只需要求出 [x^{n+1}] 项关于 y 的多项式就得到了答案。

考虑使用拉格朗日反演,多项式复合系列理论。

但注意到 \dfrac{1}{1-F(x)}\dfrac{1}{1-xF(x)y} 除了 F(x) 还有额外的 x ,不容易用经典的多项式复合直接描述。

先设W(x)=xF(x) , H(x) 满足 \dfrac{W(x)}{H(W(x))}=x, 即 \dfrac{xF(x)}{H(W(x))}=x , 则 F(x)=H(W(x))

(我们使用 H(x) 给 xF(x) 消去了 x )

又设 G(x)=\dfrac{1}{1-H(x)}\dfrac{1}{1-xy}

注 : 手动爆算可求得 G'(x)=\dfrac{y+H'(x)-yH(x)-xyH'(x)}{\big(1-H(x)\big)^2\big(1-xy\big)^2}

那么有 G(W(x))=\dfrac{1}{1-H(W(x))}\dfrac{1}{1-W(x)y}=\dfrac{1}{1-F(x)}\dfrac{1}{1-xF(x)y}

(现在我们以经典多项式复合描述了问题)

现在的问题就是 [x^{n+1}]G(W(x))。

设 P(x) 为 W(x) 的复合逆 ,不难发现 P(x)=\dfrac{x}{H(x)} , 则 \dfrac{P(x)}{x}=\dfrac{1}{H(x)}。

使用扩展拉格朗日反演可得

[x^{n+1}]H(W(x))=\dfrac{1}{n+1}[x^n]\dfrac{y+H'(x)-yH(x)-xyH'(x)}{\big(1-H(x)\big)^2\big(1-xy\big)^2}\left(\dfrac{P(x)}{x}\right)^{-(n+1)}

忽略常数 \frac{1}{n+1} 并代换

=[x^n]\dfrac{y\big(1-H(x)\big)+\big(1-xy\big)H'(x)}{\big(1-H(x)\big)^2\big(1-xy\big)^2}H(x)^{n+1}

=[x^n]\left(\dfrac{y}{\big(1-H(x)\big)\big(1-xy\big)^2}+\dfrac{H'(x)}{\big(1-H(x)\big)^2\big(1-xy\big)}\right)H(x)^{n+1}

展开关于 y 的封闭形式

=[x^n]\left(\dfrac{1}{\big(1-H(x)\big)}\sum\limits_{i=0}(i+1)x^iy^{i+1}+\dfrac{H'(x)}{\big(1-H(x)\big)^2}\sum\limits_{i=0}(xy)^i\right)H(x)^{n+1}

考虑提取 [y^m] 系数

[x^ny^m]\left(\dfrac{1}{\big(1-H(x)\big)}\sum\limits_{i=0}(i+1)x^iy^{i+1}+\dfrac{H'(x)}{\big(1-H(x)\big)^2}\sum\limits_{i=0}(xy)^i\right)H(x)^{n+1}

=[x^n]\left(\dfrac{mx^{m-1}}{\big(1-H(x)\big)}+\dfrac{H'(x)x^m}{\big(1-H(x)\big)^2}\right)H(x)^{n+1}

=[x^{n-m+1}]m\dfrac{H(x)^{n+1}}{\big(1-H(x)\big)}+[x^{n-m}]\dfrac{H'(x)H(x)^{n+1}}{\big(1-H(x)\big)^2}

分别求出 \dfrac{H(x)^{n+1}}{\big(1-H(x)\big)} 和 \dfrac{H'(x)H(x)^{n+1}}{\big(1-H(x)\big)^2} 即可,问题变为了求 H(x)。

常规的多项式复合逆肯定是不行的,注意到问题的特殊性 : W(x)=xF(x)=e^x-1

容易构造 \ln(e^x-1+1)=x , 那么有 P(x)=\ln(x+1) ,则 H(x)=\dfrac{x}{\ln(x+1)}。

( 此处有 1-H(x) 的常数项为 0 ,仍需要分式域,记得偏移次数 )

综上,我们以 O(n\log n) 或者 O(n\log^2n) 的复杂度解决了本题。

实现中使用了大力倍增快速幂,时限比较松就卡过去了。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define clr(f,n) memset(f,0,sizeof(ll)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(ll)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=135000;
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(ll *g,bool op,int n)
{
  #define ull unsigned long long
  tpre(n);
  static ull f[Maxn<<1],w[Maxn]={1};
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }      
    if (l==(1<<10))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(ll *f,ll *g,int n)
{for(int i=0;i<n;++i)f[i]=f[i]*g[i]%mod;}
void times(ll *f,ll *g,int m)
{
  int n=1;while(n<m+m)n<<=1;
  NTT(f,1,n);NTT(g,1,n);
  px(f,g,n);NTT(f,0,n);clr(f+m,n-m);
}
void invp(ll *f,int m)
{
  static ll sav[Maxn<<1],w[Maxn<<1],r[Maxn<<1];
  int n;for (n=1;n<m;n<<=1);
  w[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
    for (int i=0;i<(len>>1);i++)
      r[i]=(w[i]<<1)%mod;
    cpy(sav,f,len);
    NTT(w,1,len<<1);px(w,w,len<<1);
    NTT(sav,1,len<<1);px(w,sav,len<<1);
    NTT(w,0,len<<1);clr(w+len,len);
    for (int i=0;i<len;i++)
      w[i]=(r[i]-w[i]+mod)%mod;
  }cpy(f,w,m);clr(sav,n+n);clr(w,n+n);clr(r,n+n);
}
void dao(ll *f,int m){
  for (int i=1;i<m;i++)
    f[i-1]=f[i]*i%mod;
  f[m-1]=0;
}
void powp(ll *f,int m,int t)
{
  static ll g[Maxn<<1];
  int n=1;while(n<m+m)n<<=1;g[0]=1;
  while(t){
    NTT(f,1,n);
    if (t&1){
      NTT(g,1,n);
      px(g,f,n);NTT(g,0,n);
      clr(g+m,n-m);
    }px(f,f,n);NTT(f,0,n);
    clr(f+m,n-m);t>>=1;
  }cpy(f,g,m);clr(f+m,n-m);clr(g,n);
}
ll fac[Maxn],ifac[Maxn];
void Init(int m)
{
  fac[0]=1;
  for (int i=1;i<=m;i++)
    fac[i]=fac[i-1]*i%mod;
  ifac[m]=powM(fac[m]);
  for (int i=m;i;i--)
    ifac[i-1]=ifac[i]*i%mod;
}
int n,m;
ll F[Maxn<<1],G[Maxn<<1],S[Maxn<<1],H[Maxn<<1],H1[Maxn<<1],H2[Maxn<<1];
int main()
{
  scanf("%d",&n);Init(n+6);
  int len=1;while(len<n+n+4)len<<=1;
  m=n+2;
  for (int i=0;i<m+1;i++)G[i]=mod-ifac[i+2];
  invp(G,m+1);cpy(S,G+1,m+1);
  for (int i=0;i<m+1;i++){
    H[i]=fac[i]*ifac[i+1]%mod;
    if (i&1)H[i]=mod-H[i];
  }invp(H,m+1);cpy(H1,H,m);
  powp(H1,m,n+1);
  cpy(G,H,m+1);dao(G,m+1);
  cpy(H2,H1,m);times(H2,G,m);clr(G,len);
  for (int i=0;i<m;i++)G[i]=mod-H[i+1];
  invp(G,m);
  NTT(G,1,len);NTT(H1,1,len);NTT(H2,1,len);
  px(H1,G,len);NTT(H1,0,len);clr(H1+m,len-m);
  px(H2,G,len);NTT(H2,0,len);clr(H2+m,len-m);NTT(H2,1,len);
  px(H2,G,len);NTT(H2,0,len);clr(H2+m,len-m);
  for (int i=1,sav=powM(n+1);i<n+2;i++){
    int m=n-i+1;
    S[i]=(S[i]-(H1[n-m+2]*m+H2[n-m+2])%mod*sav)%mod;
  }S[0]=n;
  for (int i=0;i<n;i++)S[i]=(S[i]+mod)%mod*fac[i]%mod;
  for (int i=0;i<n;i++)F[i]=(i&1) ? mod-ifac[i] : ifac[i];
  clr(F+n,len-n);clr(S+n,len-n);reverse(S,S+n);
  NTT(F,1,len);NTT(S,1,len);
  px(S,F,len);NTT(S,0,len);
  clr(S+n,len-n);reverse(S,S+n);
  for (int i=0;i<n;i++)
    printf("%lld ",S[i]*fac[n]%mod*ifac[i]%mod);
  return 0;
}