[数学记录]Loj#6703. 小 Q 的序列
command_block
2020-05-03 22:35:46
戴项式入门题。做本题之前,请先跟我大喊 $\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)$$
- **希望分别求和**
三个项并不是很好,考虑把两个带 $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]$ ,记得要减去空子序列的贡献。
```cpp
#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`的常数是在复杂度瓶颈上的。规模小时暴力卷积可以有效地优化常数。
搞一搞就最优解了 : [评测记录](https://loj.ac/submission/807342)