目前已经是用 sin 和 cos 计算每个单位根的幂而不是乘法，但是压 $3$ 位或更多精度还是会炸。

void dif(cp a[]){
	constexpr int half=n/2,quarter=n/4;
	for(int i=0;i<quarter;i++){
		const cp w(cos(pi2/n*i),sin(pi2/n*i)),w3(cos(pi6/n*i),sin(pi6/n*i));
		const cp tmp1=a[i],tmp2=a[i+quarter],tmp3=a[i+half],tmp4=a[i+half+quarter];
		cp x=tmp1-tmp3,y=tmp2-tmp4;
		y=cp(y.imag(),-y.real());
		a[i]+=tmp3;
		a[i+quarter]+=tmp4;
		a[i+half]=(x-y)*w;
		a[i+half+quarter]=(x+y)*w3;
	}
	dif<half>(a);
	dif<quarter>(a+half);
	dif<quarter>(a+half+quarter);
}

（请不必理会奇怪的递归方式，总之跑的比迭代快）