如何使用 SSE 将 16 位整数除以 255？
How to divide 16-bit integer by 255 with using SSE?
我处理图像处理。
我需要将 16 位整数 SSE 向量除以 255.
我不能使用像_mm_srli_epi16()这样的移位运算符，因为255不是2的倍数。
我当然知道可以将整数转换为浮点数，执行除法然后再转换回整数。
但也许有人知道另一种解决方案...
有一个除以255的整数近似值：
inline int DivideBy255(int value)
{
return (value + 1 + (value >> 8)) >> 8;
}
所以使用 SSE2 时它看起来像：
inline __m128i DivideI16By255(__m128i value)
{
return _mm_srli_epi16(_mm_add_epi16(
_mm_add_epi16(value, _mm_set1_epi16(1)), _mm_srli_epi16(value, 8)), 8);
}
对于 AVX2：
inline __m256i DivideI16By255(__m256i value)
{
return _mm256_srli_epi16(_mm256_add_epi16(
_mm256_add_epi16(value, _mm256_set1_epi16(1)), _mm256_srli_epi16(value, 8)), 8);
}
对于 Altivec（电源）：
typedef __vector int16_t v128_s16;
const v128_s16 K16_0001 = {1, 1, 1, 1, 1, 1, 1, 1};
const v128_s16 K16_0008 = {8, 8, 8, 8, 8, 8, 8, 8};
inline v128_s16 DivideBy255(v128_s16 value)
{
return vec_sr(vec_add(vec_add(value, K16_0001), vec_sr(value, K16_0008)), K16_0008);
}
对于霓虹灯（ARM）：
inline int16x8_t DivideI16By255(int16x8_t value)
{
return vshrq_n_s16(vaddq_s16(
vaddq_s16(value, vdupq_n_s16(1)), vshrq_n_s16(value, 8)), 8);
}
GCC 使用 x 优化 x/255 是 unsigned short 到 DWORD(x * 0x8081) >> 0x17 可以进一步简化为 HWORD(x * 0x8081) >> 7 最后 HWORD((x << 15) + (x << 7) + x) >> 7.
SIMD 宏可以如下所示：
#define MMX_DIV255_U16(x) _mm_srli_pi16(_mm_mulhi_pu16(x, _mm_set1_pi16((short)0x8081)), 7)
#define SSE2_DIV255_U16(x) _mm_srli_epi16(_mm_mulhi_epu16(x, _mm_set1_epi16((short)0x8081)), 7)
#define AVX2_DIV255_U16(x) _mm256_srli_epi16(_mm256_mulhi_epu16(x, _mm256_set1_epi16((short)0x8081)), 7)
如果您希望在所有情况下都得到完全正确的结果，请遵循 Marc Glisse's comment on the question Anton linked: SSE integer division?
的建议
使用 GNU C 原生向量语法来表示向量除以给定标量，and see what it does on the Godbolt compiler explorer:
无符号除法很便宜：
typedef unsigned short vec_u16 __attribute__((vector_size(16)));
vec_u16 divu255(vec_u16 x){ return x/255; } // unsigned division
#gcc5.5 -O3 -march=haswell
divu255:
vpmulhuw xmm0, xmm0, XMMWORD PTR .LC3[rip] # _mm_set1_epi16(0x8081)
vpsrlw xmm0, xmm0, 7
ret
内部版本：
// UNSIGNED division with intrinsics
__m128i div255_epu16(__m128i x) {
__m128i mulhi = _mm_mulhi_epu16(x, _mm_set1_epi16(0x8081));
return _mm_srli_epi16(mulhi, 7);
}
如果您在前端吞吐量或 Intel CPU 上的端口 0 吞吐量方面遇到瓶颈，这仅需 2 微指令，吞吐量比@ermlg 的答案更好（但延迟更差）。 （与往常一样，当您将其用作更大函数的一部分时，它取决于周围的代码。）http://agner.org/optimize/
矢量移位仅 运行 在 Intel 芯片上的端口 0 上，因此 @ermlg 的 2 次移位 + 1 在端口 0 上增加了瓶颈。（同样取决于周围的代码）。这是 3 微指令对 2 微指令。
在 Skylake 上，pmulhuw / pmulhw 运行 在端口 0 或 1 上，因此它可以 运行 与移位并行。 （但在 Broadwell 和更早版本上，它们 运行 仅在端口 0 上，与轮班冲突。因此，Intel pre-Skylake 的唯一优势是前端的总微指令数更少，并且乱序执行保持track of.) pmulhuw 在 Intel 上有 5 个周期延迟，而在轮班时有 1 个周期延迟，但是当您可以节省 uops 以获得更多吞吐量时，OoO exec 通常可以隐藏几个周期更多的延迟。
Ryzen 在其 P0 上也只有 运行s pmulhuw，但在 P2 上移动，所以它非常适合这个。
但是有符号整数除法舍入语义与移位不匹配
typedef short vec_s16 __attribute__((vector_size(16)));
vec_s16 div255(vec_s16 x){ return x/255; } // signed division
; function arg x starts in xmm0
vpmulhw xmm1, xmm0, XMMWORD PTR .LC3[rip] ; a vector of set1(0x8081)
vpaddw xmm1, xmm1, xmm0
vpsraw xmm0, xmm0, 15 ; 0 or -1 according to the sign bit of x
vpsraw xmm1, xmm1, 7 ; shift the mulhi-and-add result
vpsubw xmm0, xmm1, xmm0 ; result += (x<0)
.LC3:
.value -32639
.value -32639
; repeated
冒着膨胀答案的风险，这里又是内在函数：
// SIGNED division
__m128i div255_epi16(__m128i x) {
__m128i tmp = _mm_mulhi_epi16(x, _mm_set1_epi16(0x8081));
tmp = _mm_add_epi16(tmp, x); // There's no integer FMA that's usable here
x = _mm_srai_epi16(x, 15); // broadcast the sign bit
tmp = _mm_srai_epi16(tmp, 7);
return _mm_sub_epi16(tmp, x);
}
在 godbolt 输出中，请注意 gcc 足够聪明，可以在内存中为 set1 和它自己为 div255 生成的常量使用相同的 16B 常量。 AFAIK，这就像字符串常量合并一样工作。
出于好奇（如果性能是个问题），这里是使用的准确性
(val + offset) >> 8 作为 (val / 255) 的替代品，用于所有 16 位值，最高为 255*255（例如，当使用 8 位混合因子混合两个 8 位值时）：
(avrg signed error / avrg abs error / max abs error)
offset 0: 0.49805 / 0.49805 / 1 (just shifting, no offset)
offset 0x7F: 0.00197 / 0.24806 / 1
offest 0x80: -0.00194 / 0.24806 / 1
所有其他偏移量都会产生更大的有符号错误和平均错误。因此，如果您可以忍受小于 0.25 的平均误差，而不是使用偏移+移位来实现小幅速度提升
// approximate division by 255 for packs of 8 times 16bit values in vals_packed
__m128i offset = _mm_set1_epi16(0x80); // constant
__m128i vals_packed_offest = _mm_add_epi16( vals_packed, offset );
__m128i result_packed = _mm_srli_epi16( vals_packed_offest , 8 );
准确版本：
#define div_255_fast(x) (((x) + (((x) + 257) >> 8)) >> 8)
当x在[0, 65536]范围内时，ERROR为0，比x/255快一倍：
http://quick-bench.com/t3Y2-b4isYIwnKwMaPQi3n9dmtQ
SIMD 版本：
// (x + ((x + 257) >> 8)) >> 8
static inline __m128i _mm_fast_div_255_epu16(__m128i x) {
return _mm_srli_epi16(_mm_adds_epu16(x,
_mm_srli_epi16(_mm_adds_epu16(x, _mm_set1_epi16(0x0101)), 8)), 8);
}
对于高于 65535 的正 x，这是另一个版本：
static inline int32_t fast_div_255_any (int32_t n) {
uint64_t M = (((uint64_t)1) << 40) / 255 + 1; // "1/255" in 24.40 fixed point number
return (M * n) >> 40; // fixed point multiply: n * (1/255)
}
更广泛（需要 64 位 mul），但仍然比 div 指令快。
我处理图像处理。 我需要将 16 位整数 SSE 向量除以 255.
我不能使用像_mm_srli_epi16()这样的移位运算符，因为255不是2的倍数。
我当然知道可以将整数转换为浮点数，执行除法然后再转换回整数。
但也许有人知道另一种解决方案...
有一个除以255的整数近似值：
inline int DivideBy255(int value)
{
return (value + 1 + (value >> 8)) >> 8;
}
所以使用 SSE2 时它看起来像：
inline __m128i DivideI16By255(__m128i value)
{
return _mm_srli_epi16(_mm_add_epi16(
_mm_add_epi16(value, _mm_set1_epi16(1)), _mm_srli_epi16(value, 8)), 8);
}
对于 AVX2：
inline __m256i DivideI16By255(__m256i value)
{
return _mm256_srli_epi16(_mm256_add_epi16(
_mm256_add_epi16(value, _mm256_set1_epi16(1)), _mm256_srli_epi16(value, 8)), 8);
}
对于 Altivec（电源）：
typedef __vector int16_t v128_s16;
const v128_s16 K16_0001 = {1, 1, 1, 1, 1, 1, 1, 1};
const v128_s16 K16_0008 = {8, 8, 8, 8, 8, 8, 8, 8};
inline v128_s16 DivideBy255(v128_s16 value)
{
return vec_sr(vec_add(vec_add(value, K16_0001), vec_sr(value, K16_0008)), K16_0008);
}
对于霓虹灯（ARM）：
inline int16x8_t DivideI16By255(int16x8_t value)
{
return vshrq_n_s16(vaddq_s16(
vaddq_s16(value, vdupq_n_s16(1)), vshrq_n_s16(value, 8)), 8);
}
GCC 使用 x 优化 x/255 是 unsigned short 到 DWORD(x * 0x8081) >> 0x17 可以进一步简化为 HWORD(x * 0x8081) >> 7 最后 HWORD((x << 15) + (x << 7) + x) >> 7.
SIMD 宏可以如下所示：
#define MMX_DIV255_U16(x) _mm_srli_pi16(_mm_mulhi_pu16(x, _mm_set1_pi16((short)0x8081)), 7)
#define SSE2_DIV255_U16(x) _mm_srli_epi16(_mm_mulhi_epu16(x, _mm_set1_epi16((short)0x8081)), 7)
#define AVX2_DIV255_U16(x) _mm256_srli_epi16(_mm256_mulhi_epu16(x, _mm256_set1_epi16((short)0x8081)), 7)
如果您希望在所有情况下都得到完全正确的结果，请遵循 Marc Glisse's comment on the question Anton linked: SSE integer division?
的建议使用 GNU C 原生向量语法来表示向量除以给定标量，and see what it does on the Godbolt compiler explorer:
无符号除法很便宜：
typedef unsigned short vec_u16 __attribute__((vector_size(16)));
vec_u16 divu255(vec_u16 x){ return x/255; } // unsigned division
#gcc5.5 -O3 -march=haswell
divu255:
vpmulhuw xmm0, xmm0, XMMWORD PTR .LC3[rip] # _mm_set1_epi16(0x8081)
vpsrlw xmm0, xmm0, 7
ret
内部版本：
// UNSIGNED division with intrinsics
__m128i div255_epu16(__m128i x) {
__m128i mulhi = _mm_mulhi_epu16(x, _mm_set1_epi16(0x8081));
return _mm_srli_epi16(mulhi, 7);
}
如果您在前端吞吐量或 Intel CPU 上的端口 0 吞吐量方面遇到瓶颈，这仅需 2 微指令，吞吐量比@ermlg 的答案更好（但延迟更差）。 （与往常一样，当您将其用作更大函数的一部分时，它取决于周围的代码。）http://agner.org/optimize/
矢量移位仅 运行 在 Intel 芯片上的端口 0 上，因此 @ermlg 的 2 次移位 + 1 在端口 0 上增加了瓶颈。（同样取决于周围的代码）。这是 3 微指令对 2 微指令。
在 Skylake 上，pmulhuw / pmulhw 运行 在端口 0 或 1 上，因此它可以 运行 与移位并行。 （但在 Broadwell 和更早版本上，它们 运行 仅在端口 0 上，与轮班冲突。因此，Intel pre-Skylake 的唯一优势是前端的总微指令数更少，并且乱序执行保持track of.) pmulhuw 在 Intel 上有 5 个周期延迟，而在轮班时有 1 个周期延迟，但是当您可以节省 uops 以获得更多吞吐量时，OoO exec 通常可以隐藏几个周期更多的延迟。
Ryzen 在其 P0 上也只有 运行s pmulhuw，但在 P2 上移动，所以它非常适合这个。
但是有符号整数除法舍入语义与移位不匹配
typedef short vec_s16 __attribute__((vector_size(16)));
vec_s16 div255(vec_s16 x){ return x/255; } // signed division
; function arg x starts in xmm0
vpmulhw xmm1, xmm0, XMMWORD PTR .LC3[rip] ; a vector of set1(0x8081)
vpaddw xmm1, xmm1, xmm0
vpsraw xmm0, xmm0, 15 ; 0 or -1 according to the sign bit of x
vpsraw xmm1, xmm1, 7 ; shift the mulhi-and-add result
vpsubw xmm0, xmm1, xmm0 ; result += (x<0)
.LC3:
.value -32639
.value -32639
; repeated
冒着膨胀答案的风险，这里又是内在函数：
// SIGNED division
__m128i div255_epi16(__m128i x) {
__m128i tmp = _mm_mulhi_epi16(x, _mm_set1_epi16(0x8081));
tmp = _mm_add_epi16(tmp, x); // There's no integer FMA that's usable here
x = _mm_srai_epi16(x, 15); // broadcast the sign bit
tmp = _mm_srai_epi16(tmp, 7);
return _mm_sub_epi16(tmp, x);
}
在 godbolt 输出中，请注意 gcc 足够聪明，可以在内存中为 set1 和它自己为 div255 生成的常量使用相同的 16B 常量。 AFAIK，这就像字符串常量合并一样工作。
出于好奇（如果性能是个问题），这里是使用的准确性 (val + offset) >> 8 作为 (val / 255) 的替代品，用于所有 16 位值，最高为 255*255（例如，当使用 8 位混合因子混合两个 8 位值时）：
(avrg signed error / avrg abs error / max abs error)
offset 0: 0.49805 / 0.49805 / 1 (just shifting, no offset)
offset 0x7F: 0.00197 / 0.24806 / 1
offest 0x80: -0.00194 / 0.24806 / 1
所有其他偏移量都会产生更大的有符号错误和平均错误。因此，如果您可以忍受小于 0.25 的平均误差，而不是使用偏移+移位来实现小幅速度提升
// approximate division by 255 for packs of 8 times 16bit values in vals_packed
__m128i offset = _mm_set1_epi16(0x80); // constant
__m128i vals_packed_offest = _mm_add_epi16( vals_packed, offset );
__m128i result_packed = _mm_srli_epi16( vals_packed_offest , 8 );
准确版本：
#define div_255_fast(x) (((x) + (((x) + 257) >> 8)) >> 8)
当x在[0, 65536]范围内时，ERROR为0，比x/255快一倍：
http://quick-bench.com/t3Y2-b4isYIwnKwMaPQi3n9dmtQ
SIMD 版本：
// (x + ((x + 257) >> 8)) >> 8
static inline __m128i _mm_fast_div_255_epu16(__m128i x) {
return _mm_srli_epi16(_mm_adds_epu16(x,
_mm_srli_epi16(_mm_adds_epu16(x, _mm_set1_epi16(0x0101)), 8)), 8);
}
对于高于 65535 的正 x，这是另一个版本：
static inline int32_t fast_div_255_any (int32_t n) {
uint64_t M = (((uint64_t)1) << 40) / 255 + 1; // "1/255" in 24.40 fixed point number
return (M * n) >> 40; // fixed point multiply: n * (1/255)
}
更广泛（需要 64 位 mul），但仍然比 div 指令快。