C test functions for 2D FFT

Question

I made a function that directly calculates the discrete Fourier transform in dimension two, as well as functions that automatically initialize a 2D array to either something random, or specific basis elements (when viewing multidimensional arrays as tensors, but that's not important). Also an approximate equality checker. This would be used to verify a more sophisticated implementation of the Fourier transform.

I am mainly looking for feedback on C best practices, proper documentation and code safety, but also on test-design. Is this the proper way to write a test unit?

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
#include <time.h>
#include <stdbool.h>
#include <assert.h>

// For a 2D FFT of size N0 x N1, both powers of two.
#define N0 8
#define N1 4

// Initializes weights
void weightInit(double complex* weights, int length)
{
    for (int j = 0; j < length; j++) {
        weights[j] = cexp(-2.0 * M_PI * j / length * I);
    }
}

// Calculates the N0 x N1 2D FFT directly. Input is overwritten by output.
void naiveFFT2D(double complex input[N0][N1])
{
    // We temporarily write the output here. This is disgusting, but
    // the only way to allocate a multidimensional array on the heap.
    double complex (*output)[N1] = malloc( sizeof(double complex[N0][N1]) );
    assert(output); // malloc returns a null pointer when the memory is not available
                    // and this is the only pointer that gets type casted to 0, so
                    // this terminates the program iff there is not sufficient memory.

    // Initialize the weights
    double complex* weightsN0 = malloc(N0 * sizeof(double complex));
    double complex* weightsN1 = malloc(N1 * sizeof(double complex));
    assert(weightsN0);
    assert(weightsN1);

    weightInit(weightsN0, N0);
    weightInit(weightsN1, N1);

    // Compute the FFT, (k0, k1) corresponds to output
    for (int k0 = 0; k0 < N0; k0++) {
        for (int k1 = 0; k1 < N1; k1++) {
            output[k0][k1] = 0;

            for (int j0 = 0; j0 < N0; j0++) {
                for (int j1 = 0; j1 < N1; j1++) {
                    output[k0][k1] += input[j0][j1] * weightsN0[j0 * k0 % N0] * weightsN1[j1 * k1 % N2];
                }
            }
        }
    }

    // Overwite input with output
    for (int k0 = 0; k0 < N0; k0++) {
        for (int k1 = 0; k1 < N1; k1++) {
            input[k0][k1] = output[k0][k1];
        }
    }

    free(weightsN0);
    free(weightsN1);
    free(output);
}

// Prints an N0 x N1 array.
void printArray(double complex arr[N0][N1])
{
    for (int j0 = 0; j0 < N0; j0++) {
        printf("\n");
        for (int j1 = 0; j1 < N1; j1++) {
            printf("%f + %fi | ", creal(arr[j0][j1]), cimag(arr[j0][j1]));
        }
    }
}

// Initializes input to array with a 1 on place (k0, k1) and 0's elsewhere.
void initBasisElement(int k0, int k1, double complex input[N0][N1])
{
    for (int j0 = 0; j0 < N0; j0++) {
        for (int j1 = 0; j1 < N1; j1++) {
            input[j0][j1] = 0.0;
        }
    }

    input[k0][k1] = 1.0;
}

// Initializes input to random array.
void initRandom(double complex input[N0][N1])
{
    srand((unsigned int) clock());

    for (int j0 = 0; j0 < N0; j0++) {
        for (int j1 = 0; j1 < N1; j1++) {
            // a + bi with 0 <= a, b < 1 random doubles
            input[j0][j1] = (double) rand() / (double) RAND_MAX + ((double) rand() / (double) RAND_MAX) * I;
        }
    }
}

// Takes N0 x N1 arrays a and b, and returns true if and only if
// a and b differ less than error at every entry.
bool testEquality(double error, double complex a[N0][N1], double complex b[N0][N1])
{
    bool equal = true;

    for (int j0 = 0; j0 < N0; j0++) {
        for (int j1 = 0; j1 < N1; j1++) {
            if (cabs(a[j0][j1] - b[j0][j1]) > error) {
                equal = false;
            }
        }
    }

    return equal;
}

int main()
{
    double complex (*input)[N1] = malloc( sizeof(double complex[N0][N1]) );;
    initBasisElement(2, 2, input);
    double complex (*b)[N1] = malloc( sizeof(double complex[N0][N1]) );;
    initRandom(b);
    naiveFFT2D(input);
    printArray(input);

    if (testEquality(0.001, input, b)) {
        printf("\n\nEquality!");
    }

    return 0;
}

Answer 1

This is not how to correctly use assert():

double complex (*output)[N1] = malloc( sizeof(double complex[N0][N1]) );
assert(output); // malloc returns a null pointer when the memory is not available
                // and this is the only pointer that gets type casted to 0, so
                // this terminates the program iff there is not sufficient memory.

Asserts are for documenting (and incidentally testing) code invariants, not for runtime checks. Remember that assert() expands to a no-op in production builds!

Instead of claiming that malloc() never fails, we need to write an actual if statement there - ideally, it should return a status value to the caller, who is in a better position to decide whether to terminate the whole program, or to do some other work first (saving the user's data, perhaps).

In passing, we can simplify the size computation, and make it more robust, by using the dereferenced variable as argument to sizeof, rather than having to write a matching type.

It's a bit inconvenient to be stuck with fixed sizes for our array. The usual way to dynamically allocate a 2-dimensional matrix of values is to malloc() storage for rows ✕ columns elements, and then access elements using an index computed as x + y ✕ rows.

I get a lot of Valgrind reports of reading off the end of allocated memory:

==1197== Invalid read of size 8
==1197==    at 0x109484: naiveFFT2D (259984.c:49)
==1197==    by 0x109953: main (259984.c:126)
==1197==  Address 0x4b79800 is 0 bytes after a block of size 64 alloc'd
==1197==    at 0x483877F: malloc (in /usr/lib/x86_64-linux-gnu/valgrind/vgpreload_memcheck-amd64-linux.so)
==1197==    by 0x1092D4: naiveFFT2D (259984.c:35)
==1197==    by 0x109953: main (259984.c:126)
==1197== 
==1197== Invalid read of size 8
==1197==    at 0x109488: naiveFFT2D (259984.c:49)
==1197==    by 0x109953: main (259984.c:126)
==1197==  Address 0x4b79808 is 8 bytes after a block of size 64 alloc'd
==1197==    at 0x483877F: malloc (in /usr/lib/x86_64-linux-gnu/valgrind/vgpreload_memcheck-amd64-linux.so)
==1197==    by 0x1092D4: naiveFFT2D (259984.c:35)
==1197==    by 0x109953: main (259984.c:126)

These should be corrected (I don't think you really intended to index element j0 * k0, for example, but don't immediately see what you actually meant). There's also leakage of memory from main(), which is easy to clean up:

==1197== 512 bytes in 1 blocks are definitely lost in loss record 1 of 2
==1197==    at 0x483877F: malloc (in /usr/lib/x86_64-linux-gnu/valgrind/vgpreload_memcheck-amd64-linux.so)
==1197==    by 0x109913: main (259984.c:122)
==1197== 
==1197== 512 bytes in 1 blocks are definitely lost in loss record 2 of 2
==1197==    at 0x483877F: malloc (in /usr/lib/x86_64-linux-gnu/valgrind/vgpreload_memcheck-amd64-linux.so)
==1197==    by 0x109937: main (259984.c:124)

I don't see why there are empty statements following the allocations in main() - are those ;; just typos?

The test looks totally flawed to me. Why are we comparing against a randomly-populated matrix? That's unlikely to be correct. And why do we always return zero (i.e. success), even when the actual and expected results are different?

If we have a function available that implements the reverse transform, a useful test would be that we can apply that and get back the original input, within rounding error. It's not sufficient as a test of both functions, but it's certainly helpful.