command_block 的博客

command_block 的博客

[数学记录]Loj#6703. 小 Q 的序列

posted on 2020-05-03 22:35:46 | under 记录 |

戴项式入门题。做本题之前,请先跟我大喊 $\texttt{stO djq Orz!}$

还有 $\texttt{stO EI Orz!}$ & $\texttt{stO dls Orz!}$

题意 : 定义一个序列的权值为$\prod\limits_{i=1}^k(a_i+i)$

给出一个长度为$n$的序列,求所有非空子序列(不要求连续)的权值和。

或者,对于长度为$1...n$的子序列分别求和。

答案对 $998244353$ 取模 , $n\leq 10^5$ , 时限$\texttt{3s}$。


参考文献 : https://codeforces.com/blog/entry/76447

题意非常自然,但是第一眼看去并没有什么直接的思路。

直接考虑一个naiveDP,设$f[i][j]$为前$i$个位置,已经选定$j$个位置的权值和。

有转移 :

$$f[i][j]=f[i-1][j]+f[i-1][j-1]*(a_i+j)$$

仍然不会,考虑套路地写出OGF,令$F_i(x)=\sum\limits_{j=0}f[i][j]x^j$

我们发现 $f[i][j-1]$ 对应的系数却是 $j$ ,这并不容易用经典的求导表示(你会发现次数错开了)。

我们考虑变换一下第二维,令$k=i-j$。

有转移 :

$$f[i][k]=f[i-1][k-1]+f[i-1][k]*(a_i+i-k)$$

设$b_i=a_i+i$ (后面别忘了这一步)

然后我们就有多项式形式的转移 :

