Пример БПФ из книги Numerical Recipes приводит к ошибке времени выполнения

Question

Пытаюсь реализовать алгоритм БПФ на C. Написал код на основе функции "four1" из книги "Численные рецепты на C". Я знаю, что использование внешних библиотек, таких как FFTW, было бы более эффективным, но я просто хотел попробовать это в качестве первого подхода. Но я получаю ошибку во время выполнения.

После попытки отладки некоторое время я решил скопировать точно такую ​​же функцию, представленную в книге, но у меня все еще та же проблема. Проблема, похоже, в следующих командах:

      tempr = wr * data[j] - wi * data[j + 1];
tempi = wr * data[j + 1] + wi * data[j];

и

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

j иногда достигает последнего индекса массива, поэтому вы не можете добавить его при индексации.

Как я уже сказал, я ничего не менял в коде, поэтому я очень удивлен, что он у меня не работает; это хорошо известный справочник по численным методам в C, и я сомневаюсь, что в нем есть ошибки. Кроме того, я нашел несколько вопросов, касающихся одного и того же примера кода, но ни один из них, похоже, не имел такой же проблемы (см., например, C: Numerical Recipies (FFT)).Что я делаю не так?

Вот код:

      #include <iostream>
#include <stdio.h>
using namespace std;

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void four1(double* data, unsigned long nn, int isign)
{
    unsigned long n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    double tempr, tempi;

    n = nn << 1;
    j = 1;
    for (i = 1; i < n; i += 2) {
        if (j > i) {
            SWAP(data[j], data[i]);
            SWAP(data[j + 1], data[i + 1]);
        }
        m = n >> 1;
        while (m >= 2 && j > m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    mmax = 2;
    while (n > mmax) {
        istep = mmax << 1;
        theta = isign * (6.28318530717959 / mmax);
        wtemp = sin(0.5 * theta);
        wpr = -2.0 * wtemp * wtemp;
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;
        for (m = 1; m < mmax; m += 2) {
            for (i = m; i <= n; i += istep) {
                j = i + mmax;

                tempr = wr * data[j] - wi * data[j + 1];
                tempi = wr * data[j + 1] + wi * data[j];
                data[j] = data[i] - tempr;
                data[j + 1] = data[i + 1] - tempi;
                data[i] += tempr;
                data[i + 1] += tempi;
            }
            wr = (wtemp = wr) * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
        }
        mmax = istep;
    }
}
#undef SWAP


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

    four1(data, 4, 1);

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

}