半在线卷积
分治fft,是一种基于分治的、可以在\(O(n\log^2n)\)的复杂度内求解一类形如卷积形式的递推式的每一项的算法。具体的来举个例子。
【模板】分治fft
题意:给定数组\(g\)的前\(n\)项(从\(1\)开始),求\(f\)的前\(n\)项。其中\(f\)满足:
答案对\(998244353\)取模。
(当然你可以多项式求逆)
首先看到这个长得像卷积一样的递推式很自然地想到多项式科技。联想我们曾经cdq分治的套路,cdq优化dp的时候也是要先递归左边,然后处理左边对右边的影响,最后递归右边。这个也一样。我们分治处理左边和右边,合并的时候做一次fft计算左边对右边的贡献。
另外由于这题有模数所以写了ntt。fft乘爆了不想调所以不写了。
分治fft分治ntt
void calc(int a[],int b[]){
ntt(a,wl,1);ntt(b,wl,1);
for(int i=0;i<wl;i++)a[i]=1ll*a[i]*b[i]%mod;
ntt(a,wl,-1);
}
void solve(int l,int r){
if(l>=r)return;
int mid=(l+r)>>1;
solve(l,mid);
int len=r-l+1;wl=1;
get(len-2);//处理长度和反转结果
for(int i=0;i<wl;i++)a[i]=b[i]=0;
for(int i=l;i<=mid;i++)a[i-l]=f[i];//f往前平移
for(int i=1;i<=r-l;i++)b[i-1]=g[i];//把g的前若干项拉过来 因为f在l时g下标最大 为r-l 所以循环到r-l就行
calc(a,b);
for(int i=mid+1;i<=r;i++)f[i]=(f[i]+a[i-l-1])%mod;
//我们一开始把f和g的下标向左平移了l+1 所以统计答案也要平移l+1
solve(mid+1,r);
}
int main(){
scanf("%d",&n);f[0]=1;n--;
for(int i=1;i<=n;i++)scanf("%d",&g[i]);
solve(0,n);
for(int i=0;i<=n;i++)printf("%d ",f[i]);
return 0;
}
刚才那算暴力实现,现在提供一个比较好的实现。而且易于封装。
参考了 qwaszx 的实现。
首先我们半在线卷积一般是分二叉,复杂度 \(O(n\log^2n)\)。接下来我们试着加几个叉。设分为 \(B\) 叉,那么每次需要计算两两间的贡献。并不需要 \(B^2\) 次 \(O(\dfrac nB\log\dfrac nB)\) 的卷积,我们只需要记录所有点值就可以做到 \(B\) 次。那么复杂度为
取 \(B=O(\log n)\),得到复杂度 \(T(n)=O(\dfrac{n\log^2n}{\log\log n})\)。
然而跑得快需要一些细节。
- 叉数:\(B\) 通常取 \(8\) 或 \(16\)。
- 小范围暴力。
- 子问题大小取为 \(2\) 的幂,不浪费 DFT 长度。
- 可以在外面预处理每一层 \(g\) 的点值。具体的,假设子问题大小为 \(d\),第 \(i\) 个块到第 \(j\) 个块 \(g\) 的贡献区间是 \([(j-i-1)d,(j-i+1)d]\),可以预处理。
- 每次多项式乘法是 \(2d-1\) 次多项式和 \(d-1\) 次多项式乘,并只取 \([d,2d-1]\) 项系数,可以循环卷积。
这样整个半在线卷积就只需要留在外面一个暴力的接口。而且暴力也基本是一样的(因为都是长的像 \(f_n=c_n\sum_{i=0}^{n-1}f_ig_{n-i}\) 这种的),所以你甚至只需要留一个 Relax 函数做接口来处理 \(c_n\),然后在外边的函数部分处理别的东西。
一份板子。
void cdq(int l,int r,int dep,int f[],int g[],void brute(int f[],int g[],int l,int r)){
if(r-l+1<=32){
brute(f,g,l,r);return;
}
if(l>=r)return;
static int tmp[300010];
int d=1<<((dep-1)*3);
for(int i=0;;i++){
int L=l+i*d,R=min(r,L+d-1);
if(i){
for(int j=0;j<(d<<1);j++)tmp[j]=0;
for(int j=0;j<i;j++){
for(int k=0;k<(d<<1);k++){
int x=1ll*f1[dep][j][k]*g1[dep][i-j][k]%mod;
tmp[k]=add(tmp[k],x);
}
}
DIT(tmp,d<<1);
for(int i=L;i<=R;i++)f[i]=add(f[i],tmp[i-L+d]);
}
cdq(L,R,dep-1,f,g,brute);
if(R==r)return;
for(int j=0;j<(d<<1);j++)f1[dep][i][j]=0;
for(int j=L;j<=R;j++)f1[dep][i][j-L]=f[j];
DIF(f1[dep][i],d<<1);
}
}
void solve(int f[],int g[],int n,void brute(int f[],int g[],int l,int r)){
if(n<=128){
brute(f,g,0,n-1);return;
}
int len=1,dep=0;
while(len<n)len<<=3,dep++;len>>=3;
int *now=buf;
for(int i=1;i<=dep;i++){
int d=1<<((i-1)*3),mn=min((n-1)/d,7);
for(int j=1;j<=mn;j++){
int l=(j-1)*d+1,r=min(n-1,(j+1)*d-1);
f1[i][j-1]=now;now+=d<<1;
g1[i][j]=now;now+=d<<1;
for(int k=0;k<(1<<d);k++)g1[i][j][k]=0;
for(int k=l;k<=r;k++)g1[i][j][k-l+1]=g[k];
DIF(g1[i][j],d<<1);
}
}
cdq(0,n-1,dep,f,g,brute);
}
好了现在问题就是暴力怎么写。一些 \(O(n^2)\) 做常见幂级数运算的方法:
- 求逆:
- \(\ln\):
- \(\exp\):