Introduction

I was curious how GCC vectorizes function that sums bytes from an array. Below is a loop-based implementation.

uint32_t sumbytes(uint8_t* array, size_t size) {

    uint32_t result = 0;

    for (size_t i=0; i < size; i++)
        result += uint32_t(array[i]);

    return result;
}

The same algorithm can be expressed with following C++ code.

#include <numeric>

uint32_t sumbytes(uint8_t* array, size_t size) {
    return std::accumulate(array, array + size, uint32_t(0));
}

When I saw the assembly generated by GCC I was sure that it's possible to make it better and faster. This text summarizes my findings.

I focus solely on Skylake performance and AVX2 code. The sources have got also implementations of SSE procedures and experiments include timings from an older CPU.

vpmovzxbw    %xmm0, %ymm2
vextracti128 $0x1, %ymm0, %xmm0
vpmovzxwd    %xmm2, %ymm1
vextracti128 $0x1, %ymm2, %xmm2
vpmovzxbw    %xmm0, %ymm0
vpmovzxwd    %xmm2, %ymm2
vpaddd       %ymm2, %ymm1, %ymm1
vpmovzxwd    %xmm0, %ymm2
vextracti128 $0x1, %ymm0, %xmm0
vpaddd       %ymm2, %ymm1, %ymm1
vpmovzxwd    %xmm0, %ymm0
vpaddd       %ymm0, %ymm1, %ymm0
vpaddd       %ymm0, %ymm3, %ymm3

32-bit accumulators

In this approach we cast directly from uint8_t to uint32_t. But casting isn't done by the mentioned VPMOVZXBD — we use only bit-ands and shifts within 32-bit lanes, that has some benefits.

1. Load the vector:

   //       |15          13 | 12          8 | 7           4 | 3           0 |
   // v   = [ a | b | c | d | e | f | g | h | i | j | k | l | m | n | o | p ]
   const __m128i v = _mm_loadu_si128(ptr);
   

2. Mask 0th bytes in each 32-bit lane

   const __m128i lobyte_mask = _mm_set1_epi32(0x000000ff);
   
   // t0  = [ 0 | 0 | 0 | d | 0 | 0 | 0 | h | 0 | 0 | 0 | l | 0 | 0 | 0 | p ]
   const __m128i t0 = _mm_and_si128(lobyte_mask, v);
   

3. Move 1st bytes into 0th position in each lane.

   const __m128i v1 = _mm_srli_epi32(v, 1*8);
   
   // t1  = [ 0 | 0 | 0 | c | 0 | 0 | 0 | g | 0 | 0 | 0 | k | 0 | 0 | 0 | o ]
   const __m128i t1 = _mm_and_si128(lobyte_mask, v1);
   

4. Move 2nd bytes into 0th position in each lane.

   const __m128i v2 = _mm_srli_epi32(v, 2*8);
   
   // t2  = [ 0 | 0 | 0 | b | 0 | 0 | 0 | f | 0 | 0 | 0 | j | 0 | 0 | 0 | n ]
   const __m128i t2 = _mm_and_si128(lobyte_mask, v2);
   

5. Move 3rd bytes into 0th position in each lane.

   // t3  = [ 0 | 0 | 0 | a | 0 | 0 | 0 | e | 0 | 0 | 0 | i | 0 | 0 | 0 | m ]
   const __m128i t3 = _mm_srli_epi32(v, 3*8);
   

6. Update the accumulator.

   accumulator = _mm_add_epi32(accumulator, t0);
   accumulator = _mm_add_epi32(accumulator, t1);
   accumulator = _mm_add_epi32(accumulator, t2);
   accumulator = _mm_add_epi32(accumulator, t3);
   

Here is an actual implementation:

uint32_t sse_sumbytes_variant2(uint8_t* array, size_t size) {

    const __m128i lobyte_mask = _mm_set1_epi32(0x000000ff);

    __m128i accumulator = _mm_setzero_si128();

    for (size_t i=0; i < size; i += 16) {
        const __m128i v = _mm_loadu_si128((__m128i*)(array + i));
        const __m128i v0 = v;
        const __m128i v1 = _mm_srli_epi32(v, 1*8);
        const __m128i v2 = _mm_srli_epi32(v, 2*8);
        const __m128i v3 = _mm_srli_epi32(v, 3*8);

        const __m128i t0 = _mm_and_si128(lobyte_mask, v0);
        const __m128i t1 = _mm_and_si128(lobyte_mask, v1);
        const __m128i t2 = _mm_and_si128(lobyte_mask, v2);
        const __m128i t3 = v3;

        accumulator = _mm_add_epi32(accumulator, t0);
        accumulator = _mm_add_epi32(accumulator, t1);
        accumulator = _mm_add_epi32(accumulator, t2);
        accumulator = _mm_add_epi32(accumulator, t3);
    }

    return uint32_t(_mm_extract_epi32(accumulator, 0)) +
           uint32_t(_mm_extract_epi32(accumulator, 1)) +
           uint32_t(_mm_extract_epi32(accumulator, 2)) +
           uint32_t(_mm_extract_epi32(accumulator, 3));
}

