山东淄博以“互联网＋”改革促老工业城市转型发展

Question

I wrote Cooley-Tukey's FFT algorithm in C. And I want to know how safe is this code? Is it just trash? How can I make it more memory safe?

#include <stdio.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.14159265358979323846

void fft(double complex *x, unsigned int N, double complex *result){

    if(N == 1){
        result[0] = x[0];
        return;
    }

    double complex sub_even[N/2];
    double complex sub_odd[N/2];

    for(unsigned int i = 0; i < N/2; i++){
        sub_even[i] = x[2 * i];
        sub_odd[i] = x[2 * i + 1];
    }

    
    double complex *even = (double complex *)malloc(sizeof(double complex) * (N/2));
    double complex *odd = (double complex *)malloc(sizeof(double complex) * (N/2));
    fft(sub_even, N/2, even);
    fft(sub_odd, N/2, odd);


    for(unsigned int k = 0; k < N / 2; k++){
        double complex twiddle = cexp(-2 * I * PI * k / N) * odd[k];
        result[k] = even[k] + twiddle;
        result[k + N/2] = even[k] - twiddle;
    }


    free(even);
    free(odd);

    
    
}

int main(){

    unsigned int N = 8;

    double complex x[N];

    for(unsigned int i = 0; i < N; i++){
        x[i] = i+1;
    }


    double complex *result = (double complex *)malloc(sizeof(double complex) * N);
    fft(x, N, result);


    for(unsigned int i = 0; i < N; i++){
        printf("%.2f + %.8fj\n", creal(result[i]), cimag(result[i]));
    }



    return 0;
}

Answer 1

double complex sub_even[N/2];
double complex sub_odd[N/2];

These are VLAs and they have the usual problem: it's easy to cause a stack overflow this way. No big deal if N=8, but that's not really a typical FFT size.

These arrays are also not necessary, they only change how the data is indexed (with increasingly large gaps) which you can do with arithmetic instead of by physically moving the data.

double complex *even = (double complex *)malloc(sizeof(double complex) * (N/2));
double complex *odd = (double complex *)malloc(sizeof(double complex) * (N/2));

Even more memory, but FFT can be implemented in-place (transforming the input into the output), or using one buffer of size N if you prefer it that way (just one, allocated before the recursion or pseudo-recursion, not a size N allocation at every level of the recursion).

As for safety, what happens if malloc fails? There's no provision for that (but it would even better to not allocate in the first place). What happens if some huge N is passed in such that sizeof(double complex) * (N/2) = 0 or in general wraps to a small number? Can't happen here really, the VLAs would explode first. But once you remove them, it becomes possible - although I think low-relevance because the only way it happens is if N is wrong (the input array cannot be sufficiently large, its size would also have wrapped), whether intentionally or due to a mistake upstream, either way this code is neither responsible for the problem nor able to do anything about it. But with all that said, it may be a better practice to at least not add an additional bad consequence.

cexp(-2 * I * PI * k / N)

Twiddles can be computed significantly cheaper than exp-ing them one by one.

These basic techniques along with transforming the recursion into iteration are shown for example on wikipedia, although it does not show an efficient way to implement the bit-reverse permutation, for which I refer you to Practically efficient methods for performing bit-reversed permutation in C++11 on the x86-64 architecture (you can do it in C too, and on different architectures).

Answer 2

We assume that N at each level is exactly divisible by 2, but we never confirm this. That can lead to users wondering why up to half of their input is ignored.

---

Consider this allocation:

double complex *result = (double complex *)malloc(sizeof(double complex) * N);

This is C, so the void* returned by malloc() can be assigned to result without any cast. Remove the pointless clutter.

Also, it's a good habit to to use sizeof with an expression rather than a type (*result instead of (double complex)) since that can avoid tediously looking up the target's type when it's far removed from the assignment:

double complex *result = malloc(sizeof *result * N);

And don't forget to ensure that result isn't a null pointer before using it (other than as argument to free()).

---

