在给 FFT 卡常，已经加上了模板递归，分裂基，DIF 和 DIT 混用，快读，三次变两次等优化（手写 complex 实测没用），请问还有什么别的能有明显优化的方法？（不考虑手动内联展开短多项式 FFT 和指令集，不需要优化后的代码，目前运行速度和 P3803 最优解差了大约 $5$ 倍。）

代码：

#include<iostream>
#include<cstdio>
#include<complex>
#include<cmath>
using namespace std;
using ll=long long;
using cp=complex<double>;
const int N=1<<21;
const double pi=acos(-1);
const double pi2=2*pi;
cp a[N];
template<const int n>
void fft(cp a[]){
	const int half=n>>1,quarter=n>>2;
	cp w(1,0),wn(cos(pi2/n),sin(pi2/n));
	for(int i=0;i<quarter;i++){
		cp w3=w*w*w;
		cp x=a[i]-a[i+half],y=a[i+quarter]-a[i+half+quarter];
		y=cp(y.imag(),-y.real());
		a[i]+=a[i+half];
		a[i+quarter]+=a[i+half+quarter];
		a[i+half]=(x-y)*w;
		a[i+half+quarter]=(x+y)*w3;
		w*=wn;
	}
	fft<half>(a);
	fft<quarter>(a+half);
	fft<quarter>(a+half+quarter);
}
template<>
void fft<2>(cp a[]){
	cp x=a[0],y=a[1];
	a[0]=x+y;
	a[1]=x-y;
}
template<>
void fft<1>(cp a[]){}
template<>
void fft<0>(cp a[]){}
void runfft(cp a[],const int& n){
	switch(n){
		case 1<<1:fft<1<<1>(a);break;
		case 1<<2:fft<1<<2>(a);break;
		case 1<<3:fft<1<<3>(a);break;
		case 1<<4:fft<1<<4>(a);break;
		case 1<<5:fft<1<<5>(a);break;
		case 1<<6:fft<1<<6>(a);break;
		case 1<<7:fft<1<<7>(a);break;
		case 1<<8:fft<1<<8>(a);break;
		case 1<<9:fft<1<<9>(a);break;
		case 1<<10:fft<1<<10>(a);break;
		case 1<<11:fft<1<<11>(a);break;
		case 1<<12:fft<1<<12>(a);break;
		case 1<<13:fft<1<<13>(a);break;
		case 1<<14:fft<1<<14>(a);break;
		case 1<<15:fft<1<<15>(a);break;
		case 1<<16:fft<1<<16>(a);break;
		case 1<<17:fft<1<<17>(a);break;
		case 1<<18:fft<1<<18>(a);break;
		case 1<<19:fft<1<<19>(a);break;
		case 1<<20:fft<1<<20>(a);break;
		case 1<<21:fft<1<<21>(a);break;
	}
}
template<const int n>
void ifft(cp a[]){
	const int half=n>>1,quarter=n>>2;
	ifft<half>(a);
	ifft<quarter>(a+half);
	ifft<quarter>(a+half+quarter);
	cp w(1,0),wn(cos(pi2/n),-sin(pi2/n));
	for(int i=0;i<quarter;i++){
		cp w3=w*w*w;
		cp tmp1=w*a[i+half],tmp2=w3*a[i+half+quarter];
		cp x=a[i],y=tmp1+tmp2;
		cp x1=a[i+quarter],y1=tmp1-tmp2;
		y1=cp(y1.imag(),-y1.real());
		a[i]+=y;
		a[i+quarter]+=y1;
		a[i+half]=x-y;
		a[i+half+quarter]=x1-y1;
		w*=wn;
	}
}
template<>
void ifft<2>(cp a[]){
	cp x=a[0],y=a[1];
	a[0]=x+y;
	a[1]=x-y;
}
template<>
void ifft<1>(cp a[]){}
template<>
void ifft<0>(cp a[]){}
void runifft(cp a[],const int& n){
	switch(n){
		case 1<<1:ifft<1<<1>(a);break;
		case 1<<2:ifft<1<<2>(a);break;
		case 1<<3:ifft<1<<3>(a);break;
		case 1<<4:ifft<1<<4>(a);break;
		case 1<<5:ifft<1<<5>(a);break;
		case 1<<6:ifft<1<<6>(a);break;
		case 1<<7:ifft<1<<7>(a);break;
		case 1<<8:ifft<1<<8>(a);break;
		case 1<<9:ifft<1<<9>(a);break;
		case 1<<10:ifft<1<<10>(a);break;
		case 1<<11:ifft<1<<11>(a);break;
		case 1<<12:ifft<1<<12>(a);break;
		case 1<<13:ifft<1<<13>(a);break;
		case 1<<14:ifft<1<<14>(a);break;
		case 1<<15:ifft<1<<15>(a);break;
		case 1<<16:ifft<1<<16>(a);break;
		case 1<<17:ifft<1<<17>(a);break;
		case 1<<18:ifft<1<<18>(a);break;
		case 1<<19:ifft<1<<19>(a);break;
		case 1<<20:ifft<1<<20>(a);break;
		case 1<<21:ifft<1<<21>(a);break;
	}
}
inline char gtc(){
	static char buf[100000],*p=0,*top=0;
	if(p==top)top=(p=buf)+fread(buf,1,100000,stdin);
	return *p++;
}
inline int read(){
	int x;
	char c;
	while((c=gtc())<'0');
	x=c^'0';
	while((c=gtc())>='0')x=(x<<1)+(x<<3)+(c^'0');
	return x;
}
int main(){
	int n,m;
	n=read();
	m=read();
	for(int i=0;i<=n;i++)a[i].real(read());
	for(int i=0;i<=m;i++)a[i].imag(read());
	int least=max(n+1,m+1)<<1,s=1;
	while(s<least)s<<=1;
	runfft(a,s);
	for(int i=0;i<s;i++)a[i]*=a[i];
	runifft(a,s);
	least=n+m;
	const double inv=0.5/s;
	for(int i=0;i<=least;i++)printf("%ld ",lround(a[i].imag()*inv));
	return 0;
}