The assembly code of the main loop is:

vpsrld $0x8, %ymm1, %ymm3
vpsrld $0x10, %ymm1, %ymm2
vpsrld $0x18, %ymm1, %ymm0
vpand  %ymm3, %ymm4, %ymm3
vpaddd %ymm5, %ymm0, %ymm0
vpand  %ymm2, %ymm4, %ymm2
vpand  %ymm1, %ymm4, %ymm1
vpaddd %ymm3, %ymm0, %ymm0
vpaddd %ymm2, %ymm0, %ymm0
vpaddd %ymm1, %ymm0, %ymm0

This approach is two times faster than the compiler's version. There are two reasons of that:

1. We have 10 instructions instead of 13.
2. All of these instruction are simple and they are fast. Let's look at the instructions' parameters for the Skylake architecture taken from Agner Fog's site.

First of all, most instructions emitted by compiler have latency 3 cycles. But latency could be compensated (especially in unrolled loops). Worse is that all instructions except VPADD are dispatched on the same execution port #5. So, it's likely that there will be congestion on the port.

Instructions used by our approach can be dispatched to different ports (and thus executed in parallel); moreover, all have latency just 1 cycle.

16-bit accumulators

Using 16-bit accumulators instead of 32-bit accumulators doubles capacity, so in AVX2 code we're updating 16 instead of 8 accumulators in a single iteration. But this approach limits the number of additions that can be performed without overflow to max(uint16_t)/max(uint8_t) = 65535/255 = 257.

The overcome this we need two nested loops. The outer loop keeps 32-bit accumulators, as the initial version does. The inner loop processes 256 * 8 bytes, i.e. iterates 64 times. When the inner loop completes, 16-bit accumulators are extended to 32-bit and the outer accumulators are updated.

uint32_t avx2_16bit_sumbytes_variant2(uint8_t* array, size_t size) {

    const __m256i mask_lobytes = _mm256_set1_epi16(0x00ff);

    __m256i accumulator = _mm256_setzero_si256();

#define SIZE (8 * 256)
    for (size_t i=0; i < size; i += SIZE) {

        __m256i accu_16bit = _mm256_setzero_si256();

        for (size_t j=0; j < SIZE; j += 32) {
            const __m256i v = _mm256_loadu_si256((__m256i*)(array + i + j));

            const __m256i t0 = _mm256_and_si256(v, mask_lobytes);
            const __m256i t1 = _mm256_srli_epi16(v, 8);

            accu_16bit = _mm256_add_epi16(accu_16bit, t0);
            accu_16bit = _mm256_add_epi16(accu_16bit, t1);
        }

        const __m256i t0 = _mm256_and_si256(accu_16bit, _mm256_set1_epi32(0x0000ffff));
        const __m256i t1 = _mm256_srli_epi32(accu_16bit, 16);

        accumulator = _mm256_add_epi32(accumulator, t0);
        accumulator = _mm256_add_epi32(accumulator, t1);
    }
#undef SIZE

    return uint32_t(_mm256_extract_epi32(accumulator, 0)) +
           uint32_t(_mm256_extract_epi32(accumulator, 1)) +
           uint32_t(_mm256_extract_epi32(accumulator, 2)) +
           uint32_t(_mm256_extract_epi32(accumulator, 3)) +
           uint32_t(_mm256_extract_epi32(accumulator, 4)) +
           uint32_t(_mm256_extract_epi32(accumulator, 5)) +
           uint32_t(_mm256_extract_epi32(accumulator, 6)) +
           uint32_t(_mm256_extract_epi32(accumulator, 7));
}

acc_lo += input;
if (carry_flag)
    acc_hi += 1;

uint8_t satadd(uint8_t a, uint8_t b) {

    if (uint32_t(a) + uint32_t(b) <= 255) // this check is done *somehow* in the hardware
        return a + b;
    else
        return 255;
}