$$F_i(x)=xF_{i-1}(x)+b_iF_{i-1}(x)-xF'_{i-1}(x)$$

  • 希望分别求和

    三个项并不是很好,考虑把两个带 $x$ 的项合并。

    求导又要和不求导的合并,还带个 $-$ 号,使用坚挺的 $e^{-x}$ 吧!

    $$F_i(x)e^{-x}=xF_{i-1}(x)e^{-x}+b_iF_{i-1}(x)e^{-x}-xF'_{i-1}(x)e^{-x}$$

    $$=x(F_{i-1}(x)e^{-x}-F'_{i-1}(x)e^{-x})+b_iF_{i-1}(x)e^{-x}$$

    回忆 : $\big(F(x)*G(x)\big)'=F'(x)*G(x)+F(x)*G'(x)$

    那么就有

    $$(F(x)e^{-x})'=F'(x)e^{-x}-F(x)e^{-x}$$

    前面一切精妙的构造就是为了凑出这个玩意。

    $$\text{原式}=x(F_{i-1}(x)e^{-x})'+b_iF_{i-1}(x)e^{-x}$$

    为了避免啰嗦,定义 $H_i(x)=F_i(x)e^{-x}$ 。

    则有更为简洁的转移:

    $$H_i(x)=xH'_{i-1}(x)+b_iH_{i-1}$$

    现在可以考虑系数,得到 :

    $$H[i][j]=j*H[i-1][j]+b_i*H[i-1][j]=(j+b_i)H[i-1][j]$$

    出人意料地简洁! 由于次数是对齐的,现在我们能观察到 :

    $$H[n][j]=H[0][j]*\prod_{i=1}^n(j+b_i)$$

    定义$P(x)=\prod_{i=1}^n(x+b_i)$,然后跑个多点求值就好。复杂度$O(n\log^2n)$。

    这样可以做到对长度为 $1...n$ 的子序列分别求和。

  • 只希望整体求和

如果我们只希望整体求和,可以直接大力搞那个三项式。

$$F_i(x)=(x-x{\small\frac{d}{dx}}+b_i)F_{i-1}(x)$$

仍然考虑计算$P(x)=\prod_{i=1}^n(x+b_i)$,然后将$(x-x\frac{d}{dx})$代入进去。

$$F_n(x)=\sum\limits_{i=0}^nP[i](x-x{\small\frac{d}{dx}})^i$$

让我们考虑$(x-x\frac{d}{dx})^k$对应的幂级数究竟是什么。我们只要求得其系数和就能解决问题。

设$Q_k(x)=(x-x\frac{d}{dx})^k$,考虑递推式可得

$$Q_k[n]=Q_{k-1}[n-1]-nQ_{k-1}[n]$$

这和第二类斯特林数的递推式十分相似,考虑第二类斯特林数递推的组合意义。

“ $n$ 个有区别的球放入 $m$ 个无区别的盒子里,无空盒的方案数 ”

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

意思是,新加入一个球可以自成一盒,也可以投入前面的任何一个盒子里。

现在就相当于,每往盒子里投多一个数,就要把贡献取反。

这就能得到$Q_k(x)=\sum\limits_{i=0}^k\begin{Bmatrix}k\\i\end{Bmatrix}(-1)^{k-i}x^i$,而系数和就是把 $x$ 去掉。

考虑计算 $\sum\limits_{i=0,k=0}\begin{Bmatrix}k\\i\end{Bmatrix}(-1)^{k-i}x^k$ ,那么 $[x^k]$ 项的系数就是 $(x-x\frac{d}{dx})^k$ 的系数和。

我们就是要求这个多项式,设为 $H(x)$。

直接算并不容易,考虑前置 $i$ 并且 EGF : $\sum\limits_{k=0}\begin{Bmatrix}k\\i\end{Bmatrix}(-1)^{k-i}\dfrac{x^k}{k!}$

$=(-1)^k\sum\limits_{k=0}\begin{Bmatrix}k\\i\end{Bmatrix}\dfrac{(-x)^k}{k!}$

  • 回忆 : $\frac{1}{k!}(e^x-1)^k=\sum\limits_{k=0}\begin{Bmatrix}k\\i\end{Bmatrix}\dfrac{x^k}{k!}$

原式 $=\frac{1}{k!}(1-e^{-x})^k$

然后枚举 $i$ 并求和 :

$H(x)=\sum\limits_{i=0}\dfrac{(1-e^{-x})}{i!}=\exp(1-e^{-x})$

于是乎跑个 $\exp$ 就好了,注意要还原。

最后的答案就是 $\sum\limits_{i=0}^nP[i]H[i]$ ,记得要减去空子序列的贡献。

#include<algorithm>
#include<cstring>
#include<cstdio>
#include<ctime>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=135000;
inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,int 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],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(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn],w[Maxn>>1]={1};
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<4)+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(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int G[Maxn],s1[Maxn],s2[Maxn<<1];
void cdq(int l,int r)
{
  if (l+1==r){
    if (l>0)G[l]=powM(l)*G[l]%mod;
    return ;
  }if (l+1==r)return ;
  int mid=(l+r)>>1,n=r-l;
  cdq(l,mid);
  cpy(s1,G+l,n/2);
  if (n>=64){
    clr(s1+(n/2),n/2);NTT(s1,1,n);
    px(s1,s2+n,n);NTT(s1,0,n);
    for (int i=n/2;i<n;i++)
    G[l+i]=(G[l+i]+s1[i])%mod;
  }else {
    for (int i=n/2;i<n;i++)
      for (int j=0;j<(n>>1);j++)
        G[l+i]=(G[l+i]+1ll*s1[j]*s2[i-j])%mod;
  }cdq(mid,r);
}
void exp(int *F,int m)
{
  G[0]=1;
  for (int i=0;i<m;i++)F[i]=1ll*i*F[i]%mod;
  int n=1;
  for (;(n>>1)<m;n<<=1){
    cpy(s1,F,n);NTT(s1,1,n);
    cpy(s2+n,s1,n);
  }for (int i=0;i<64;i++)s2[i]=F[i];
  cdq(0,n>>=1);cpy(F,G,m);
}
int fac[Maxn],ifac[Maxn];
void Init(int lim)
{
  fac[0]=1;
  for (int i=1;i<=lim;i++)
    fac[i]=1ll*fac[i-1]*i%mod;
  ifac[lim]=powM(fac[lim]);
  for (int i=lim;i;i--)
    ifac[i-1]=1ll*ifac[i]*i%mod;
}
int _s[Maxn<<1],*tp=_s;
struct Data{int *f,c;}t[Maxn];
int n,F[Maxn];
int main()
{
  Init(n=read());
  for (int i=1,x;i<=n;i++){
    x=read();
    t[i].f=tp;tp+=(t[i].c=2);
    t[i].f[0]=x+i;t[i].f[1]=1;
  }
  int tot=n;
  while(tot>1){
    for (int i=2+(tot&1);i<=tot;i+=2){
      int c1=t[i-1].c,c2=t[i].c;
      cpy(s1,t[i-1].f,c1);cpy(s2,t[i].f,c2);
      if (1ll*c1*c2<=26ll*max(c1,c2)){
        clr(F,c1+c2-1);
        for (int i=0;i<c1;i++)
          for (int j=0;j<c2;j++)
            F[i+j]=(F[i+j]+1ll*s1[i]*s2[j])%mod;
        cpy(t[i-1].f,F,c1+c2-1);
      }else {
        int n=1;for (;n<c1+c2;n<<=1);
        clr(s1+c1,n-c1);clr(s2+c2,n-c2);
        NTT(s1,1,n);NTT(s2,1,n);
        px(s1,s2,n);NTT(s1,0,n);
        cpy(t[i-1].f,s1,c1+c2-1);
      }t[i-1].c=c1+c2-1;
      t[(i+1)/2]=t[i-1];
    }tot=(tot+1)/2;
  }
  for (int i=0;i<=n;i++)
    F[i]=(i&1) ? ifac[i] : mod-ifac[i];
  F[0]++;exp(F,n+1);
  int ans=0;
  for (int i=0;i<=n;i++)
    ans=(ans+1ll*F[i]*fac[i]%mod*_s[i])%mod;
  printf("%d",(ans-1+mod)%mod);
  return 0;
}

要注意的是,分治FFT和半在线卷积时,小规模NTT的常数是在复杂度瓶颈上的。规模小时暴力卷积可以有效地优化常数。

搞一搞就最优解了 : 评测记录