如何检查 SSE 中 16 位整数乘法的溢出?

How to check overflow for multiplication of 16 bit integers in SSE?

我想在SSE(类似Izhikevich spiking neuron model的程序)中实现一个简单的功能。它应该使用 16 位有符号整数(8.8 定点),并且需要在某些积分步骤中检查溢出情况,并设置 SSE 掩码(如果发生溢出):

// initialized like following:
short I = 0x1BAD; // current injected to neuron
short vR = 0xF00D; // some reset threshold when spiked (negative)

// step to be vectorized:
short v0 = vReset;
for(;;) {

    // v0*v0/16 likely overflows => use 32 bit (16.16)
    short v0_sqr = ((int)v0)*((int)v0) / (1<<(8+4)); // not sure how "(v0*v0)>>(8+4)" would affect sign..
     // or   ((int)v0)*((int)v0) >> (8+4); // arithmetic right shift
     // original paper used v' = (v0^2)/25 + ...

    short v1 = v0_sqr + v0 + I;
    int m; // mask is set when neuron fires
    if(v1_overflows_during_this_operation()) { // "v1 > 0x7FFF" - way to detect?
        m=0xFFFFFFFF;
    else
        m=0;
    v0 = ( v1 & ~m ) | (vR & m );
}

但是我没有找到_mm_mul_epi16()指令,检查乘法的高字。为什么以及如何在 SSE 中实施这样的任务 v1_overflows_during_this_operation()

与 32x32 => 64 不同,没有加宽 16x16 -> 32 SSE 乘法指令。

相反,_mm_mulhi_epi16 and _mm_mulhi_epu16 只为您提供完整结果的带符号或不带符号的上半部分。

(和 _mm_mullo_epi16,它确实打包 16x16 => 16 位低半截断乘法,这对于有符号或无符号都是相同的)。

您可以使用 _mm_unpacklo/hi_epi16 将 low/high 的两半交织成一对具有 32 位元素的向量,但这会很慢。但是,是的,您可以 _mm_srai_epi32(v, 8+4) 算术右移 12,然后重新打包,可能使用 _mm_packs_epi32(带符号的饱和度回到 16 位)。然后我想检查饱和度?


您的用例不寻常。 _mm_mulhrs_epi16 为您提供高 17 位,四舍五入,然后截断为 16 位。 (见说明)。这对于某些定点算法很有用,在这些算法中输入被缩放以将结果放在上半部分,并且您希望四舍五入包括下半部分而不是截断。

您可能实际上使用 _mm_mulhrs_epi16_mm_mulhi_epi16 作为保持最精确度的最佳选择,也许通过在平方前左移您的 v0 到高点一半会给你 (v0*v0) >> (8+4).

So do you think it is easier not to allow result to overflow, and just to generate mask with _mm_cmpge_epi16(v1, vThreshold) as author does in the original paper?

该死!获得另一位或两位精度可能会使性能损失 2 倍,因为您必须计算另一个乘法结果来检查溢出,或者有效地扩大到 32 位(将每个向量的元素数量减半),如上所述。

通过比较结果,v0 = ( v1 & ~m ) | (vR & m ); 成为 SSE4.1 混合:_mm_blendv_epi8


如果您的 vThreshold 在顶部有 2 个未设置的位,您就有空间左移而不会丢失任何最高有效位。由于 mulhi 给你 (v0*v0) >> 16,所以你可以这样做:

// losing the high 2 bits of v0
__m128i v0_lshift2   = _mm_slli_epi16(v0, 2);    // left by 2 before squaring
__m128i v0_sqr_asr12 = _mm_mulhi_epi16(v0_lshift2, v0_lshift2);
__m128i v1 = _mm_add_epi16(v0, I);
        v1 = _mm_add_epi16(v1, v0_sqr_asr12);

    // v1 = ((v0<<2)* (int)(v0<<2))) >> 16) + v0 + I

    // v1 = ((v0*(int)v0) >> 12) + v0 + I

平方前左移 2 与平方后左移 4 相同(完整的 32 位结果)。它将我们想要的16位准确地放入了高16位。

但是,如果您的 v0 非常接近全范围,以至于您在左移时可能溢出。

否则可以在乘

之前丢掉v0的6个低位

使用算术右移向 -Infinity 舍入会损失 6 位精度,但不可能溢出。

// losing the low 6 bits of v0
__m128i v0_asr6 = _mm_srai_epi16(v0, 6);
__m128i v0_sqr_asr12 = _mm_mullo_epi16(v0_asr6, v0_asr6);
__m128i v1 = _mm_add_epi16(v0, I);
        v1 = _mm_add_epi16(v1, v0_sqr_asr12);

    // v1 =  (v0>>6) * (int)(v0>>6)) + v0 + I

    // v1 ~= ((v0*(int)v0) >> 12) + v0 + I

我认为您会以这种方式损失更多的精度,因此最好将 vThreshold 设置得足够小,以便您有足够的开销来使用高半乘法。这种方式包括可能更差的舍入。

pmulhrsw 舍入而不是截断可能会更好,如果我们可以有效地设置它的话。但我不认为我们可以,因为右移 1 是奇数。我认为我们需要进行 2 个单独的输入,一个 v0_lshift2 和一个仅左移 1.