参考这样一段普通的 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 函数可以直接复制粘贴得到,编写用时可能还不到前面的 51,而 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 加位逆序置换长,但实际上两部分代码是高度相似的,实际编写难度并不比后者高。
那么这样的代码长度(码量?)究竟应该怎么算呢?