《数值食谱》书中的 FFT 示例导致运行时错误

Question

我正在尝试在C 上实现FFT 算法。我根据《C 中的数字食谱》一书中的函数“four1”编写了一段代码。我知道使用 FFTW 等外部库会更有效，但我只是想尝试将此作为第一种方法。但我在运行时遇到错误。

\n

经过一段时间的尝试调试后，我决定复制书中提供的完全相同的功能，但我仍然遇到同样的问题。问题似乎出在以下命令中：

\n

tempr = wr * data[j] - wi * data[j + 1];\ntempi = wr * data[j + 1] + wi * data[j];\n

\n

和

\n

data[j + 1] = data[i + 1] - tempi;\n\n

\n

j 有时与数组的最后一个索引一样高，因此在索引时不能加 1。

\n

正如我所说，我没有对代码进行任何更改，所以我很惊讶它对我不起作用；这是一本著名的 C 数值方法参考书，我怀疑其中是否有错误。另外，我发现了一些关于相同代码示例的问题，但它们似乎都没有相同的问题（例如，请参见C: Numerical Recipies (FFT) ）。我究竟做错了什么？

\n

这是代码：

\n

#include <iostream>\n#include <stdio.h>\nusing namespace std;\n\n#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr\nvoid four1(double* data, unsigned long nn, int isign)\n{\n    unsigned long n, mmax, m, j, istep, i;\n    double wtemp, wr, wpr, wpi, wi, theta;\n    double tempr, tempi;\n\n    n = nn << 1;\n    j = 1;\n    for (i = 1; i < n; i += 2) {\n        if (j > i) {\n            SWAP(data[j], data[i]);\n            SWAP(data[j + 1], data[i + 1]);\n        }\n        m = n >> 1;\n        while (m >= 2 && j > m) {\n            j -= m;\n            m >>= 1;\n        }\n        j += m;\n    }\n    mmax = 2;\n    while (n > mmax) {\n        istep = mmax << 1;\n        theta = isign * (6.28318530717959 / mmax);\n        wtemp = sin(0.5 * theta);\n        wpr = -2.0 * wtemp * wtemp;\n        wpi = sin(theta);\n        wr = 1.0;\n        wi = 0.0;\n        for (m = 1; m < mmax; m += 2) {\n            for (i = m; i <= n; i += istep) {\n                j = i + mmax;\n\n                tempr = wr * data[j] - wi * data[j + 1];\n                tempi = wr * data[j + 1] + wi * data[j];\n                data[j] = data[i] - tempr;\n                data[j + 1] = data[i + 1] - tempi;\n                data[i] += tempr;\n                data[i + 1] += tempi;\n            }\n            wr = (wtemp = wr) * wpr - wi * wpi + wr;\n            wi = wi * wpr + wtemp * wpi + wi;\n        }\n        mmax = istep;\n    }\n}\n#undef SWAP\n\n\nint main()\n{\n    // Testing with random data\n    double data[] = {1, 1, 2, 0, 1, 3, 4, 0};\n\n    four1(data, 4, 1);\n\n    for (int i = 0; i < 7; i++) {\n        cout << data[i] << " ";\n    }\n\n}\n

\n

Answer 1

C 语言数值食谱的前 2 个版本使用了不常见的（对于 C）约定，即数组是从 1 开始的。（这可能是因为 Fortran（基于 1）版本首先出现，并且在不考虑约定的情况下完成了对 C 的翻译。）

您应该阅读1.2 节 科学计算的一些 C 约定，特别是有关向量和一维数组的段落。除了试图证明他们基于 1 的决策的合理性之外，本节还解释了如何适当地调整指针以匹配他们的代码。

在你的情况下，这应该有效 -

int main()
{
    // Testing with random data
    double data[] = {1, 1, 2, 0, 1, 3, 4, 0};
    double *data1based = data - 1;

    four1(data1based, 4, 1);

    for (int i = 0; i < 7; i++) {
        cout << data[i] << " ";
    }

}

然而，正如 @Some 程序员家伙在评论中提到的，本书提倡的解决方法是未定义的行为，即data1based数组边界之外的点data。

虽然这种方式在实践中效果很好，但另一种非 UB 解决方法是改变你的解释以符合他们的惯例 -

int main()
{
    // Testing with random data
    double data[] = { -1 /*dummy value*/, 1, 1, 2, 0, 1, 3, 4, 0};

    four1(data, 4, 1);

    for (int i = 1; i < 8; i++) {
        cout << data[i] << " ";
    }

}

我会非常警惕这种情况会传染并广泛感染您的代码。

第三版默认认识到了这个“错误”，并且作为支持 C++ 和标准库集合的一部分，转而使用从零开始的数组的 C 和 C++ 约定。