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

· · 个人记录

戴项式入门题。做本题之前,请先跟我大喊 \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

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

直接考虑一个naive的DP,设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)

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

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!}(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的常数是在复杂度瓶颈上的。规模小时暴力卷积可以有效地优化常数。

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