我目前正在研究需要FFT进行卷积的问题,但是当我从我的存档中引入FFT模板时,我发现输出存在问题。寻找关于FFT模板的帮助
例如:
输入:(0,0)(0,0)(4166667,0)(1,0)
正确的输出:(4166668,0)(-4166667, 1)(4166666,0)(-4166667,-1)
模板输出:(4166668,0)(-4166667,-1)(4166666,0)(-4166667,)
代码:
#define MAXN
#define ld long double
#define op operator
struct base {
typedef ld T; T re, im;
base() :re(0), im(0) {}
base(T re) :re(re), im(0) {}
base(T re, T im) :re(re), im(im) {}
base op + (const base& o) const { return base(re + o.re, im + o.im); }
base op - (const base& o) const { return base(re - o.re, im - o.im); }
base op * (const base& o) const { return base(re * o.re - im * o.im, re * o.im + im * o.re); }
base op * (ld k) const { return base(re * k, im * k); }
base conj() const { return base(re, -im); }
};
base w[MAXN]; //omega lookup table
int rev[MAXN]; //reverse lookup table
void build_rev(int k) {
static int rk = -1;
if(k == rk)return ; rk = k;
for(int i = 1; i < (1<<k); i++) {
int j = rev[i-1], t = k-1;
while(t >= 0 && ((j>>t)&1)) { j ^= 1 << t; --t; }
if(t >= 0) { j ^= 1 << t; --t; }
rev[i] = j;
}
}
void fft(base *a, int k) {
build_rev(k); int n = 1 << k;
for(int i = 0; i < n; i++) if(rev[i] > i) swap(a[i], a[rev[i]]);
for(int l = 2, lll = 1; l <= n; l += l, lll += lll) {
if(w[lll].re == 0 && w[lll].im == 0) {
ld angle = PI/lll;
base ww(cosl(angle), sinl(angle));
if(lll > 1) for(int j = 0; j < lll; ++j) {
if(j & 1) w[lll + j] = w[(lll+j)/2] * ww;
else w[lll + j] = w[(lll+j)/2];
} else w[lll] = base(1, 0);
}
for(int i = 0; i < n; i += l)
for(int j = 0; j < lll; j++){
base v = a[i + j], u = a[i + j + lll] * w[lll + j];
a[i + j] = v + u; a[i + j + lll] = v - u;
}
}
}
//ideone compiled example: http://ideone.com/8PTjW5
我试图检查位反向统一表的根但我没有发现在这两个部分的任何问题。我也检查了一些在线资料来验证步骤,但没有什么奇怪的东西给我看。
会有人介意帮我找出这个模板有什么问题吗?
在此先感谢。
编辑:我决定在最后依靠另一个模板,感谢所有人的回复。
我怀疑,除了你之外,任何人都可以调试这段代码。它非常密集,没有任何评论。 – OutOfBound
只要您是唯一读过代码的人,使用变量名称如'lll'就没有问题:P – user463035818
我建议将单位向量转换(或反向转换)为只有一个非零零的元素,方式你可能会更容易看到什么系数是关闭的 – user463035818