Vectorized 16-bit addition in Standard C

Question

The idea is to add a given 16-bit number N to each element of an array of 16-bit integers of arbitrary length, taking advantage of 64-bit integer types and instructions to perform the additions 4 at a time. The function accepts an arbitrary nonzero length. As far as I can tell there is no undefined behavior, but I would like to confirm that this is in fact portable.

_Static_assert(~0==-1, "Unsupported sign representation");
#include <stdint.h>
#include <stddef.h>
typedef union Short2 {
    int16_t a[2];
    int32_t i;
} Short2;
typedef union Short4 {
    Short2 a[2];
    int64_t i;
} Short4;
void add16(int16_t *b, const int16_t a, size_t c) {
    // c == 0 causes unpredictable behavior
    if ((uintptr_t)b&2) { // Handle first element if unaligned; remaining elements will be 4 byte aligned
        *b++ += a;
        --c;
    }
    if (c >= 2) {
        Short4 i;
        const int n = a < 0; // If `a` < 0, simply concatenating 4 copies together leads to unexpected results; the positive version must be copied and concatenated instead, and the resulting 64 bit word negated
        i.a[1] = i.a[0] = (Short2){.a = {a-n, a-n}};
        i.i += n;
        if ((uintptr_t)b&4) { // Handle next pair of elements if unaligned; remaining elements will be 8 byte aligned
            (*(Short2 *)b).i += i.a[0].i;
            b += 2;
            c -= 2;
        }
        while (c >= 4) { // Handle all remaining complete batches of 4 elements
            (*(Short4 *)b).i += i.i;
            b += 4;
            c -= 4;
        }
        if (c&2) { // Handle remaining pair of 2 elements if needed
            (*(Short2 *)b).i += i.a[0].i;
            b += 2;
        }
    }
    if (c&1) // Handle remaining element if needed
        *b += a;
}

Code to test the function (for convenience, no need to review this):

#include <stdio.h>
#include <stdlib.h>
int main(void) {
    size_t s = rand()&127;
    int16_t *arr = malloc(s*sizeof(int16_t));
    for (int i = 0; i < s; ++i) {
        printf("%d, ", arr[i] = rand());
    }
    putchar('\n');
    add16(arr, 3, s);
    for (int i = 0; i < s; ++i) {
        printf("%d, ", arr[i]);
    }
    putchar('\n');
    free(arr);
    return 0;
}

Answer 1

You Have an XY-Problem

Let’s take a look at a very naïve loop to do this in C:

void vec_add16( int16_t *const b, const int16_t x, const size_t n ) {
  for ( size_t i = 0; i < n; ++i ) {
    b[i] += x;
  }
}

Compiling this on Clang 15.0.0 with -std=c17 -march=x86-64-v3 -O3 (Godbolt), the inner loop becomes:

.LBB0_8:                                # =>This Inner Loop Header: Depth=1
        vpaddw  ymm1, ymm0, ymmword ptr [rdi + 2*rcx]
        vpaddw  ymm2, ymm0, ymmword ptr [rdi + 2*rcx + 32]
        vpaddw  ymm3, ymm0, ymmword ptr [rdi + 2*rcx + 64]
        vpaddw  ymm4, ymm0, ymmword ptr [rdi + 2*rcx + 96]
        vmovdqu ymmword ptr [rdi + 2*rcx], ymm1
        vmovdqu ymmword ptr [rdi + 2*rcx + 32], ymm2
        vmovdqu ymmword ptr [rdi + 2*rcx + 64], ymm3
        vmovdqu ymmword ptr [rdi + 2*rcx + 96], ymm4
        vpaddw  ymm1, ymm0, ymmword ptr [rdi + 2*rcx + 128]
        vpaddw  ymm2, ymm0, ymmword ptr [rdi + 2*rcx + 160]
        vpaddw  ymm3, ymm0, ymmword ptr [rdi + 2*rcx + 192]
        vpaddw  ymm4, ymm0, ymmword ptr [rdi + 2*rcx + 224]
        vmovdqu ymmword ptr [rdi + 2*rcx + 128], ymm1
        vmovdqu ymmword ptr [rdi + 2*rcx + 160], ymm2
        vmovdqu ymmword ptr [rdi + 2*rcx + 192], ymm3
        vmovdqu ymmword ptr [rdi + 2*rcx + 224], ymm4
        sub     rcx, -128
        add     r9, -2
        jne     .LBB0_8

The compiler vectorizes this code and unrolls the loop sixty-four times per iteration, using the packed 128-bit vector instructions and preloading the next thirty-two elements. In other words, a modern optimizer can see what you’re doing and already optimizes the loop into better code than you would get by tricking it into using a poor-man’s SIMD on 64-bit instructions.

