4 ms·
GCC's vectorized loop is using the usual "round-up" method, but we can easily implement the distributive method there by mixing 32 and 64-bit arithmetic: ;
by pbsd 12y ago
GCC's vectorized loop is using the usual "round-up" method, but we can easily implement the distributive method there by mixing 32 and 64-bit arithmetic:
; xmm7 = 0 magic 0 magic
; xmm6 = 0
; xmm5 = 0 (accumulator)
division_loop:
vmovdqu xmm4, [rax] ; xmm4 = w3 w2 w1 w0
vpunpckldq xmm0, xmm4, xmm6 ; xmm0 = 0 w1 | 0 w0
vpunpckhdq xmm2, xmm4, xmm6 ; xmm2 = 0 w3 | 0 w2
vpmuludq xmm1, xmm0, xmm7 ; xmm1 = w1 * magic | w0 * magic
vpmuludq xmm3, xmm2, xmm7 ; xmm3 = w3 * magic | w2 * magic
vpaddq xmm1, xmm1, xmm0 ; xmm1 = w1 * magic + w1 | w0 * magic + w0
vpaddq xmm3, xmm3, xmm2 ; xmm3 = w3 * magic + w3 | w2 * magic + w2
vpsrlq xmm1, xmm1, XX ; XX = shift constant
vpsrlq xmm3, xmm3, XX ; XX = shift constant
vshufps xmm1, xmm1, xmm3, 10001000b ; back to w3 w2 w1 w0
vpaddd xmm5, xmm5, xmm1 ; add to accumulator
add rax, 16
cmp rax, rbp
jnz division_loop
The shifts can be further tweaked to allow to use a blend instead of a shuffle, in case there's a penalty for changing from integer to floating-point domains. Saturating increment also works well, despite there not being a PADDUSD instruction; it can be done in two cycles:
; xmm6 = -1 -1 -1 -1
vpcmpeqd xmm3, xmm4, xmm6 ; xmm3[i] = (xmm4[i] == -1 ? -1 : 0)
vpsubd xmm4, xmm4, xmm6 ; xmm4[i] += 1
vpaddd xmm4, xmm4, xmm3 ; xmm4[i] += xmm3[i]
- nkurz 12y agoNice work! Seems like it extends fine to AVX2 for doing 8 divides at a time in less than 1 cycle per division if for some reason you need to divide a whole array by the same constant. Though rare, this is impressive considering the tradition of expensive division. I don't know if you will have a domain penalty for the vshufps, but I think the answer is yes, and that the blend would be better. The blend can also run on p015, versus only p5 for the shufps, so it might pipeline better. I hadn't realized until now that vshufps had different semantics than vpshufd and is able to combine from two different registers. If there is no penalty, this could be quite a useful instruction. Any reason for using 'w' for your 32-bit elements rather than 'd' for double-word as used in the instruction names?
- pbsd 12y agoAccording to Agner Fog, data domain penalties with SHUFPS are nonexistent, or at least pretty hard to trigger, on Haswell. Previous chips might complain, though. Problem with blends is that they require SSE4.1+, unlike the current version that is fully SSE2 (modulo AVX encoding niceties). I chose 'w' for no real reason other than being thinking of 'word' at the time.
- nkurz 12y agoI just tested, and I can now confirm that there is no domain switching penalty for using vshufps versus vpshufd on integer data for XMM and YMM on Haswell, nor for XMM on Sandy Bridge (no YMM comparison possible there because vpshufd requires AVX2).