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

Question

我想在SSE(像Izhikevich尖峰神经元模型这样的程序)中实现一个简单的函数.它应该使用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()指令,检查乘法的高位字.为什么以及如何v1_overflows_during_this_operation()在SSE中实施此类任务？

Answer 1

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

相反，_mm_mulhi_epi16and_mm_mulhi_epu16只给出完整结果的有符号或无符号上半部分。

（并且_mm_mullo_epi16，它执行打包 16x16 => 16 位低半截断乘法，这对于有符号或无符号相同）。

您可以使用_mm_unpacklo/hi_epi16将低/高半部分交错成一对具有 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)。

那么您认为不让结果溢出，而只是像_mm_cmpge_epi16(v1, vThreshold)作者在原始论文中那样生成掩码更容易吗？

当然好！获得另外一两位精度可能会导致性能损失 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。