关于代码实现的一些问题
  • 板块灌水区
  • 楼主zhangbo1000
  • 当前回复1
  • 已保存回复1
  • 发布时间2024/4/5 10:11
  • 上次更新2024/4/5 11:17:31
查看原帖
关于代码实现的一些问题
760291
zhangbo1000楼主2024/4/5 10:11

参考这样一段普通的 FFT:

void rev(const int& n,cp a[]){
	for(int i=0,j=0;i<n;i++){
		if(i>j)swap(a[i],a[j]);
		for(int l=(n>>1);(j^=l)<l;l>>=1);
	}
}
template<const int n>
void fft(cp a[]){
	constexpr int half=n>>1,quarter=n>>2;
	cp w(1,0);
	constexpr cp 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;
	}
    rev(n,a);
}

单论字数,这样的代码可以算是相当长,但实际上最长的 runfft 函数可以直接复制粘贴得到,编写用时可能还不到前面的 15\dfrac{1}{5},而 fft 函数虽然短,但是推公式需要一定时间,还需要注意细节,还需要特化边界。

再比如,如果我们要实现 DIF 和 DIT 混用以省去位逆序置换,于是写出了如下代码:

using cp=complex<double>;
const double pi2=2*acos(-1);
void dit(cp a[],const int& n){
	if(n<=1)return;
	const int half=n/2;
	dit(a,half);
	dit(a+half,half);
	cp w(1,0),wn(cos(pi2/n),-sin(pi2/n));
	for(int i=0;i<half;i++){
		cp x=a[i],y=w*a[i+half];
		a[i]=x+y;
		a[i+half]=x-y;
		w*=wn;
	}
}
void dif(cp a[],const int& n){
	if(n<=1)return;
	const int half=n/2;
	cp w(1,0),wn(cos(pi2/n),sin(pi2/n));
	for(int i=0;i<half;i++){
		cp x=a[i],y=a[i+half];
		a[i]=x+y;
		a[i+half]=w*(x-y);
		w*=wn;
	}
	dif(a,half);
	dif(a+half,half);
}

虽然代码比单用 DIT/DIF 加位逆序置换长,但实际上两部分代码是高度相似的,实际编写难度并不比后者高。

那么这样的代码长度(码量?)究竟应该怎么算呢?

2024/4/5 10:11
加载中...