如何使用SSE将16位整数除以255?

13

我处理图像处理。

我需要将16位整数SSE向量除以255。

由于255不是2的幂的倍数,因此我无法使用移位运算符如_mm_srli_epi16()。

当然,我知道可以将整数转换为浮点数,进行除法操作,然后再将其转换回整数。

但也许有人知道另一种解决方案...


1
这个有帮助吗? - Anton Savin
1
通常情况下,您会除以256(四舍五入而不是截断)-是否有什么原因必须是255而不是256? - Paul R
1
也许对你来说这个问题很有趣。当你将来需要处理非常数整数除法时,转换为浮点数也是一种快速的选择。 - Youka
1
@Paul-R 如果我除以256,这将导致平均亮度降低,当执行此操作多次时。 - Claudio
@Claudio:没有看到你的其它代码,很难说。不过如果你正在使用完整的16位范围,然后缩小到8位,那么比例因子应该是256。否则,如果你取例如完整刻度uint16_t = 65535,然后除以255,你会得到257。我猜你知道自己在做什么-我只是感到困惑... - Paul R
我更新了我的答案,提供了一个无符号版本,其吞吐量比已接受的答案更好。之前我没有意识到你在使用无符号,但我认为你是基于srli而不是srai - Peter Cordes
5个回答

16

有一种整数近似除以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(Power):

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);
}

对于 NEON(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);
}

1
这对于 value == 65535 和所有负数都是错误的(因此既不适用于有符号 16 位整数也不适用于无符号 16 位整数)。 - Anton Savin
1
我知道它在 alpha 混合方面完美运作。但我不排除在其他情况下可能会出现错误。 - ErmIg
@AntonSavin:通过除以255来实现精确的无符号除法比使用pmulhuw更有效,尽管延迟较差。我几年前的答案只考虑了有符号数,而没有考虑无符号数。(由于某种原因,似乎没有人区分这两者,有时使用int,有时使用SIMD逻辑右移。) - Peter Cordes

10
如果您想在所有情况下得到完全正确的结果,请按照Anton链接问题中Marc Glisse的评论中的建议进行操作:SSE整数除法? 使用GNU C原生矢量语法来表示矢量除以给定标量,并在Godbolt编译器浏览器上查看其运行结果

无符号除法是廉价的:

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);
}

仅使用2个微操作,如果您的瓶颈在前端吞吐量或Intel CPU上的端口0吞吐量,则此方法的吞吐量更高(但延迟较差),与@ermlg的答案相比。 (始终取决于在将其用作更大功能的一部分时周围代码的情况。)http://agner.org/optimize/ 在英特尔芯片上,向量移位仅在端口0上运行,因此@ermlg的2次移位+1次加法在端口0上成为瓶颈。(再次取决于周围代码)。而且它是3个微操作,而这个方法只有2个。
在Skylake上,pmulhuw / pmulhw在端口0或1上运行,因此可以与移位并行运行。(但在Broadwell及更早版本中,它们仅在端口0上运行,与移位冲突。因此,在英特尔Skylake之前,唯一的优点是前端和无序执行的总微操作数更少,以便跟踪。)在英特尔上,pmulhuw的延迟为5个周期,与移位的1个周期相比较高,但是当您可以为更高的吞吐量节省微操作时,OoO exec通常可以隐藏更多的延迟。
Ryzen也仅在其P0上运行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足够聪明,使用相同的16B常量存储在内存中用于set1和生成的div255常量。据我所知,这类似于字符串常量合并。

3

GCC会优化x/255,其中xunsigned 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)

3

准确版本:

#define div_255_fast(x)    (((x) + (((x) + 257) >> 8)) >> 8)

当x在[0, 65536]的范围内时,ERROR为0。它比快两倍:

enter image description here

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位乘法),但仍比div指令快。


另外,我认为你的意思是 [0, 65535] 或者 [0, 65536),因为 65536 = 2^16 超出了 epu16 的范围。 - Peter Cordes
@PeterCordes,已更新,[0, 65536]适用于C版本。 - skywind3000
在普通的C标量中,您不需要特殊技巧。编译器已经知道如何使用乘法逆元:为什么GCC在实现整数除法时使用奇怪的数字进行乘法?。只需使用“unsigned”,您就可以获得高效的汇编代码,或者像您的带符号版本一样的汇编代码,适用于任何“int32_t”:https://godbolt.org/z/cTzYzb。您的“fast_div_255_any”似乎处理了符号位,因此它实际上适用于任何x,而不仅仅是正数,对吗?我想这是将其转换为SIMD(`_mm_mul_epu32`(`pmuludq`))的有用步骤。 - Peter Cordes
@PeterCordes,div_255_fast比fast_div_255_any便宜得多(gcc可能会使用它)。如果您的数字在0到65536之间,它仍然比编译器快。 - skywind3000
如果你是编译针对Silvermont或Atom的代码,那么速度会明显更快,因为64位imul运算较慢。而且,如果你做了什么像使用int32_t x这样的傻事来进行x/255的有符号除法并迫使编译器处理它,那也会更慢。 - Peter Cordes
显示剩余5条评论

2

出于好奇(并且如果性能是一个问题),以下是使用 (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 );

通过饱和加法,即使对于接近环绕点的值,您也可以使用此功能而不会产生灾难性的结果。 _mm_adds_epi16epu16 用于有符号/无符号饱和。 您确定要使用 srli(逻辑右移)而不是算术右移 srai 吗? 还是您正在为无符号执行此操作? - Peter Cordes

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接