We're using much more memory than we need to, by allocating at every step of the recursion. We have this lovely result array that we don't use until after the recursive call returns, so that can be used as scratch space. Or re-work to allocate a scratch space at top level and pass that through the recursive calls.

Answer 3

Memory

Seem strange to use VLAs and multiple malloc() allocations - and not checking for success (that is unsafe).

Consider instead 1 allocation and add a check:

// Allocate to the size of the reference object, not the type.
// It is easier to review and maintain.
double complex *dc = malloc(sizeof dc[0] * N/2 * 4);
if (dc == NULL) {
  handle_OOM();  // Out of memory
  ...
} else {
  double complex *sub_even = &dc[0];
  double complex *sub_odd = &dc[N/2];
  ...
  double complex *even = &dc[N/2 * 2];
  double complex *odd = &dc[N/2 * 3];
  ...
  free(dc);
}

Code should check if N is not a power of 2.

Pedantic code would check if N is too large.

if (N > SIZE_MAX/2/sizeof dc[0]) {
  Handle_bad_size();
}

Use const and restrict

When able, use const. It makes for safer maintenance, wider applicability and potential optimizations.

When able, use restrict. When refenced data does not overlap, it allows for potential optimizations. It also lets the user to know data should not overlap. As is, x and result could overlap and code works correctly. Yet I suspect later code re-org may not support that. Do you really want the user to be allowed to provide overlapped data?

// void fft(double complex *x, unsigned int N, double complex *result){
void fft(const double complex * restrict x, unsigned int N, double complex *restrict result){

Printing

Avoid %f (fixed point) for general floating point printing. It it verbose when values are large and uninformative when values are small. Consider %g or even %a.

IMO, this improves output as small values do not report like a zero value.

//printf("%.2f + %.8fj\n", creal(result[i]), cimag(result[i]));
printf("%g + %gj\n", creal(result[i]), cimag(result[i]));

Formatting

Save time, use an auto formatter for code.

Sizing

fft() uses unsigned. size_t is the Goldilocks type for array sizing and indexing. Neither too narrow nor too wide for a general purpose fft().

I'd expect fft() calls to pass in a size_t for size. unsigned may narrow a size_t and thus create a pesky compiler warning.

Design

Perhaps allocate the result space?

// void fft(double complex *x, unsigned int N, double complex *result){
double complex *fft_alloc(size_t n, const double complex x[n]) {

Answer 4

Code reformatted without semantic (or syntactic) changes for better readability:

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define PI 3.14159265358979323846

void fft(double complex *x, unsigned int N, double complex *result) {
  if (N == 1) {
    result[0] = x[0];
    return;
  }

  double complex sub_even[N / 2];
  double complex sub_odd[N / 2];

  for (unsigned int i = 0; i < N / 2; i++) {
    sub_even[i] = x[2 * i];
    sub_odd[i] = x[2 * i + 1];
  }

  double complex *even = (double complex *)malloc(sizeof(double complex) * (N / 2));
  double complex *odd = (double complex *)malloc(sizeof(double complex) * (N / 2));
  fft(sub_even, N / 2, even);
  fft(sub_odd, N / 2, odd);

  for (unsigned int k = 0; k < N / 2; k++) {
    double complex twiddle = cexp(-2 * I * PI * k / N) * odd[k];
    result[k] = even[k] + twiddle;
    result[k + N / 2] = even[k] - twiddle;
  }

  free(even);
  free(odd);
}

int main() {
  unsigned int N = 8;

  double complex x[N];

  for (unsigned int i = 0; i < N; i++) {
    x[i] = i + 1;
  }

  double complex *result = (double complex *)malloc(sizeof(double complex) * N);
  fft(x, N, result);

  for (unsigned int i = 0; i < N; i++) {
    printf("%.2f + %.8fj\n", creal(result[i]), cimag(result[i]));
  }

  return 0;
}

Answer 5

You are defining PI and also importing math.h.

M_PI is already defined in math.h, so no need to define it again.