Introduction

This is the second part of SIMDized sum of all bytes in the array. The first part describes summing unsigned bytes, here we're going to experiment with summing of signed bytes.

The baseline C implementation is:

int32_t sumbytes(int8_t* array, size_t size) {

    int32_t result = 0;

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

    return result;
}

And the C++ implementation:

#include <numeric>

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

vpmovsxbw       %xmm1, %ymm2
vextracti128    $0x1, %ymm1, %xmm1
vpmovsxwd       %xmm2, %ymm3
vextracti128    $0x1, %ymm2, %xmm2
vpmovsxbw       %xmm1, %ymm1
vpaddd          %ymm0, %ymm3, %ymm3
vpmovsxwd       %xmm2, %ymm0
vpaddd          %ymm3, %ymm0, %ymm2
vpmovsxwd       %xmm1, %ymm0
vextracti128    $0x1, %ymm1, %xmm1
vpaddd          %ymm2, %ymm0, %ymm0
vpmovsxwd       %xmm1, %ymm1
vpaddd          %ymm0, %ymm1, %ymm0

Explicit casting

There's instruction VPMOVSXBD which casts directly from int8_t to int32_t. It get lower 64 bits of input register, thus to convert all bytes from an AVX2 register the instruction has to be called 4 times and some shifting is needed.

Below is sample implementation

