Example on FFT from Numerical Recipes book results in runtime error-CodePudding

I am trying to implement the FFT algorithm on C. I wrote a code based on the function "four1" from the book "Numerical Recipes in C". I know that using external libraries such as FFTW would be more efficient, but I just wanted to try this as a first approach. But I am getting an error at runtime.

After trying to debug for a while, I have decided to copy the exact same function provided in the book, but I still have the same problem. The problem seems to be in the following commands:

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

and

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

the j is sometimes as high as the last index of the array, so you cannot add one when indexing.

As I said, I didn´t change anything from the code, so I am very surprised that it is not working for me; it is a well-known reference for numerical methods in C, and I doubt there are errors in it. Also, I have found some questions regarding the same code example but none of them seemed to have the same issue (see C: Numerical Recipies (FFT), for example). What am I doing wrong?

Here is the code:

#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] << " ";
    }

}

CodePudding user response：

The first 2 editions of Numerical Recipes in C use the unusual (for C) convention that arrays are 1-based. (This was probably because the Fortran (1-based) version came first and the translation to C was done without regard to conventions.)

You should read section 1.2 Some C Conventions for Scientific Computing, specifically the paragraphs on Vectors and One-Dimensional Arrays. As well as trying to justify their 1-based decision, this section does explain how to adapt pointers appropriately to match their code.

In your case, this should work -

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] << " ";
    }

}

However, as @Some programmer dude mentions in the comments the workaround advocated by the book is undefined behaviour as data1based points outside the bounds of the data array.

Whilst this way well work in practice, an alternative and non-UB workaround would be to change your interpretation to match their conventions -

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] << " ";
    }

}

I'd be very wary of this becoming contagious though and infecting your code too widely.

The third edition tacitly recognised this 'mistake' and, as part of supporting C and standard library collections, switched to use the C & C conventions of zero-based arrays.