2024-05-08

Vectors in C: an all-you-can-eat buffet of so-so choices

Let's say you want to model a physical system. Often you'll want to represent spacial vectors (members of the space $\mathbb{R}^N$) somehow. In fact there are many use case for Cartesian 2-, 3-, or 4-vectors, and for Lorentz 4-vectors as well as less common uses for other dimensionalities, but let's stick with Cartesian 3-vectos for the sake of definiteness because I want to focus on the programming tradeoffs you face if you want to use C (or a C interface1) for this purpose. Nor do we want to worry about manually optimising for the SIMD module of our chips (compilers are smart these days) or worse still laying things out for the benefit of the GPU.2

Aside: I mostly program in C++ where there are some better options, but I get to mess with a lot of legacy code, so the consequences of someone else making this choice for a code-base I work on are still with me. I might get around to a followup post on ways to adapt an existing legacy library for cleaner inter-operation with new C++ code, but that's for another day: first we must understand the root problem.

C offers us an obvious choice with two painful drawbacks (or perhaps it's one underlying drawback that rear's it's ugly head in two contexts) and a clever way to avoid that issue at the cost of having your soul slowly nibbled to death by syntactic ducks. Nice, huh?

Arrays

The obvious choice is the built-in array type: double vec[3];, right?

The underlying problem is that array are only sort-of first-class types. Consider this code:

#include <stdio.h>

void pass_ptr(double* p)
{
    printf("Passed pointer: %lu\n", sizeof(p)/sizeof(double));
}

void pass_ary(double a[3])
{
    printf("Passed 'array': %lu\n", sizeof(a)/sizeof(double));
}

void pass_c99(double a[static 3]) // Syntax added in c99
{
    printf("Staticly sized: %lu\n", sizeof(a)/sizeof(double));
}

int main()
{
    double vec[3];
    
    printf("In local scope: %lu\n", sizeof(vec)/sizeof(double));

    pass_ptr(vec);
    pass_ary(vec);
    
    return 0;
}
Each of the printf statements is executed once on the same variable, but they generate two results: one says 3 and the others say 1.

This is a classic trap for C newbies. Arrays are not, as is sometimes said, "just pointers" because the symbol table knows how big they are when declared at static or automatic scope. But most things that you can do with them drop that knowledge at which point all that is left is a pointer to the start.3

The other manifestation of the limitation is that you can't assign or perform operations on arrays. That is, this is not legal code:

double v1[3] = {1, 2, 3}
      double v2[3];
      v2 = v3;       // Error! Even when the compiler *does* know the sizes!
and as a result you end up writing your library functions with signatures like cross(double *result, const double *v1, const double *v2) which makes you write particualrly clunky code to use the library:
double v1[3] = {1, 2, 3}
      double v2[3] = {4, 5, 6};
      double cp[3];
      crossProduct(cp, v1, v2);
Ugh. And you get to do it over and over again.

Array-in-struct

Oddly it is amazingly easy to solve both these problems: you just wrap the array declaration in a structure declaration (and typedef it for convenience):

typedef struct {
    double a[3];  // a for "array"
} vector;
This takes up exactly as much memory as before, still knows how many elements are involved when you pass it to a function, and can be assigned! Yeah! It's like magic.

Of course, to access an elements you now write vec.a[2] instead of vec[2], but that's a small price to pay. Right? It's not like your soul will die a little bit each time or anything.

Seasoning with unions

Once you're drunk that CoolAid there is no reason not to go a little further: maybe sometimes you'll want to talk about the coordinates of these things, right? So you make that possible, too:4

typedef struct {
    union {
        double a[3];               // a for "array"
        struct {double x, y, z} c; // c for "coordinate"
    }
} vector;

You still need to write a full set of library routines, but now they can have signatures like5 vector cross(const vector v1, const vector v2) and you can call them like double v1[3] = {1, 2, 3} double v2[3] = {4, 5, 6}; double cp[3] = crossProduct(v1, v2); which is much better.


1 Even if you're confident that you won't be writing in C, you may find it necessary to deal with the limits of the language while planning a binary interface.

2 Those are great tools and worthy of your attention if you have a computationaly demanding task, but beyond the scope of this post.

3 Note that as far as the compiler is concered the declaration of pass_ary is identical to that of pass_ptr: the array-like notation is accepted as syntactic sugar only and the compiler pays no attention to the 3. The _c99 variant is a little more subtle but it still doesn't really know the size of the array, it just assumes a minimum for the sake of optimization (and compilers can complain if static analysis show the assumption is violated). Some folks like to use the array form because it makes the declarion express intent to the human reader (though, like comments, it can lie). Others are not so enamored of it, with Linus famously coming out strongly against it.

4 Type punning with unions this way is strictly forbidden in C++ (because lifetime-model and invariant-enforcement, that's why; and don't give me any guff about POD types either, sonny, the divine gave you memcpy and std::bit_cast for a reason!), but C programmers are rugged, self-reliant individualists who carry Colt-45 six pointers (colt45****** shootin_iron;) on their hips and ain't afraid of no Endian no how.

5 Of course, you might pass pointers for the in-parameters to save a little copying at the cost of writing & all over the place. C sure seems have that same issue come up a lot, eh?.

No comments:

Post a Comment