半在线卷积小记
command_block
2020-05-05 22:29:14
**广告** : [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan)
- [P4721 【模板】分治 FFT](https://www.luogu.com.cn/problem/P4721)
**题意** : 给出多项式$F(x)$,求$G$使得:
- $G[0]=1$
- $G[n]=\sum\limits_{i=1}^nG[n-i]F[i]$
常识数据范围。
这个题可以直接求逆做，由于和主题无关我们就不介绍了。
此题中，某个位置的系数依赖于前面的计算，而又会向后面贡献。
这和一些数据结构优化 `DP` 题有点类似，实际上,这就是强制中序遍历转移的 `CDQ` 分治。
考虑每次计算跨越隔板的贡献。
式子是卷积的形式，且此时没有依赖性，可以直接用 `FFT` 计算，复杂度是 $O(n\log^2n)$ 的。
来模拟一下，令`F=<0,1,3,2>`
`G=< 1, 0, 0, 0 >`
`G=<[1 |0] 0, 0 >`
使用`<1>`卷上$F$的一部分 : `<0,1>` 得到`<*,1>`,将其加到对应的位置。
`*` 表示我们不在意那个位置的具体值。
`G=<[1 |1] 0, 0 >`
`G=<[1, 1 |0, 0]>`
使用 `<1,1>` 卷上 $F$ 的一部分 : `<0,1,3,2>` 得到 `<*,*,4,5>` ,将其加到对应的位置。
`G=<[1, 1 |4, 5]>`
`G=< 1, 1 [4 |5]>`
使用 `<4>` 卷上 $F$ 的一部分 : `<0,1>` 得到 `<*,4>` ,将其加到对应的位置。
注意，并非和后半部的 `<3,2>` 进行卷积。
`G=< 1, 1, 4, 9 >`
这就是答案了。
可以先填充成 $2$ 的幂,这样会比较好写。
由于我们不关心前半部分的结果，我们可以把循环卷积的次数定为 $O(len)$ 而非 $O(1.5len)$ ,这样就算溢出也只会损坏我们不关心的前半部分。这能够较为可观地优化常数，耗时为原来的`60~75%`。
又能发现，我们只会使用 $F$ 的某个前缀的 `DFT` 结果，可以预处理出来。耗时进一步减少到原来的`70~80%`。
```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#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};
for (int i=0;i<n;i++)f[i]=g[tr[i]];
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<<17))
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)return ;
int mid=(l+r)>>1,n=r-l;
cdq(l,mid);
cpy(s1,G+l,n/2);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;
cdq(mid,r);
}
int n,m,F[Maxn];
int main()
{
m=read();G[0]=1;
for (int i=1;i<m;i++)F[i]=read();
for (n=1;(n>>1)<m;n<<=1){
cpy(s1,F,n);NTT(s1,1,n);
cpy(s2+n,s1,n);
}cdq(0,n>>=1);
for (int i=0;i<m;i++)
printf("%d ",G[i]);
return 0;
}
```
- 多项式 $\exp$
在有标号计数的推导中我们经常要用到多项式 $\exp$ ，然而在考场上写 $O(n\log n)$ 的牛顿迭代 $\exp$ 性价比并不高。
这时就要请我们的半在线卷积 `Exp` 上场了,复杂度是 $O(n\log^2 n)$ 的，但由于常数小并不会比牛顿迭代慢。
( ~~我才不会告诉你们这玩意比我使劲卡常的牛顿迭代`Exp`快好多~~ )
设 $B(x)=e^{F(x)}$ ，两边求导可得 $B'(x)=B(x)F'(x)$
提取第 $n$ 次项系数可得 : $(n+1)B[n+1]=\sum\limits_{i=0}^{n}(i+1)F[i+1]B[n-i]$
$B[n]=\frac{1}{n}\sum\limits_{i=1}^{n}iF[i]B[n-i]$
这就回到了经典问题,常数$\frac{1}{n}$可以在分治边界乘上。
```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#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=1050000;
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<<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(int *g,bool op,int n)
{
tpre(n);
static ull f[Maxn<<1],w[Maxn<<1]={1};
for (int i=0;i<n;i++)f[i]=g[tr[i]];
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<<17))
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 F[Maxn],G[Maxn],s1[Maxn],s2[Maxn];
void cdq(int l,int r)
{
if (l+1==r){
if (l>0)G[l]=powM(l)*G[l]%mod;
return ;
}int mid=(l+r)>>1,n=r-l;
cdq(l,mid);
cpy(s1,G+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
cpy(s2,F,n);NTT(s2,1,n);
px(s1,s2,n);NTT(s1,0,n);
for (int i=n/2;i<n;i++)
G[l+i]=(G[l+i]+s1[i])%mod;
cdq(mid,r);
}
int n,m;
int main()
{
m=read();G[0]=1;
for (int i=0;i<m;i++)F[i]=1ll*read()*i%mod;
for (n=1;n<m;n<<=1);
cdq(0,n);
for (int i=0;i<m;i++)
printf("%d ",G[i]);
return 0;
}
```
- **更加真实的半在线卷积**
已知 $F[0]=G[0]=0,\ F[1]=1$
令 $G[n]=F[0...n]$ 的异或和。
$F[n]=\sum\limits_{i=1}^{n-1}G[i]F[n-i]$
前面的问题都是事先给定 $G$ ,此时 $G[n]$ 要在 $F[0...n]$ 都得到之后才能算出。
前面我们的策略是 : $F[l,mid]*G[0,r-l)\rightarrow F(mid,r)$
当 $l=0$ 时是这样的 : $F[0,mid]*G[0,r-l)\rightarrow F(mid,r)$
$G(mid,r-l)$ 还没有被正确计算出来。直接这样做是错误的。
正确的策略是 :
- $l=0$ 时
$F[0,mid]*G[0,mid]\rightarrow F(mid,r)$
这样,暂时只有 $F[mid+1]$ 是正确的。余者没有考虑 $F[0,mid]*G(mid,r-l)$ 的贡献。
- $l>0$ 时
$\left.\begin{matrix}G[l,mid]*F[0,r-l) \\ F[l,mid]*G[0,r-l) \end{matrix}\right\}\rightarrow F(mid,r)$
容易发现，此时左侧涉及到的所有系数都已经计算完毕。
下半部是经典的，上半部就是在补充 $l=0$ 时漏掉的贡献。
- 应用 : [[数学记录]P5900 无标号无根树计数](https://www.luogu.com.cn/blog/command-block/post-shuo-xue-ji-lu-p5900-wu-biao-hao-wu-gen-shu-ji-shuo)
- **更加真实的半在线卷积 II**
已知 $F[0]=1$ ,给出 $G$。
$F[n]=\sum\limits_{i=1}^{n}G[i]F^2[n-i]$
此时，当 $F$ 多算得 $F[n]$ 一位时，并不能马上快速更新 $F^2[n]$。
使用半在线卷积计算 $F^2$ ，然后并行计算即可。
- 应用 : [[数学记录]Uoj#50. 【UR #3】链式反应](https://www.luogu.com.cn/blog/command-block/post-50-ur-3-lian-shi-fan-ying)