int32_t avx2_sumsignedbytes(int8_t* array, size_t size) {

    __m256i accumulator = _mm256_setzero_si256();

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

        const __m128i lo = _mm256_extracti128_si256(v, 0);
        const __m128i hi = _mm256_extracti128_si256(v, 1);

        const __m256i t0 = _mm256_cvtepi8_epi32(lo);
        const __m256i t1 = _mm256_cvtepi8_epi32(hi);
        const __m256i t2 = _mm256_cvtepi8_epi32(_mm_bsrli_si128(lo, 8));
        const __m256i t3 = _mm256_cvtepi8_epi32(_mm_bsrli_si128(hi, 8));

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

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

Shift-based casting

In this approach we also cast directly from 8 into 32-bit numbers, but we use 32-bit shifts.

To convert 3rd bytes in all 32-bit words of registers we simply do arithmetic shift right by 24 bits right. This shift repeats the most significant (a sign bit) of our 8-bit value.

However, to do the same for other bytes in a 32-bit we need two shifts. The first one left, which places byte at 3rd position. Then arithmetic shift right is used to extend the type.

Algorithm is:

1. Load the input vector.

   // v   = [  5 | -1 |  2 | -3 |  7 |  1 |  2 |  3 | -6 | -1 | -3 |  8 | -7 | -12|  3 |  2 ]
   const __m256i v = _mm256_loadu_si256(ptr);
   

2. Extend 3rd bytes.

   // v0  = [                 5 |                 7 |                -6 |                -7 ]
   const __m256i v0 = _mm256_srai_epi32(v, 3*8);
   

3. Extend 2nd bytes.

   // v1  = [ -1 |  2 | -3 |  0 |  1 |  2 |  3 |  0 | -1 | -3 |  8 |  0 | -12|  3 |  2 |  0 ] >>> 24
   //     = [                -1 |                 1 |                -1 |               -12 ]
   const __m256i v1 = _mm256_srai_epi32(_mm256_slli_epi32(v, 1*8), 3*8);
   

4. Extend 1st bytes.

   // v2  = [  2 | -3 |  0 |  0 |  2 |  3 |  0 |  0 | -3 |  8 |  0 |  0 |  3 |  2 |  0 |  0 ] >>> 24
   //     = [                 2 |                 2 |                -3 |                 3 ]
   const __m256i v2 = _mm256_srai_epi32(_mm256_slli_epi32(v, 2*8), 3*8);
   

2. Extend 0th bytes.

   // v3  = [ -3 |  0 |  0 |  0 |  3 |  0 |  0 |  0 |  8 |  0 |  0 |  0 |  2 |  0 |  0 |  0 ] >>> 24
   //     = [                -3 |                 3 |                 8 |                 2 ]
   const __m256i v3 = _mm256_srai_epi32(_mm256_slli_epi32(v, 3*8), 3*8);
   

4. Update the accumulator.

   accumulator = _mm256_add_epi32(accumulator, v0);
   accumulator = _mm256_add_epi32(accumulator, v1);
   accumulator = _mm256_add_epi32(accumulator, v2);
   accumulator = _mm256_add_epi32(accumulator, v3);
   

Sample implementation:

int32_t avx2_sumsignedbytes_variant2(int8_t* array, size_t size) {

    __m256i accumulator = _mm256_setzero_si256();

    for (size_t i=0; i < size; i += 32) {
        const __m256i v = _mm256_loadu_si256((__m256i*)(array + i));
        const __m256i v0 = _mm256_srai_epi32(v, 3*8);
        const __m256i v1 = _mm256_srai_epi32(_mm256_slli_epi32(v, 1*8), 3*8);
        const __m256i v2 = _mm256_srai_epi32(_mm256_slli_epi32(v, 2*8), 3*8);
        const __m256i v3 = _mm256_srai_epi32(_mm256_slli_epi32(v, 3*8), 3*8);

        accumulator = _mm256_add_epi32(accumulator, v0);
        accumulator = _mm256_add_epi32(accumulator, v1);
        accumulator = _mm256_add_epi32(accumulator, v2);
        accumulator = _mm256_add_epi32(accumulator, v3);
    }

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

AVX2-specific instruction VPSADBW — naive approach

AVX2 has got 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 VPSADBW yields four numbers.

And while the instruction is perfect for summing unsigned bytes, in case of signed bytes it not that great. VPSADBW has to be used twice: for positive and negative elements of input vector, as it works only on unsigned bytes.

Below is the outline of algorithm:

1. Load the input vector.

   // v   = [  5 | -1 |  2 | -3 |  7 |  1 |  2 |  3 | -6 | -1 | -3 |  8 | -7 | -12|  3 |  2 ]
   const __m256i v = _mm256_loadu_si256(ptr);
   

2. Find mask for negative numbers.

   // m   = [ 00 | ff | 00 | ff | 00 | 00 | 00 | 00 | ff | ff | ff | 00 | ff | ff | 00 | 00 ]
   const __m256i m  = _mm256_cmplt_epi8(v, zero);
   

3. Left positive elements.

   // pos = [  5 |  0 |  2 |  0 |  7 |  1 |  2 |  3 |  0 |  0 |  0 |  8 |  0 |  0 |  3 |  2 ]
   const __m128i pos = _mm256_andnot_si256(m, v);
   

4. Sum the positive elements.

   // t0  = [                14 |                13 |                 8 |                 5 ]
   const __m256i t0 = _mm256_sad_epu8(pos, zero);
   

5. Get the absolute value.

   // va  = [  5 |  1 |  2 |  3 |  7 |  1 |  2 |  3 |  6 |  1 |  3 |  8 |  7 |  12|  3 |  2 ]
   const __m256i va = _mm256_abs_epi8(v);
   

6. Left only the absolute values of negative elements.

   // neg = [  0 |  1 |  0 |  3 |  0 |  0 |  0 |  0 |  6 |  1 |  3 |  0 |  7 |  0 |  0 |  0 ]
   const __m256i neg = _mm256_and_si256(m, va);
   

7. Sum the negative elements.

   // t1  = [                 4 |                 0 |                10 |                 7 ]
   const __m256i t1 = _mm256_sad_epu8(neg, zero);
   

8. Update the accumulators.

   positive = _mm256_add_epi32(positive, t0);
   negative = _mm256_sub_epi32(negative, t1);
   

Below is the actual implementation.

#define _mm256_cmplt_epi8(a, b) _mm256_cmpgt_epi8(b, a)

int32_t avx2_sadbw_sumsignedbytes(int8_t* array, size_t size) {

    const __m256i zero = _mm256_setzero_si256();
    __m256i positive = zero;
    __m256i negative = zero;

    for (size_t i=0; i < size; i += 32) {
        const __m256i v  = _mm256_loadu_si256((__m256i*)(array + i));
        const __m256i m  = _mm256_cmplt_epi8(v, zero);
        const __m256i va = _mm256_abs_epi8(v);

        // sum just positive numbers
        const __m256i t0 = _mm256_sad_epu8(_mm256_andnot_si256(m, v), zero);

        // sum just negative numbers
        const __m256i t1 = _mm256_sad_epu8(_mm256_and_si256(m, va), zero);

        positive = _mm256_add_epi32(positive, t0);
        negative = _mm256_sub_epi32(negative, t1);
    }

    const __m256i accumulator = _mm256_add_epi32(positive, negative);

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

AVX2-specific instruction VPSADBW — faster algorithm

Matthieu Darbois came up with a way batter method that utilizing VPSADBW. The input range of bytes is changed from int8_t range, i.e. -128 .. +128, into uint8_t range by adding 128. This addition sets the input range to 0 .. 255, and thanks to that single invocation of VPSADBW is sufficient to add unsigned bytes. The input adjustment is done with a cheap binary operation, thus the procedure is as fast as its unsigned counterpart.

Due to the input adjustments, the final result is increased by 128 times the number of bytes. This value must be extracted to get the real singed sum.

Below is the actual implementation.

int32_t avx2_sadbw_variant_sumsignedbytes(int8_t* array, size_t size) {

    const __m256i zero   = _mm256_setzero_si256();
    const __m256i addend = _mm256_set1_epi8(-128);
    __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_xor_si256(v, addend);
        const __m256i t1 = _mm256_sad_epu8(t0, zero);

        accumulator = _mm256_add_epi32(accumulator, t1);
    }

    return int32_t(_mm256_extract_epi32(accumulator, 0)) +
           int32_t(_mm256_extract_epi32(accumulator, 2)) +
           int32_t(_mm256_extract_epi32(accumulator, 4)) +
           int32_t(_mm256_extract_epi32(accumulator, 6)) -
           128 * int32_t(size);
}

AVX2-specific instruction VPMADDUBSW

Matthieu Darbois also proposed to use instruction VPMADDUBSW (_mm256_madd_epi16) to perform horizontal additions of bytes.

The instruction itself is more complex, as it multiplies vertically bytes and then adds the pairs of intermediate signed 16-bit results, forming a new vector of 16-bit values. If another vector is filled with 1, then efficiently we add pair of bytes horizontally.

The inner loop of algorithm uses VPMADDUBSW to obtain 16-bit counters which update a 16-bit local accumulator. To avoid overflow, after a fixed number of iterations (128), the local accumulator is extended to 32-bit values and added to the global, 32-bit accumulator.

Below is the actual implementation.

int32_t avx2_maddubs_sumsignedbytes(int8_t* array, size_t size) {

    const __m256i zero   = _mm256_setzero_si256();
    const __m256i one_8  = _mm256_set1_epi8(1);
    const __m256i one_16 = _mm256_set1_epi16(1);
    __m256i accumulator = zero;

    for (size_t j = size / (32 * 128); j != 0; --j) {
        __m256i local_accumulator = zero;
        for (size_t i = 0; i < 128 * 32; i += 32) {
            const __m256i v   = _mm256_loadu_si256((__m256i*)(array + i));
            const __m256i t0 = _mm256_maddubs_epi16(one_8, v);
            local_accumulator = _mm256_add_epi16(local_accumulator, t0);
        }
        const __m256i t1 = _mm256_madd_epi16(one_16, local_accumulator);
        accumulator = _mm256_add_epi32(accumulator, t1);
        array += 128 * 32;
    }

    size_t remainder = size - ((size / (32 * 128)) * (32 * 128));
    if (remainder) {
        __m256i local_accumulator = zero;
        for (size_t i = 0; i < remainder; i += 32) {
            const __m256i v   = _mm256_loadu_si256((__m256i*)(array + i));
            const __m256i t0 = _mm256_maddubs_epi16(one_8, v);
            local_accumulator = _mm256_add_epi16(local_accumulator, t0);
        }
        const __m256i t1 = _mm256_madd_epi16(one_16, local_accumulator);
        accumulator = _mm256_add_epi32(accumulator, t1);
    }

    const __m128i accumulator128 = _mm_add_epi32(_mm256_extracti128_si256(accumulator, 0), _mm256_extracti128_si256(accumulator, 1));

    return int32_t(_mm_extract_epi32(accumulator128, 0)) +
           int32_t(_mm_extract_epi32(accumulator128, 1)) +
           int32_t(_mm_extract_epi32(accumulator128, 2)) +
           int32_t(_mm_extract_epi32(accumulator128, 3));
}

Experiments

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

Conclusions

• On Skylake and SkylakeX architectures VSADBW method is fastest, 2.5 times faster than compiled code.
• On Haswell VPMADDUBSW method significantly outperforms other approaches, it's faster 1.7 times than single VPSADBW method. However, on Skylake and newer architectures both methods are comparable.
• The results from GCC presented in this article before 2019-02-04 were worse for SkylakeX — to my surprise the flag -O3 doesn't enable -funroll-loops. Thanks to Richard Startin and Travis Downs for pointing this out.

Acknowledgements

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

Source code

Source code is available on github.