Fortunately, the compiler sees through that, too, mostly ignores what you told it, and optimizes your version to have a nearly-identical inner loop:

.LBB1_10:                               # =>This Inner Loop Header: Depth=1
        vpaddq  ymm1, ymm0, ymmword ptr [rdi + rcx]
        vpaddq  ymm2, ymm0, ymmword ptr [rdi + rcx + 32]
        vpaddq  ymm3, ymm0, ymmword ptr [rdi + rcx + 64]
        vpaddq  ymm4, ymm0, ymmword ptr [rdi + rcx + 96]
        vmovdqu ymmword ptr [rdi + rcx], ymm1
        vmovdqu ymmword ptr [rdi + rcx + 32], ymm2
        vmovdqu ymmword ptr [rdi + rcx + 64], ymm3
        vmovdqu ymmword ptr [rdi + rcx + 96], ymm4
        vpaddq  ymm1, ymm0, ymmword ptr [rdi + rcx + 128]
        vpaddq  ymm2, ymm0, ymmword ptr [rdi + rcx + 160]
        vpaddq  ymm3, ymm0, ymmword ptr [rdi + rcx + 192]
        vpaddq  ymm4, ymm0, ymmword ptr [rdi + rcx + 224]
        vmovdqu ymmword ptr [rdi + rcx + 128], ymm1
        vmovdqu ymmword ptr [rdi + rcx + 160], ymm2
        vmovdqu ymmword ptr [rdi + rcx + 192], ymm3
        vmovdqu ymmword ptr [rdi + rcx + 224], ymm4
        add     rcx, 256
        add     rbx, -2
        jne     .LBB1_10

So, in the end, the hand-tuned version performs only slightly worse: the initialization is longer, but that’s just constant time.

You might get shorter, simpler code to analyze if you also tell the compiler -Os or -fno-unroll-loops.

Although the main lesson here is that this kind of trick only worked in the last century and is pessimal today, there are a few other tips, useful in general, that jump out at me.

Offsets are size_t (or Maybe ptrdiff_t), not int!

On many systems today, int is a signed 32-bit type with a maximum value of 2,147,483,647. There are some very old or embedded systems where it’s only 32,767, and one compiler in the ’90s made int 64-bit but long 32-bit (which the Standards Committee prohibited the next time it met).

Arrays can have more elements than that! Even a 32-bit system could, in theory, fit an array with INT_MAX+1 16-bit elements into memory and overflow the int, counter. This is undefined behavior! Which is bad for you because most C compilers will take that as permission to silently break your program.

If you want a signed type for this, so that the difference of pointers or indices will never wrap around, the type you want is ptrdiff_t from <stddef.h>.

I have no idea why this is such a common mistake of C programmers. It was even more of an obvious bug to use int array indices before the 32-bit era, when it would break on just a few dozen kilobytes of data. It doesn’t seem to go back to K&R, which sometimes used long (You can forgive Dennis Ritchie for not worrying about a computer having more than 2 GiB of memory back in 1973) and sometimes unsigned. Is it just because int has the shortest name?

Beware of Endian-Ness

You assume little-endian code. Since x86 won on the desktop, most other architectures that used to use big-endian order went bi-endian, but this code will most likely still break on a big-endian machine. (Or if someone tried to compile this on a VAX for laughs.)

Use Static Single Assignments Where Possible

It’s rare to want to modify a function argument, although sometimes coders get cute and write a loop that either increments the base pointer or decrements n. Generally, though, you’ll make it harder for yourself to write bugs, and easier for the optimizer, if you make most of your variables const and assign them their unchanging value as soon as they are created.

In this case, the array pointer cannot be a const int16_t*, because it overwrites the input array in place, but it can be int16_t* const because it’s never modified within the function. What’s especially important in a loop like

for ( size_t i = 0; i < n; ++i )

is to make sure that the compiler knows nothing outside the loop will be modifying either i or n. This code accomplishes that by declaring const size_t n and declaring i within the scope of the loop. If the compiler is sure nothing can, it can perform the optimizations we saw above. If, however, either of those variables could “escape” and be used or modified while the loop is running, the compiler suddenly needs to make the loop do what you literally wrote, which is slower than what it could do if it doesn’t ever need the actual values of i or n.

Give Your Variables Better Names

You currently name them a, b, c, and so forth, apparently in the order that you added them to the code. This is a terrible way to name your variables. In my re-implementation, I kept the one-letter names, but at least used letters with somewhat conventional meanings: i for a loop index, b for a base pointer, n for an array size, and x for an input value. Your IDE will autocomplete longer variable names for you anyway, so there’s no reason to avoid them.