Base64 encoding with SIMD instructions

Base64 decoding could also be vectorized, although the speedup is not very impressive, merely 35%. Read more ...

Base64 encoding with SIMD instructions

An SSE code is more than 2 times faster on Core i7, and around 70% faster on Core i5. Read more...

SIMD-ized searching in unique constant dictionary

The problem: there is a ordered dictionary containing only unique keys. Dictionary is read only, and keys are 32-bit (SSE) or 64-bit (AVX2). Read more

SIMD: detecting a bit pattern

The problem: there are 64-bit values with some data bits and some metadata bits; metadata includes a k-bit field describing a "type" (k >= 0). Type field is located in a lower 32-bits.

Procedure processes two "types", one denoted with code 3 and another with 5. When all items are of type 3 then we can use a fast AVX2 path, if there are some types 5, we have to call an additional function (a virtual method, to be precise). Read more ...

AVX512: ternary functions evaluation

Intel's version of SIMD offers following 2-argument (binary) boolean functions: and, or, xor, and not. There isn't a single argument not, this function can be expressed with xor reg, ones, however this require additional, pre-set register.

AVX512F will come with very interesting instruction called vpternlog. Read more ...

Not everything in AVX2 is 256-bit

AVX2 has added support for 256-bit arguments for many operations on packed integers, although not all. Some instructions accept the 256-bit registers, but operates on 128-bit lanes rather whole register.

There are three major groups of instructions: packing (narrowing conversion), unpacking (interleave) and permutations; below is a full list of instructions (with intrinsics):

• valignr (_mm256_alignr_epi8)
• vpslldq (_mm256_bslli_epi128)
• vpsrldq (_mm256_bsrli_epi128)
• vmpsadbw (_mm256_mpsadbw_epu8)
• vpacksswb (_mm256_packs_epi16)
• vpackssdw (_mm256_packs_epi32)
• vpackuswb (_mm256_packus_epi16)
• vpackusdw (_mm256_packus_epi32)
• vperm2i128 (_mm256_permute2x128_si256)
• vpermq (_mm256_permute4x64_epi64)
• vpermpd (_mm256_permute4x64_pd)
• vpshufd (_mm256_shuffle_epi32)
• vpshufb (_mm256_shuffle_epi8)
• vpshufhw (_mm256_shufflehi_epi16)
• vpshuflw (_mm256_shufflelo_epi16)
• vpslldq (_mm256_slli_si256)
• vpsrldq (_mm256_srli_si256)
• vpunpckhwd (_mm256_unpackhi_epi16)
• vpunpckhdq (_mm256_unpackhi_epi32)
• vpunpckhqdq (_mm256_unpackhi_epi64)
• vpunpckhbw (_mm256_unpackhi_epi8)
• vpunpcklwd (_mm256_unpacklo_epi16)
• vpunpckldq (_mm256_unpacklo_epi32)
• vpunpcklqdq (_mm256_unpacklo_epi64)
• vpunpcklbw (_mm256_unpacklo_epi8)

For me the most surprising are packing instructions (vpack*) as they require additional shuffling (after or before the instruction) if we want to keep order of values. In some cases the order is crucial.

SSE: Generating mask where n leading (trailing) bytes are set

Informal specification:

__m128i mask_lower(const unsigned n) {
    
    assert(n < 16);
    switch (n) {
        case 0: return {0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
        case 1: return {0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
        case 2: return {0xff, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
        // ...
        case 14: return {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x00};
        case 15: return {0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};
    }
}

__m128i mask_higher(const unsigned n) {
    
    assert(n < 16);
    return ~mask_lower(15 - n);
}

Read more ...

Conversion int and double to string - comparison

Milo Yip has compared different itoa and dtoa implementations on Core i7, including my itoa algorithm 2, that use SSE2 instructions.

Results for itoa are interesting: SSE2 version is not as good as it seemed to be. Tricky branchlut algorithm is only 10% slower, moreover is perfectly portable. One obvious drawback of this method is using lookup-table - in real environment where is a big pressure on cache, memory access could be a bottleneck.

Speeding up searching in linked list

Sounds crazy, but it's possible in some cases. Here are experiments results - 3 times faster isn't so bad.

list                          :      0.780s, speedup  1.00
array list (4)                :      0.703s, speedup  1.11
array list (8)                :      0.515s, speedup  1.51
SIMD array list (4)           :      0.365s, speedup  2.14
SIMD array list (8)           :      0.258s, speedup  3.03

Conversion number to hexadecimal representation

Conversion numbers to hexadecimal representation - SWAR, plain SSE, and draft of BMI2 implementation.

Article SSSE3: printing hex values describes the same topic but is limited to exploit PSHUFB.

Conversion numbers to binary representation

New article Conversion numbers to binary representation — SIMD & SWAR versions.

Few years ago I've described an MMX variant of SIMD algorithm. But the text was in Polish, so audience was limited.

Scalar version of SSE move mask instruction

SSE instruction PMOVMSKB gathers all most significant bits from bytes and stores them as a single 16-bit value; similar action is performed by MOVMSKPD and MOVMSKPS.

Such operation could be easily done using scalar multiplication. Read more ...

SIMD-friendly Rabin-Karp modification

Rabin-Karp algorithm uses a weak hash function to locate possible substring positions. This modification uses merely equality of first and last char of searched substring, however equality of chars can be done very fast in parallel, even without SIMD instruction. Read more ...

Integer log 10 of unsigned integer - SIMD version

Fast calculate ceil(log10(x)) of some unsigned number is described on Bit Twiddling Hacks, this text show the SIMD solution for 32-bit numbers.

Algorithm:

1. populate value in XMM registers. Since maximum value of this function is 10 we need three registers.

movd   %eax, %xmm0          // xmm0 = packed_dword(0, 0, 0, x)
pshufd $0, %xmm0, %xmm0 \n" // xmm0 = packed_dword(x, x, x, x)
movapd %xmm0, %xmm1
movapd %xmm0, %xmm2

2. compare these numbers with sequence of powers of 10.

// powers_a = packed_dword(10^1 - 1, 10^2 - 1, 10^3 - 1, 10^4 - 1)
// powers_c = packed_dword(10^5 - 1, 10^6 - 1, 10^7 - 1, 10^8 - 1)
// powers_c = packed_dword(10^9 - 1, 0, 0, 0)
pcmpgtd powers_a, %xmm0
pcmpgtd powers_b, %xmm1
pcmpgtd powers_c, %xmm2

result of comparisons are: 0 (false) or -1 (true), for example:

xmm0 = packed_dword(-1, -1, -1, -1)
xmm1 = packed_dword( 0, 0, -1, -1)
xmm2 = packed_dword( 0, 0, 0, 0)

3. calculate sum of all dwords

psrld $31, %xmm0       // xmm0 = packed_dword( 1, 1, 1, 1) - convert -1 to 1
psubd %xmm1, %xmm0     // xmm0 = packed_dword( 1, 1, 2, 2)
psubd %xmm2, %xmm0     // xmm0 = packed_dword( 1, 1, 2, 2)

// convert packed_dword to packed_word
pxor %xmm1, %xmm1
packssdw %xmm1, %xmm0 // xmm0 = packed_word(0, 0, 0, 0, 1, 1, 2, 2)

// max value of word in xmm0 is 3, so higher
// bytes are always zero
psadbw %xmm1, %xmm0   // xmm0 = packded_qword(0, 6)

4. save result, i.e. the lowest dword

movd %xmm0, %eax      // eax = 6

Sample program is available.

Mask for zero/non-zero bytes

The description of Determine if a word has a zero byte from "Bit Twiddling Hacks" says about haszero(v): "the result is the high bits set where the bytes in v were zero".

Unfortunately this is not true. High bits are also set for ones followed zeros, i.e. haszero(0xff010100) = 0x00808080. Of course the result is still valid (non-zero if there were any zero byte), but if we want to iterate over all zeros or find last zero index, this could be a problem.

It's possible to create exact mask:

uint32_t nonzeromask(const uint32_t v) {
    // MSB are set if any of 7 lowest bits are set
    const uint32_t nonzero_7bit = (v & 0x7f7f7f7f) + 0x7f7f7f7f;

    return (v | nonzero_7bit) & 0x80808080;
}

uint32_t zeromask(const uint32_t v) {
     // negate MSBs
     return (nonzeromask(v)) ^ 0x80808080;
}

Function nonzeromask requires 4 simple instructions, and zeromask one additional xor.

Extensions to x86 ISA are useless

Intel announced new extension to SSE: instructions accelerating calculating hashes SHA-1 and SHA256. As everything else added recently to x86 ISA, these new instructions address special cases of "something". Number of instructions, encoding modes, etc. is increasing, but do not help in general.

Let see what sha1msg1 xmm1, xmm2 does (type of arguments is packed dword):

 result[0] := xmm1[0] xor xmm1[2]
 result[1] := xmm1[1] xor xmm1[3]
 result[2] := xmm1[2] xor xmm2[0]
 result[3] := xmm1[3] xor xmm2[1]

1. Logical operation "xor" is hardcoded. Why we can't use "or", "and", "not and"? These operations are already present in ISA.
2. Indices to xmm1 and xmm2 are hardcoded too. Instruction pshufd accepts immediate argument (1 byte) to select permutation, why sha1msg1 couldn't be feed with 2 bytes allowing programmer to select any permutations of arguments?
3. Sources of operators are also hardcoded. Why not use another immediate (1 byte) to select sources, for example 00b = xmm1/xmm1, 01b = xmm1/xmm2, 10b = xmm2/xmm1, 11b = xmm2/xmm2.

Such generic instruction would be saved as generic_op xmm1, xmm2, imm_1, imm_2, imm_3 and execute following algorithm:

 for i := 0 to 3 do
  arg1_indice := imm_1[2*i:2*i + 1]
  arg2_indice := imm_2[2*i:2*i + 1]

  if imm_3[2*i] = 1 then
   arg1 := xmm1
  else
   arg1 := xmm2
  end if

  if imm_3[2*i + 1] = 1 then
   arg2 := xmm2
  else
   arg2 := xmm1
  end if

  result[i] := arg1[arg1_indice] op arg2[arg2_indice]
 end for

Then sha1msg1 is just a special case:

 generic_xor xmm1, xmm2, 0b11100100, 0b01001110, 0b01010000

Maybe this example is "too generic", too complex, and would be hard to express in hardware. I just wanted to show that we will get shine new instructions useful in few cases. Compilers can vectorize loops and make use of SSE, but SHA is used in drivers, OS and is encapsulated in libraries --- sha1msg1 and friends will never appear in ordinary programs.

Set of great articles

Lockless Inc publish a lot of interesting, advanced articles. There are some low-level things, algorithms, threading, and many more. Definitely worth to read!

Branchless set mask if value greater or how to print hex values

Suppose we need to get mask when nonnegative argument is greater then some constant value; in other words, we want to evaluate following expression:

if x > const_n then
   mask := 0xffffffff;
else
   mask := 0x00000000;

Portable branchless solution:

• choose magic number M := (1 << (k-1)) - 1 - n, where k is a bit position, for example 31 if we operate on 32-bit words
• calculate R := x + M
• k-th bit of R is set if x > n
• fill mask with this bit - see note Fill word with selected bit

The key to understand this trick is binary form of M: 0111..1111zzzz, where z is 0 or 1 depending on n value. When x is greater then n, then x + M has form 1000..000zzzz, because carry bit propagate through series of ones to k-th position of result.

Real world example - branchless converting hex digit to ASCII (M=0x7ffffff6 for k=31 and n=9).

; input:    eax - hex digit
; output:   eax - ASCII letter (0-9, A-F or a-f)
; destroys: ebx

        andl 0xf, %eax
        leal 0x7ffffff6(%eax), %ebx     ; MSB(ebx)=1 when eax >= 10
        sarl $31, %ebx                  ; ebx - mask
        andl  $7, %ebx                  ; ebx = 7 when eax >= 10 (for A-F letters)
        ;andl $39, %ebx                 ; ebx = 39 when eax >= 10 (for a-f letters)
        leal '0'(%eax, %ebx), %eax      ; eax = '0' + eax + ebx => ASCII letter

It is also possible to convert 4 hex digits in parallel using similar algorithm, but input data have to be correctly prepared. Moreover generating mask requires 3 instructon and one extra register (in scalar version just one arithmetic shift). I guess it wont be fast on x86, maybe this approach would be good for SIMD code, where similar code transforms more bytes at once.

; input: eax - four hex digits in form [0a0b0c0d]
; output: eax - four ascii letters
; destroys: ebx, ecx

        leal 0x76767676(%eax), %ebx        ; MSB of each byte is set when corresponding eax byte is >= 10
                                           ; (here: 0x7f - 9 = 0x76)
        andl $0x80808080, %ebx
        movl %ebx, %ecx
        shrl    $7, %ebx
        subl %ebx, %ecx                    ; ecx - byte-wise mask
        ;andl $0x07070707, %ecx            ; for ASCII letters A-F
        andl $0x27272727, %ecx             ; for ASCII letters a-f
        leal 0x30303030(%eax, %ecx), %eax  ; ecx - four ascii letters

See also: SSSE3: printing hex values (weird use of PSHUFB instruction)