我正在尝试为以下函数编写AXV版本：
void
hashids_shuffle(char *str, size_t str_length, char *salt, size_t salt_length)
{
size_t i, j, v, p;
char temp;
if (!salt_length) {
return;
}
for (i = str_length - 1, v = 0, p = 0; i > 0; --i, ++v) {
v %= salt_length;
p += salt[v];
j = (salt[v] + v + p) % i;
temp = str[i];
str[i] = str[j];
str[j] = temp;
}
}
我正在尝试对v %= salt_length;进行矢量化
我想初始化一个包含从0到str_length的数字的向量，以便使用SVML的_mm_rem_epu64来计算每个循环迭代的v
如何正确初始化矢量？
仅仅询问如何初始化向量基本上就是询问教程。在谷歌上搜索英特尔的一些内部函数使用指南。我将假装这个问题不是微不足道的，回答如何有效地实现这个函数。这绝对不是一个初学者想要矢量化的函数。
有关文档的链接，请参阅x86标记wiki，尤其是英特尔的内部指南。请参阅sse标签wiki，以获取一个链接，该链接指向使用sse内部函数进行SIMD编程的非常好的介绍，以及如何有效地使用SIMD，以及其他链接。
内容摘要：
- 展开/矢量化免费处理
v % salt_length - 如果
v++; v %= loop_invariant;不是2的幂或编译时间常数，如何对其进行矢量化。包括关于使用_mm_set_epi8或为此目的初始化矢量的其他方法的标题问题的答案 - 如何开始对这样的复杂循环进行矢量化：从展开开始查找串行依赖项
-
完整函数的未测试版本，它对除
% i和交换之外的所有内容进行矢量化。(也就是说，按照您的要求，对所有廉价的操作进行矢量化)。-
(v + salt[v] + p)(及其之前的所有内容)向量化为两条vpaddw指令。循环外用于矢量化p的前缀和设置很棘手，但我最终也对其进行了矢量化。 -
函数的绝大多数运行时间将在
j元素向量上的标量内部循环中，div上的瓶颈(或SVML可以做的任何事情)，和/或具有非常大字符串的缓存未命中。
-
由于与伪随机索引的交换会产生不可预测的串行依赖关系，因此整个循环无法轻松向量化。使用AVX512聚集+混洗+分散，使用AVX512CD查找冲突位掩码可能是可能的，但这必须是一个单独的问题。我不确定高效地做到这一点有多难，或者你是否经常重复多次向量洗牌，只在一个不冲突的元素上取得进展。
由于salt_length = sizeof(size_t)是一个编译时间常数，并且比向量长度小2的幂，因此v++和v%=salt_length根本不需要循环内的任何代码，并且作为有效展开循环以并行执行多个v值的副作用，它是免费的。
(使用依赖于平台的salt大小意味着32位构建将无法处理使用64位salt创建的数据。即使是x32 ABI也有32位size_t，因此更改为uint64_t似乎是有意义的，除非您永远不需要在机器之间共享salt哈希。)
在标量循环中，v遵循0 1 2 3 0 1 2…的重复模式。。。(或0..7，对于64位)。在矢量代码中，我们可能对32字节矢量中的4B元素同时执行8个v值，或者对2B元素同时执行16次迭代。
因此v成为循环不变的常向量。有趣的是，salt[v]也是如此，所以我们永远不需要在循环中进行任何盐表查找实际上，v+salt[v]可以预先计算标量和矢量。
标量版本应该预先计算v+salt[v]，并也按4或8展开，删除LUT查找，以便所有内存/缓存吞吐量都可用于实际交换。编译器可能不会为您这样做，所以您可能需要手动展开并编写额外的代码来处理最后奇数个字符串字节。(在不展开的情况下，您仍然可以预先计算v+salt[v]的查找表，该查找表的类型足够宽，无法环绕)。
即使只是确保salt_length在编译时是已知的，也可以获得更好的代码。v %= compile_time_constant比divinsn便宜，当它是2的幂时非常便宜。(它刚好变成v &= 7)。如果标量版本可以内联，或者如果您使用salt_length = sizeof(size_t)而不是将其作为函数arg传递，编译器可能会为您这样做。
如果您还不知道salt_length:，即在您透露有关salt_length的关键信息之前@harold的建议：
由于我们一开始就知道v < salt_length，所以我们最多只需要一个v -= salt_length就可以将其包装回正确的范围并保持不变。这被称为"强度减少"运算，因为减法比除法更弱(也更便宜)。
// The scalar loop would benefit from this transformation, too.
// or better, unroll the scalar loop by 8 so everything becomes constant
v++;
if( v >= salt_length)
v-= salt_length;
要对进行矢量化这个：假设我们只知道salt_length <= 16，所以我们可以使用32个uint8_t值的矢量。(我们可以使用pshufb对salt[v]LUT查找进行矢量化)。
// untested // Vectorizing v++; v %= unknown_loop_invariant_value;
if (!salt_length) return;
assert(salt_length <= 16); // so we can use pshufb for the salt[v] step
__m256i vvec = _mm256_setr_epi8( // setr: lowest element first, unlike set
0%salt_length, 1%salt_length, 2%salt_length, 3%salt_length,
4%salt_length, 5%salt_length, 6%salt_length, 7%salt_length,
8%salt_length, 9%salt_length, 10%salt_length, 11%salt_length,
12%salt_length, 13%salt_length, 14%salt_length, 15%salt_length,
16%salt_length, 17%salt_length, 18%salt_length, 19%salt_length,
20%salt_length, 21%salt_length, 22%salt_length, 23%salt_length,
24%salt_length, 25%salt_length, 26%salt_length, 27%salt_length,
28%salt_length, 29%salt_length, 30%salt_length, 31%salt_length);
__m256i v_increment = _mm256_set1_epi8(32 % salt_length);
__m256i vmodulus = _mm256_set1_epi8(salt_length);
// salt_lut = _mm256_set1_epi64x(salt_byval); // for known salt length. (pass it by val in a size_t arg, instead of by char*).
// duplicate the salt into both lanes of a vector. Garbage beyond salt_length isn't looked at.
__m256i salt_lut = _mm256_broadcastsi128_si256(_mm_loadu_si128(salt)); // nevermind that this could segfault if salt is short and at the end of a page.
//__m256i v_plus_salt_lut = _mm256_add_epi8(vvec, salt_lut); // not safe with 8-bit elements: could wrap
// We could use 16-bit elements and AVX512 vpermw (or vpermi2w to support longer salts)
for (...) {
vvec = _mm256_add_epi8(vvec, v_increment); // ++v;
// if(!(salt_length > v)) { v-= salt_length; }
__m256i notlessequal = _mm256_cmpgt_epi8(vmodulus, vvec); // all-ones where salt_length > v.
// all-zero where salt_length <= v, where we need to subtract salt_length
__m256i conditional_sub = _mm256_and_si256(notlessequal, vmodulus)
vvec = _mm256_sub_epi8(vvec, conditional_sub); // subtract 0 or salt_length
// salt[v] lookup:
__m256i saltv = _mm256_shuffle_epi8(salt_lut, vvec); // salt[v]
// then maybe pmovzx and vextracti128+pmovzx to zero-extend to 16-bit elements? Maybe vvec should only be a 16-bit vector?
// or unpack lo/hi with zeros (but that behaves differently from pmovzx at the lane boundary)
// or have vvec already holding 16-bit elements with the upper half of each one always zero. mask after the pshufb to re-zero,
// or do something clever with `vvec`, `v_increment` and `vmodulus` so `vvec` can have `0xff` in the odd bytes, so pshufb zeros those elements.
}
当然，如果我们知道salt_length是2的幂，我们应该屏蔽掉每个元素中除相关低位之外的所有低位：
vvec = _mm256_add_epi8(vvec, _mm256_set1_epi8(salt_length)); // ++v;
vvec = _mm256_and_si256(vvec, _mm256_set1_epi8(salt_length - 1)); // v &= salt_length - 1; // aka v%=salt_length;
注意到我们从错误的元素大小开始，是因为我们意识到一次只矢量化一行是个坏主意，因为现在我们必须更改所有已经编写的代码，以使用更宽的元素，有时需要不同的策略来做同样的事情。
当然，你确实需要从一个粗略的大纲开始，无论是心理上的还是书面上的，你可以分别做每一步。正是在仔细思考的过程中，你才能看到不同的部分是如何结合在一起的。
对于复杂循环，有用的第一步可能是尝试手动展开标量循环。这将有助于找到串行依赖关系，以及通过展开来简化的事情
(stuff) % i：最难的部分
我们需要足够宽的元素来保持i的最大值，因为i是而不是2的幂，并且不是常数，所以模运算需要功。任何更宽的范围都是浪费，会降低我们的吞吐量。如果我们可以对循环的其余部分进行矢量化，那么对于不同范围的str_length，使用不同版本专门化函数可能是值得的。(或者可能使用64b元素循环，直到i<=UINT32_MAX，然后循环，直到i<=UINT16_MAX，等等)。如果您知道不需要处理大于4GiB的字符串，那么只需使用32位数学运算就可以加快常见情况的速度。(64位除法比32位除法慢，即使高位都为零)。
实际上，我认为我们需要与最大p一样宽的元素，因为它会永远累积(直到它在原始标量代码中以2^64结束)。与常数模不同，我们不能只使用p%=i来控制它，即使模是分布的。CCD_ 38。即使对齐到16也无济于事：12345%32！=12345%48%32
这将很快使p过大，无法重复条件减法i(直到条件掩码全部为假)，即使对于相当大的i值也是如此。
通过已知的整数常数进行模运算有一些技巧(请参阅http://libdivide.com/)，但AFAIK为附近的除数(即使是像16这样的两步幂)计算乘法模逆并不比为完全独立的数容易。因此，我们不能只为i值的下一个向量调整常数。
小数定律也许让我们有必要剥离最后一对向量迭代，具有乘法模逆的预先计算的向量所以CCD_ 43可以用矢量来完成。一旦我们接近字符串的末尾，它可能在一级缓存中很热，所以我们在div上完全陷入瓶颈，而不是交换加载/存储。为此，我们可能会使用序言来达到16的倍数的i值，因此当我们接近i=0时，最后几个向量总是具有相同的i值排列。否则，我们将有一个i值范围内的常数LUT，并简单地从中进行未对齐的加载。这意味着我们不必旋转salt_v和p。
转换为FP可能会很有用，因为最近的Intel CPU(尤其是Skylake)具有非常强大的FP分割硬件，具有显著的流水线(吞吐量：延迟比)。如果我们能通过正确的四舍五入选择得到准确的结果，那就太好了。(float和double可以精确地表示高达其尾数大小的任何整数。)
我想值得一试英特尔的_mm_rem_epu16(用i值的向量，用set1(16)向量递减)。如果他们使用FP来获得准确的结果，那就太好了。如果它只是对标量进行解包并进行整数除法，那么将浪费时间将值返回到向量中。
无论如何，最简单的解决方案当然是用标量循环迭代向量元素。在您使用AVX512CD进行交换之前，这似乎是合理的，但如果交换在一级缓存中都很热，那么它可能比交换慢一个数量级。
(未测试)函数的部分矢量化版本：
以下是Godbolt编译器资源管理器上的代码，包含完整的设计注释注释，包括我在计算SIMD前缀sum算法时绘制的图表。我最终记得，在@ZBoson的浮点SSE前缀和答案中，我看到了一个更窄的版本，作为一个构建块，但直到我自己重新发明了它之后。
// See the godbolt link for full design-notes comments
// comments about what makes nice ASM or not.
#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h>
#include <assert.h>
static inline
__m256i init_p_and_increment(size_t salt_length, __m256i *p_increment, __m256i saltv_u16, __m128i saltv_u8)
{ // return initial p vector (for first 16 i values).
// return increment vector by reference.
if (salt_length == 4) {
assert(0); // unimplemented
// should be about the same as length == 8. Can maybe factor out some common parts, like up to psum2
} else {
assert(salt_length == 8);
// SIMD prefix sum for n elements in a vector in O(log2(n)) steps.
__m128i sv = _mm256_castsi256_si128(saltv_u16);
__m128i pshift1 = _mm_bslli_si128(sv, 2); // 1 elem (uint16_t)
__m128i psum1 = _mm_add_epi16(pshift1, sv);
__m128i pshift2 = _mm_bslli_si128(psum1, 4); // 2 elem
__m128i psum2 = _mm_add_epi16(pshift2, psum1);
__m128i pshift3 = _mm_bslli_si128(psum2, 8); // 4 elem
__m128i psum3 = _mm_add_epi16(pshift3, psum2); // p_initial low 128. 2^3 = 8 elements = salt_length
// psum3 = the repeating pattern of p values. Later values just add sum(salt[0..7]) to every element
__m128i p_init_low = psum3;
__m128i sum8_low = _mm_sad_epu8(saltv_u8, _mm_setzero_si128()); // sum(s0..s7) in each 64-bit half
// alternative:
// sum8_low = _mm_bsrli_si128(p_init_low, 14); // has to wait for psum3 to be ready: lower ILP than doing psadbw separately
__m256i sum8 = _mm256_broadcastw_epi16(sum8_low);
*p_increment = _mm256_slli_epi16(sum8, 1); // set1_epi16(2*sum(salt[0..7]))
__m128i p_init_high = _mm_add_epi16(p_init_low, _mm256_castsi256_si128(sum8));
__m256i p_init = _mm256_castsi128_si256(p_init_low);
p_init = _mm256_inserti128_si256(p_init, p_init_high, 1);
// not supported by gcc _mm256_set_m128i(p_init_high, psum3);
return p_init;
}
}
void
hashids_shuffle_simd(char *restrict str, size_t str_length, size_t salt_byval)
{
//assert(salt_length <= 16); // so we can use pshufb for the salt[v] step for non-constant salt length.
// platform-dependent salt size seems weird. Why not uint64_t?
size_t salt_length = sizeof(size_t);
assert(str_length-1 < UINT16_MAX); // we do p + v + salt[v] in 16-bit elements
// TODO: assert((str_length-1)/salt_length * p_increment < UINT16_MAX);
__m128i saltv_u8;
__m256i v, saltv;
if(salt_length == 4) {
v = _mm256_set1_epi64x(0x0003000200010000); // `v%salt_length` is 0 1 2 3 repeating
saltv_u8 = _mm_set1_epi32( salt_byval );
saltv = _mm256_cvtepu8_epi16( saltv_u8 ); // salt[v] repeats with the same pattern: expand it to 16b elements with pmovzx
} else {
assert(salt_length == 8);
v = _mm256_cvtepu8_epi16( _mm_set1_epi64x(0x0706050403020100) );
saltv_u8 = _mm_set1_epi64x( salt_byval );
saltv = _mm256_cvtepu8_epi16( saltv_u8 );
}
__m256i v_saltv = _mm256_add_epi16(v, saltv);
__m256i p_increment;
__m256i p = init_p_and_increment(salt_length, &p_increment, saltv, saltv_u8);
for (unsigned i=str_length-1; i>0 ; /*i-=16 */){
// 16 uint16_t j values per iteration. i-- happens inside the scalar shuffle loop.
p = _mm256_add_epi16(p, p_increment); // p += salt[v]; with serial dependencies accounted for, prefix-sum style
__m256i j_unmodded = _mm256_add_epi16(v_saltv, p);
// size_t j = (v + saltv[v] + p) % i;
//////// scalar loop over 16 j elements, doing the modulo and swap
// alignas(32) uint16_t jbuf[16]; // portable C++11 syntax
uint16_t jbuf[16] __attribute__((aligned(32))); // GNU C syntax
_mm256_store_si256((__m256i*)jbuf, j_unmodded);
const int jcount = sizeof(jbuf)/sizeof(jbuf[0]);
for (int elem = 0 ; elem < jcount ; elem++) {
if (--i == 0) break; // in fact returns from the whole function.
// 32-bit division is significantly faster than 64-bit division
unsigned j = jbuf[elem] % (uint32_t)i;
// doubtful that vectorizing this with Intel SVML _mm_rem_epu16 would be a win
// since there's no hardware support for it. Until AVX512CD, we need each element in a gp reg as an array index anyway.
char temp = str[i];
str[i] = str[j];
str[j] = temp;
}
}
}
这编译成了看起来不错的asm，但我还没有运行它
Clang做了一个相当合理的内部循环。这是-fno-unroll-loops的可读性。至于性能，请忽略这一点，尽管这并不重要，因为循环开销不是瓶颈。
# The loop part of clang3.8.1's output. -march=haswell -fno-unroll-loops (only for human readability. Normally it unrolls by 2).
.LBB0_6: # outer loop # in Loop: Header=BB0_3 Depth=1
add esi, 1
.LBB0_3: # first iteration entry point # =>This Loop Header: Depth=1
vpaddw ymm2, ymm2, ymm1 # p += p_increment
vpaddw ymm3, ymm0, ymm2 # v+salt[v] + p
vmovdqa ymmword ptr [rsp], ymm3 # store jbuf
add esi, -1
lea r8, [rdi + rsi]
mov ecx, 1
.LBB0_4: # inner loop # Parent Loop BB0_3 Depth=1
# gcc's version fully unrolls the inner loop, leading to code bloat
test esi, esi # if(i==0) return
je .LBB0_8
movzx eax, word ptr [rsp + 2*rcx - 2] # load jbuf
xor edx, edx
div esi
mov r9b, byte ptr [r8] # swap
mov al, byte ptr [rdi + rdx]
mov byte ptr [r8], al
mov byte ptr [rdi + rdx], r9b
add esi, -1
add r8, -1
cmp rcx, 16 # silly clang, not macro-fusing cmp/jl because it wants to use a weird way to increment.
lea rcx, [rcx + 1]
jl .LBB0_4 # inner loop
jmp .LBB0_6 # outer loop
[免责声明：考虑32位应用程序-其中size_t为unsigned int]
对齐分配可以通过使用Aligned_malloc函数来完成，您可以在windows和linux中找到这些函数。
以这种方式分配字符串(到64字节边界)将允许您通过对所有字节使用对齐的loads _mm256-load_si256将数据直接加载到_mm256i寄存器中。
如果字符串未正确对齐，则可以在开头使用未对齐的load_mm256_load_si256加载第一个字节。
您执行的第一个模运算(v%=salt_length)是用一个常量操作数完成的，您可以在avx寄存器中的循环之前使用_mm256_set1_epi32:对其进行初始化
__m256i mod = _mm256_set2_epi32(salt_length);
对于下一个，您可以使用_mm256_set_epi32，它使用提供的所有值初始化寄存器(注意的逆序)。
请注意，如果您使用的是AVX2或AVX512(而不仅仅是AVX——您的问题对指令集有点困惑)，您也可以使用collect指令加载数据，该指令从第二个参数中指定的索引处的向量加载数据。
在余数运算中可以与AVX512一起使用的最小数字是8位，内部为：_mm512_rem_epu8
然而，要用0到N的值初始化它的输入，您需要使用_mm512_set_epi32并将8位整数打包为32位整数，因为似乎没有一个内在的64个整数，每个8位。代码看起来像：
const __m512i sseConst = _mm512_set_epi32(
(63<<24) | (62<<16) | (61<<8) | 60,
(59<<24) | (58<<16) | (57<<8) | 56,
... etc ...);
如果你不喜欢这样的打字或害怕打字错误，可以考虑为此编写一个代码生成器。
如果您不使用malloc()分配__m512i类型，那么它的使用应该自动为您处理对齐。否则，在"aligned malloc"中搜索编译器。您需要的对齐方式是64字节(等于512位)。
当你需要向量中的下一个N整数时，你可以做：
const __m512i inc = _mm512_set1_epi32((N<<24) | (N<<16) | (N<<8) | N);
然后可以添加inc和sseConst(内在_mm512_add_epi32)。