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.

MSVC 2013 Release code

Today I had a quite long session with debugger and release code, you know: no debugger symbols and optimized code. I spotted two pieces of assembly that forced me to check if compilation mode was really set to release.

00000000000000D9  mov         byte ptr [rcx+1Fh],dl 
00000000000000DC  mov         byte ptr [rdx+rcx],0 
00000000000000E0  movzx       edx,byte ptr [rcx+1Fh] 
00000000000000E4  movzx       eax,dl 

First dl is stored in memory. It's ok. Then edx is used as an offset. Also ok. Seems that compiler knows that highest bits of edx are zeros i.e. edx = dl. Not so fast - edx is reloaded with it's original value and then eax is populated with the same value. These two movzx could be replaced with single mov eax, edx.

And another:

00000000000000BF  movaps      xmm6,xmmword ptr [foo] 
00000000000000C4  movdqa      xmmword ptr [foo],xmm6 

Yes, load & store the same value. Completely useless! Moreover, xmm6 isn't used in following code. (It's worth noting that load is done by an FP-unit and store by integer unit, inter-unit transfers cost one additional cycle on some processor.)

Above instruction sequences were produced by the newest compiler from Microsoft (MSVC 2013, update 3).

Penalties of errors in SSE floating point calculations

SSE provides not widely known control register, called MXCSR. This register plays three roles:

1. controls calculations:
• flag "flush to zero" (described later)
• flag "denormals are zeros" (described later)
• rounding mode (not covered in this text)
2. allow to mask/unmask floating-point exceptions
3. save information about floating-point errors - these flags are sticky, i.e. the programmer is responsible for clearing them.

Floating point errors may cause significant slowdown, some of the flags can silence errors makes a program faster. Read more ...

x86 - ISA where 80% of instructons are unimportant

Few years ago I counted instructions from many Linux binaries --- 2014 is good year to repeated this experiment and see that nothing has changed.

I use 32-bit Debian, my installation has been updated few months ago. All files from /usr/bin and all *.so files from /usr/lib was disassembled with objdump (5050 files were processed). Instructions were grouped simply by mnemonic name, taking into account all addressing and encoding modes would be overkill. I've published script that does the job.

Short summary

• Number of distinct decoded instructions is around 650. Since objdump use AT&T syntax, same opcode is seen under different mnemonics, for example mov is saved as movw, movb, movl depending on argument size.
• Total number of x86 instructions is around 750. Read: one hundred instructions never appeared in the binaries.
• There are 81 instructions used just once. For example quite useful CMPPD.
• There are 22 instructions used twice. For example MFENCE --- no one aware of memory ordering?
• There are 15 instructions used three times. For example BTC, but bit manipulating operations are useless.
• 81 plus 22 plus 15 is 118. Another hundred of useless stuff.

Lets look at top 15 rows from detailed results, i.e. instruction with frequency grater than 1%:

• The total count of these instructions is 87.84% of all instructions (almost all, isn't it?).
• The most frequently used instruction is data transfer (mov/movl) --- 42%
• Control flow instructions (call/ret/jmp) --- 13%.
• Conditions (cmp/test/condition jumps: je/jne) --- 10%.
• Basic arithmetic (add/sub/lea) --- 12%
• Simple stack operations (push/pop) --- 6%

Very interesting observation is that conditions are mostly based on je/jne, i.e. jump if zero/jump if not zero.

First FPU instruction appear at 28-th position. First integer SSE appear at 167-th position. First SSE instruction operating on packed floats appear at 315-th position.

Detailed results

Whole table as txt file.

instruction	count	%
mov	5934098	37.63%
call	1414355	8.97%
lea	1071501	6.79%
movl	760677	4.82%
push	655921	4.16%
jmp	611540	3.88%
add	560517	3.55%
je	490250	3.11%
test	475899	3.02%
pop	441608	2.80%
sub	366228	2.32%
cmp	326379	2.07%
jne	264110	1.67%
nop	242356	1.54%
ret	238569	1.51%
xor	148194	0.94%
movzbl	122730	0.78%
and	88863	0.56%
xchg	66885	0.42%
cmpl	64907	0.41%
movzwl	64589	0.41%
movb	57247	0.36%
or	52138	0.33%
shl	50908	0.32%
cmpb	50152	0.32%
jle	41083	0.26%
leave	39923	0.25%
fldl	37428	0.24%
fstpl	37368	0.24%
shr	36503	0.23%
jbe	32866	0.21%
ja	32333	0.21%
sar	30917	0.20%
flds	29672	0.19%
subl	27636	0.18%
setne	27626	0.18%
testb	27420	0.17%
addl	25906	0.16%
imul	25569	0.16%
jg	24796	0.16%
fstp	24349	0.15%
fxch	23464	0.15%
js	21550	0.14%
fstps	21248	0.13%
sbb	16607	0.11%
inc	16200	0.10%
lock	16049	0.10%
jae	14825	0.09%
sahf	14765	0.09%
dec	14276	0.09%
fnstsw	14026	0.09%
sete	13902	0.09%
movw	13895	0.09%
adc	13640	0.09%
jb	12467	0.08%
jl	11700	0.07%
repz	11178	0.07%
fldcw	11110	0.07%
jge	11019	0.07%
movswl	10816	0.07%
fildl	8852	0.06%
cmpw	7601	0.05%
jns	7490	0.05%
fldz	7331	0.05%
fmul	7229	0.05%
out	7203	0.05%
not	7028	0.04%
movsbl	6720	0.04%
in	6503	0.04%
fld	6309	0.04%
faddp	6254	0.04%
fstl	5760	0.04%
fucom	5753	0.04%
neg	5725	0.04%
fucompp	5354	0.03%
rep	5059	0.03%
fmuls	5039	0.03%
pushl	4430	0.03%
jp	4424	0.03%
fnstcw	4400	0.03%
fld1	4176	0.03%
fmulp	4133	0.03%
orl	3927	0.02%
fadds	3789	0.02%
movq	3779	0.02%
fistpl	3709	0.02%
cltd	3597	0.02%
fmull	3313	0.02%
stos	3298	0.02%
lret	3183	0.02%
scas	3103	0.02%
lods	3066	0.02%
cwtl	3064	0.02%
fadd	2852	0.02%
fucomp	2678	0.02%
orb	2481	0.02%
fildll	2418	0.02%
andl	2379	0.02%
setb	2337	0.01%
andb	2263	0.01%
552 rows more...

I accidentally created an infinite loop

I needed to iterate through all values of 32-bit unsigned integer, so I wrote:

#include <stdint.h>
for (uint32_t i=0; i <= UINT32_MAX; i++) {
   // whatever
}

Is it ok? No, because value of uint32_t will never exceed UINT32_MAX = 0xffffffff. Of course we can use larger types, like uint64_t, but on 32-bit machines this requires some additional instructions. For example gcc 4.7 compiled following code:

void loop1(void (*fun)()) {
 for (uint64_t i=0; i <= UINT32_MAX; i++) {
  fun();
 }
}

to:

00000000 :
   0:   57                      push   %edi
   1:   bf 01 00 00 00          mov    $0x1,%edi
   6:   56                      push   %esi
   7:   31 f6                   xor    %esi,%esi
   9:   53                      push   %ebx
   a:   8b 5c 24 10             mov    0x10(%esp),%ebx
   e:   66 90                   xchg   %ax,%ax
  10:   ff d3                   call   *%ebx
  12:   83 c6 ff                add    $0xffffffff,%esi
  15:   83 d7 ff                adc    $0xffffffff,%edi
  18:   89 f8                   mov    %edi,%eax
  1a:   09 f0                   or     %esi,%eax
  1c:   75 f2                   jne    10 
  1e:   5b                      pop    %ebx
  1f:   5e                      pop    %esi
  20:   5f                      pop    %edi
  21:   c3                      ret    

TBH, I have no idea, why such weird sequence has been generated (add, adc, or, jnz). The simplest and portable solution is to detect wrap-around 32-bit value after increment:

uint32_t i=0;
while (1) {
 // loop body

 i += 1;
 if (i == 0) // wrap-around
  break;
}

In assembly code it's even simpler, because CPU sets the carry flag:

     xor %eax, %eax
loop:
     ; loop body

     add $1, %eax
     jnc loop

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.

The most complex instruction in x86 ISA

This instruction is... FBSTP. The instruction is advanced converter from integer part of floating point number to BCD format. Sounds weird? Read more.

itoa on steroids

Function itoa implemented with SSE instructions is about 3 times faster then scalar CPU version. Read more.

SSSE3 population count vs hardware

Peter Kankowski compared speed of SSE4.2 instructions crc32 and popcnt against software implementations. Hardware CRC32 is significantly faster, but population count is slightly slower than my SSSE3 popcount!

Speedup reversing table of bytes

With help of BSWAP instruction or SSE instructions (PSHUFD, PSHUFLW, PSHUFHW) or SSSE3 instruction (PSHUFB) reversing table can be faster. Speedup depends on three factors:

• table size: larger=faster
• table address: aligned=faster/much faster (15.5 speedup - possible! see chart)
• CPU type

Read full article

Determining if an integer is a power of 2

Method from Bit Twiddling Hacks: (x != 0) && (x & (x-1) == 0). GCC compiles this to following code:

; input/ouput: eax
; destroys: ebx

        test    %eax,  %eax     ; x == 0?
        jz      1f

        leal -1(%eax), %ebx     ; ebx := x-1
        test    %eax,  %ebx     ; ZF  := (eax & ebx == 0)

        setz     %al
        movzx    %al, %eax       ; eax := ZF
1:

We can use also BSF and BSR instructions, that determines position of first and last bit=1. If number is power of 2, then just one bit is set, and thus these positions are equal. BSx sets also ZF flag if input value is zero.

; input/output: eax
; destroys: ebx, edx

        bsf     %eax, %ebx      ; ebx := LSB's position if eax != 0, ZF = 1 if eax = 0
        jz      1f
        bsr     %eax, %edx      ; edx := MSB's position

        cmp     %ebx, %edx      ; ZF  := (ebx = edx)

        setz    %al
        movzx   %al, %eax       ; eax := ZF
1:

Brenchless conditional exchange

Suppose we have to exchange (or just move) two registers A and B:

1. C := A xor B
2. C := 0 if condition is not true
3. A := A xor C
4. B := B xor C

If C is 0, then  A and B left unchanged, else A and B are swapped. If only conditional move from B to A is needed, then step 4th have to be skipped.

Here is a sample x86 code, where condition is value of CF:

sbb edx, edx ; part of step 2. - edx = 0xffffff if CF=1, 0x000000 otherwise
mov ecx, eax
xor ecx, ebx ; step 1
and ecx, edx ; completed step 2. - now C is 0 or (A xor B)
xor eax, ecx ; step 3
xor ebx, ecx ; step 4

Branchless moves are possible in Pentium Pro and higher with instructions cmovcc.

See also XOR linked list.

Branchless signum

Problem: calculate value of sign(x):

• -1 when x < 0
• 0 when x = 0,
• +1 when x > 0.

My solution do not involve any hardware specific things like ALU flags nor special instructions - just plain AND, OR, shifts.

 ; input: eax = X

 movl %eax, %ebx
 sarl $31, %eax  // eax = -1 if X less then zero, 0 otherwise

 andl $0x7fffffff, %ebx
 addl $0x7fffffff, %ebx // MSB is set if any lower bits were set
 shrl $31, $ebx  // eax = +1 if X greater then zero, 0 otherwise

 orl %ebx, %eax  // eax = result

C99 implementation:

int32_t sign(int32_t x) {
 int32_t y;
 y = (x & 0x7fffffff) + 0x7fffffff;
 return (x >> 31) | ((uint32_t)y >> 31);
}

Fill word with selected bit

This is continuation of subproblem from previous post: we have a word (byte, dword, whatever) and want to fill it with selected bit.

1. The most general algorithm:

1. mask bit
   [10111010] => [00010000]
   
2. clone word
   a=[00010000], b=[00010000]
3. shift bit in first word to MSB, and to LSB in second word
   a=[10000000], b=[00000001]
4. subtract c = a - b
   c=[01111111]
5. add missing MSB c = c OR a
   c=[11111111]

2. If arithmetic shifts are supported:

1. shift bit to MSB
   a=[10000000]
2. arithmetic shift right
   a=[11111111]

3. On processor 386+ we can copy selected bit from reg to carry flag (CF) with bt reg, reg and then clone CF with sbb reg, reg.

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>

uint32_t fill1(uint32_t a, int bit) {
uint32_t b;

b = a = a & (1 << bit);

a <<= 31 - bit;
b >>= bit;

return (a - b) | a;
}

uint32_t fill2(uint32_t a, int bit) {
a <<= 31 - bit;
return (int32_t)(a) >> 31;
}

uint32_t fill386(uint32_t a, int bit) {
uint32_t result;             
__asm__ __volatile__ (
"bt  %1, %0\n"
"sbb %0, %0\n"
: "=r" (result)
: "r" (bit), "0" (a)
);
return result;
}

int main(int argc, char* argv[]) {
uint32_t x, i;
if (argc > 1)
x = (unsigned long)strtol(argv[1], NULL, 0);

printf("input = %08x\n", x);
for (i=0; i < 32; i++)
printf("bit %2d: fill1 = %08x, fill2 = %08x, fill386 = %08x\n",
i,
fill1(x, i),
fill2(x, i),
fill386(x, i)
);

return EXIT_SUCCESS;
}

Transpose bits in byte using SIMD instructions

Method presented here allows to get any bit permutation, transposition is just one of possible operations. Lookup-based approach would be faster, but algorithm is worth to (re)show.

Algorithm outline for 8-byte vector (with SSE instruction it is possible to get 2 operations in parallel):

1. fill vector with given byte
   [11010001] =>
   [11010001|11010001|11010001|11010001|11010001|11010001|11010001|11010001]
   
2. leave one bit per byte
   [10000000|01000000|00000000|00010000|00000000|00000000|00000000|00000001]
   
3. perform desired transposition ("move" bits around)
   [00000001|000000010|00000000|00001000|00000000|00000000|00000000|10000000]
   
4. perform horizontal OR of all bytes
   [10001011]
   

Here is my old MMX code (polish text); below SSE/SSE5 implementation details.

Ad 1. Series of punpcklbw/punpcklwb/shufps or pshufb if CPU supports SSSE3.

# 1.1
movd      %eax, %xmm0
shufps    $0x00, %xmm0, %xmm0
punpcklbw %xmm0, %xmm0
punpcklwd %xmm0, %xmm0

# 1.2
pxor      %xmm1, %xmm1
movd       %eax, %xmm0
pshufb    %xmm1, %xmm0

Ad 2. Simple pand with mask packed_qword(0x8040201008040201).

pand  MASK1, %xmm0

Ad 3. If plain SSE instructions are supported this step require some work. First each bit is populated to fill whole byte (using pcmpeq - we get negated result), then mask bits on desired positons.

SSE5 has powerful instruction protb that can do perform rotation of each byte with independent amount - so in this case just one instruction is needed.

# SSE
pcmpeq  %xmm1, %xmm0
pandn   MASK2, %xmm0

# SSE5
protb    ROT, %xmm0, %xmm0

Ad 4. Since bits are placed on distinct positions, we can use instruction psadbw, that calculate horizontal sums of bytewide differences from two registers (separately for low and high registers half). If one register is full of zeros, we get sum of bytes from other register.

psadbw  %xmm1, %xmm0
movd    %xmm0, %eax

Depending on instruction set three (SSE) or two (SSE5) additional tables are needed.