分治FFT

分治FFT

概述

严格的分治FFT为“一般形式”，而我们把所有带有FFT的分治都称为分治FFT，所以这个名称并没有什么意义。
分治FFT代码比较复杂，低阶的比赛里应该不会出现。

一般形式

严格的分治FFT负责解决类似(displaystyle b_k=sum^{k-1}_{i=1}{b_ia_{k-i}})的数列，其中整个式子的形状符合卷积，且所求数列的每一项都依赖于前面的项。

已知数组(a)，求数组(b)，其中(displaystyle b_k=sum^{k-1}_{i=1}{b_ia_{k-i}})

考虑分治，对于区间([l,r))，先计算([l,m))的答案。对于(kgeq m)，将(b_k)拆成(displaystyle sum^{m-1}_{i=0}b_ia_{k-i})和(displaystyle sum^{k-1}_{i=m}b_ia_{k-i})两部分，后者可以给数组下标添加偏移量之后递归计算，而前者可以通过普通的卷积计算。

总时间复杂度为(Theta(n log^2 n))。

实现

假设我们求出了(l o mid)的答案，要求这些点对(mid+1 o r)的影响，那么对右半边点(x)的贡献为:(displaystyle w_x=sum_{i=l}^{mid} f[i] * g[x-i])这部分可以利用卷积来快速计算。计算完以后，答案直接加到答案数组就可以了。需要注意的是，如果要求左边点对右边点的影响，首先整个区间以左对该区间的贡献应该先求出。所以分治过程为先分治左边，在求出中间，然后在递归右边。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAXN=1e6+10;
const LL MOD=998244353,G=3,iG=332748118;
inline LL fpm(LL base,LL p){
	LL ret=1;
	while(p){
		if(p&1)
			ret=ret*base%MOD;
		base=base*base%MOD;
		p>>=1;
	}
	return ret;
}
void exgcd(LL a,LL b,LL &x,LL &y){
	if(!b){
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,y,x);
	y-=x*(a/b);
}
inline LL inv(LL x){
	LL ii,jj;
	exgcd(x,MOD,ii,jj);
	return (ii%MOD+MOD)%MOD;
}
int lim,rev[MAXN],N,L;
inline void prec(int l,int r){
	lim=1;
	L=0;
	while(lim<r-l){
		L++;
		lim<<=1;
	}
	for(int i=0;i<lim;i++)
		rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
}
inline void NTT(LL *a,int type){
	for(int i=0;i<lim;i++)
		if(i<rev[i])
			swap(a[i],a[rev[i]]);
	for(int hf=1;hf<lim;hf<<=1){
		int len=hf<<1;
		LL Wn=fpm(type==1?G:iG,(MOD-1)/len);
		for(int j=0;j<lim;j+=len){
			LL w=1;
			for(int k=0;k<hf;k++){
				LL t1=a[j+k],t2=a[j+k+hf]*w%MOD;
				a[j+k]=(t1+t2)%MOD;
				a[j+k+hf]=(t1-t2)%MOD;
				w=w*Wn%MOD;
			}
		}
	}
}
LL f[MAXN],g[MAXN],A[MAXN],B[MAXN];
inline void cdq(int l,int r){
	if(l+1>=r||l>=N)
		return;
	int mid=(l+r)>>1;
	cdq(l,mid);
	prec(l,r);
	memcpy(A,f+l,((r-l)/2)<<3);
	memcpy(B,g,(r-l)<<3);
	memset(A+(r-l)/2,0,((r-l)/2)<<3);
	NTT(A,1);
	NTT(B,1);
	for(int i=0;i<r-l;i++)
		A[i]=A[i]*B[i]%MOD;
	NTT(A,-1);
	int ii=inv(r-l);
	for(int i=0;i<r-l;i++)
		A[i]=A[i]*ii%MOD;
	for(int i=(r-l)>>1;i<r-l;i++)
		f[i+l]=(f[i+l]+A[i])%MOD;
	cdq(mid,r);
}
int main(){
	scanf("%d",&N);
	for(int i=1;i<N;i++)
		scanf("%lld",g+i);
	prec(0,N);
	f[0]=1;
	cdq(0,lim);
	for(int i=0;i<N;i++)
		printf("%lld ",(f[i]+MOD)%MOD);
	return 0;
}

例题

P4721 【模板】分治 FFT

见“实现”

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAXN=1e6+10;
const LL MOD=998244353,G=3,iG=332748118;
inline LL fpm(LL base,LL p){
	LL ret=1;
	while(p){
		if(p&1)
			ret=ret*base%MOD;
		base=base*base%MOD;
		p>>=1;
	}
	return ret;
}
void exgcd(LL a,LL b,LL &x,LL &y){
	if(!b){
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,y,x);
	y-=x*(a/b);
}
inline LL inv(LL x){
	LL ii,jj;
	exgcd(x,MOD,ii,jj);
	return (ii%MOD+MOD)%MOD;
}
int lim,rev[MAXN],N,L;
inline void prec(int l,int r){
	lim=1;
	L=0;
	while(lim<r-l){
		L++;
		lim<<=1;
	}
	for(int i=0;i<lim;i++)
		rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
}
inline void NTT(LL *a,int type){
	for(int i=0;i<lim;i++)
		if(i<rev[i])
			swap(a[i],a[rev[i]]);
	for(int hf=1;hf<lim;hf<<=1){
		int len=hf<<1;
		LL Wn=fpm(type==1?G:iG,(MOD-1)/len);
		for(int j=0;j<lim;j+=len){
			LL w=1;
			for(int k=0;k<hf;k++){
				LL t1=a[j+k],t2=a[j+k+hf]*w%MOD;
				a[j+k]=(t1+t2)%MOD;
				a[j+k+hf]=(t1-t2)%MOD;
				w=w*Wn%MOD;
			}
		}
	}
}
LL f[MAXN],g[MAXN],A[MAXN],B[MAXN];
inline void cdq(int l,int r){
	if(l+1>=r||l>=N)
		return;
	int mid=(l+r)>>1;
	cdq(l,mid);
	prec(l,r);
	memcpy(A,f+l,((r-l)/2)<<3);
	memcpy(B,g,(r-l)<<3);
	memset(A+(r-l)/2,0,((r-l)/2)<<3);
	NTT(A,1);
	NTT(B,1);
	for(int i=0;i<r-l;i++)
		A[i]=A[i]*B[i]%MOD;
	NTT(A,-1);
	int ii=inv(r-l);
	for(int i=0;i<r-l;i++)
		A[i]=A[i]*ii%MOD;
	for(int i=(r-l)>>1;i<r-l;i++)
		f[i+l]=(f[i+l]+A[i])%MOD;
	cdq(mid,r);
}
int main(){
	scanf("%d",&N);
	for(int i=1;i<N;i++)
		scanf("%lld",g+i);
	prec(0,N);
	f[0]=1;
	cdq(0,lim);
	for(int i=0;i<N;i++)
		printf("%lld ",(f[i]+MOD)%MOD);
	return 0;
}