acc_lo = 0;
acc_hi = 64;

for (size_t i=0; i < 64; i++) {
    const uint8_t tmp = satadd(acc_lo, input);
    acc_lo += input;

    if (acc_lo == tmp) {
        acc_hi -= 1;
    }
}

uint32_t sse_8bit_sumbytes(uint8_t* array, size_t size) {

    __m128i accumulator = _mm_setzero_si128();

    for (size_t i=0; i < size; i += 256) {

        __m128i accu_8bit_lo = _mm_setzero_si128();
        __m128i accu_8bit_hi = _mm_set1_epi8(16);

        for (size_t j=0; j < 256; j += 16) {
            const __m128i v = _mm_loadu_si128((__m128i*)(array + i + j));

            const __m128i t0 = _mm_add_epi8(v, accu_8bit_lo);
            const __m128i t1 = _mm_adds_epu8(v, accu_8bit_lo);

            const __m128i not_carry = _mm_cmpeq_epi8(t0, t1);

            accu_8bit_lo = t0;
            accu_8bit_hi = _mm_add_epi8(accu_8bit_hi, not_carry);
        }

        const __m128i sum_lo = _mm_sad_epu8(accu_8bit_lo, _mm_setzero_si128());
        const __m128i sum_hi = _mm_slli_epi32(_mm_sad_epu8(accu_8bit_hi, _mm_setzero_si128()), 8);

        accumulator = _mm_add_epi32(accumulator, sum_lo);
        accumulator = _mm_add_epi32(accumulator, sum_hi);
    }

    return uint32_t(_mm_extract_epi32(accumulator, 0)) +
           uint32_t(_mm_extract_epi32(accumulator, 1)) +
           uint32_t(_mm_extract_epi32(accumulator, 2)) +
           uint32_t(_mm_extract_epi32(accumulator, 3));
}

AVX2-specific instruction VPSADBW

AVX2 defines instruction VPSADBW (_mm256_sad_epu8) that calculates Sum Of Absolute Differences (SAD) of unsigned bytes. Single SAD function works on eight-element subvectors (64-bit slices), and stores the results on corresponding 64-bit elements of the result vector; in case of AVX2 the VPSADB yields four numbers.

Below is a sample C++ code that naively implements SAD function.

uint64_t SAD(uint8_t a[], uint8_t b[]) {

    uint64_t result = 0;

    for (i=0; i < 8; i++) {
        if (a[i] > b[i])
            result += a[i] - b[i];
        else
            result += b[i] - a[i];
    }

    return reuslt;
}

If one of argument vectors is full of zeros, then SAD calculates the sum of bytes from another argument — which exactly what we need. Moreover, since the instruction calculates 64-bit sums, bigger arrays can be processed without overflow.

Below is the simplest, non-unrolled C++ code.

uint32_t avx2_sadbw_sumbytes(uint8_t* array, size_t size) {

    const __m256i zero = _mm256_setzero_si256();
    __m256i accumulator = zero;

    for (size_t i=0; i < size; i += 32) {
        const __m256i v  = _mm256_loadu_si256((__m256i*)(array + i));
        const __m256i t0 = _mm256_sad_epu8(v, zero);

        accumulator = _mm256_add_epi32(accumulator, t0);
    }

    return uint32_t(_mm256_extract_epi32(accumulator, 0)) +
           uint32_t(_mm256_extract_epi32(accumulator, 2)) +
           uint32_t(_mm256_extract_epi32(accumulator, 4)) +
           uint32_t(_mm256_extract_epi32(accumulator, 6));
}

Although the VPSADBW instruction has latency of 3 cycles, it has low throughput of 1 cycle. In case of an unrolled procedure the latency is almost ridden of. The performance boost over the compiler code is almost 6 times.

Experiments

The procedures were run three times and minimum values were noted.

Conclusions

• On Skylake and SkylakeX architectures SADBW method is unbeatable, it's 5-6 times faster than compiled code.
• On Westmere we can be even faster than the SADB method (but not too much). Anyways, we're faster 2-3 times than compiled code.
• It's worth to note that GCC compiles plain C and C++ procedures into the same code. It is another proof that C++ templates are not that bad (of course in term of performance).

See also

There is also second part, dealing with signed bytes

Acknowledgements

Big thanks to Daniel Lemire who provides access to Skylake/SkylakeX machines, where I can run benchmarks.

Source code

Source code is available on github.