如何在 sse 中实现有符号定点数学中向零的衰减？

Posted 2023-02-16

【中文标题】如何在 sse 中实现有符号定点数学中向零的衰减？

【英文标题】：How to implement decay towards zero in signed fixed point math, in sse?

【发布时间】：2018-10-03 15:34:22

【问题描述】：

有许多类似衰减的物理事件（例如体摩擦或电荷泄漏），通常在像x' = x * 0.99这样的迭代器中建模，这通常非常易于用浮点算术编写。

但是，我需要在 sse 中以 16 位“8.8”签名定点方式执行此操作。为了在典型 ALU 上有效实现，提到的公式可以重写为 x = x - x/128; 或 x = x - (x>>7) 其中 >> 是“算术”，符号扩展右移。

我卡在这里，因为_mm_sra_epi16() 产生完全违反直觉的行为，这很容易通过以下示例验证：

#include <cstdint>
#include <iostream>
#include <emmintrin.h>

using namespace std;

int main(int argc, char** argv) 
    cout << "required: ";
    for (int i = -1; i < 7; ++i) 
        cout << hex << (0x7fff >> i) << ", ";
    
    cout << endl;
    cout << "produced: ";
    __m128i a = _mm_set1_epi16(0x7fff);
    __m128i b = _mm_set_epi16(-1, 0, 1, 2, 3, 4, 5, 6);
    auto c = _mm_sra_epi16(a, b);
    for (auto i = 0; i < 8; ++i) 
        cout << hex << c.m128i_i16[i] << ", ";
    
    cout << endl;
    return 0;

输出如下：

required: 0, 7fff, 3fff, 1fff, fff, 7ff, 3ff, 1ff,
produced: 0, 0, 0, 0, 0, 0, 0, 0,

它只适用于所有的第一个转变，就像它实际上是_mm_sra1_epi16 函数一样，意外地命名为sra 并给__m128i 第二个参数提供了一个无缘无故的有趣子句。所以这不能在 SSE 中使用。

另一方面，我听说除法算法非常复杂，因此 _mm_div_epi16 在 SSE 中不存在，也无法使用。 该做什么以及如何实现/矢量化这种流行的“衰减”技术？

【问题讨论】：

_mm_sra_epi16 使用第二个源向量的低 64 位作为适用于所有元素的 64 位移位计数。这不是白痴，但每个元素的移位计数需要 AVX2（用于 32/64 位元素）或 AVX512BW 用于_mm_srav_epi16 或 64 位算术右移，这对于您尝试使用它的方式是有意义的。 （但移位计数是无符号的，所以-1 也将移出所有位）。

是_mm_sra_epi16()“白痴”，因为它不是一个可变的转变？它基于最低通道将所有元素移动相同的量。如果你想要一个 16 位宽度的可变移位，你需要 AVX512。

@Mysticial 为什么_mm_add_epi16() 没有“基于最低车道”添加相同的数量？我认为它的行为不足并且违反了一致的接口，从而破坏了移植/矢量化。实际上，该指令应该命名为_mm_sra1_epi16()。

@PeterCordes 那么你认为这个任务没有解决方案吗？

为什么您认为每个元素需要不同的班次计数？

【参考方案1】：

x -= x>>7 使用 SSE2 实现是微不足道的，使用恒定的移位计数来提高效率。如果 AVX 可用，这将编译为 2 条指令，否则需要 movdqa 在破坏性右移之前复制 v。

__m128i downscale(__m128i v)
    __m128i dec = _mm_srai_epi16(v, 7);
    return _mm_sub_epi16(v, dec);

GCC 甚至自动对其进行矢量化 (Godbolt)。

void foo(short *__restrict a) 
    for (int i=0 ; i<10240 ; i++) 
        a[i] -= a[i]>>7;  // inner loop uses the same psraw / psubw
    

---

与float 不同，定点在整个范围内具有恒定的绝对精度，而不是恒定的相对精度。因此，对于小的正数，v>>7 将为零，您的递减将停止。 （负输入下溢到-1，因为算术右移向 -infinity 舍入。）

如果移位可能下溢到 0 的小输入，您可能需要与 _mm_set1_epi16(1) 进行 OR 以确保递减量不为零。对大输入的影响可以忽略不计。但是，这最终会使缩减链从 0 变为 -1。 （然后回到 0，因为 -1 | 1 == -1 在 2 的补码中）。

__m128i downscale_nonzero(__m128i v)
    __m128i dec = _mm_srai_epi16(v, 7);
    dec = _mm_or_si128(dec, _mm_set1_epi16(1));
    return _mm_sub_epi16(v, dec);

如果从负数开始，序列将是 -large，对数直到 -128，线性直到 -4、-3、-2、-1、0、-1、0、-1，...

---

您的代码全为零，因为_mm_sra_epi16 使用第二个源向量的低 64 位作为适用于所有元素的 64 位移位计数。 Read the manual。因此，您将所有位从每个 16 位元素中移出。

这不是白痴，但每个元素的移位计数需要 AVX2（用于 32/64 位元素）或 AVX512BW 用于 _mm_srav_epi16 或 64 位算术右移，这对于您尝试的方式是有意义的用它。 （但移位计数是无符号的，所以-1 也将移出所有位）。

确实，该指令应该命名为_mm_sra1_epi16()

是的，这是有道理的。但请记住，当这些被命名时，AVX2 _mm_srav_* 还不存在。此外，该特定名称并不理想，因为 1 和 i 在视觉上并不是最明显的。 （i 是立即数，对于 psraw xmm1, imm16 形式而不是 asm 指令的 psraw xmm1, xmm2/m128 形式：http://felixcloutier.com/x86/PSRAW:PSRAD:PSRAQ.html）。

另一种有意义的方式是 MMX/SSE2 asm 指令有两种形式：立即数（当然，所有元素的计数相同）和向量。向量版本不是强迫您将计数广播到所有元素，而是采用向量寄存器底部的标量计数。我认为预期的用例是在 movd xmm0, eax 之后。

---

如果您需要不使用 AVX512 的每个元素变量的移位计数，请参阅有关模拟它的各种问答，例如Shifting 4 integers right by different values SIMD.

一些变通方法使用乘以 2 的幂进行可变左移，然后右移将数据放在需要的位置。 （但是您需要以某种方式准备好1<<n SIMD 向量，因此如果将同一组计数用于许多向量，或者特别是如果它是编译时常量，则此方法有效。

对于 16 位元素，您可以只使用一个 _mm_mulhi_epi16 来执行运行时变量右移计数，而不会造成精度损失或范围限制。 mulhi(x*y) 与 (x*(int)y) >> 16 完全相同，因此您可以使用 y=1<<14 在该元素中右移 16-14 = 2。

【讨论】：

我知道了。刚刚看到一个任务x -= x >> y，在 Haswell 之前无法有效地矢量化。很高兴他们添加了指令，但很惊讶为什么通用 sra 没有与 add、mul、and、or 基本指令一起添加。添加的是sra1，它被错误地命名为sra，并且用例更少。

我认为我做错了，我错过了另一种在这种情况下使用的衰减算法。至于位，它不需要收敛到零，因为 x*=.99 永远不会为零。

@xakepp35 你不是唯一一个注意到 SIMD 指令不一致并且遗漏很多情况的人。大多数问题已在 AVX512 中得到修复。迟到了，但至少英特尔有人注意到了。

AFAICT 现在，AVX512 中的大多数“明显但缺失”的功能要么没人使用，要么难以在硬件中实现。所以他们终于做对了。

我解决了！ _mm_sub_epi16(x, _mm_srai_epi16(x, a))几乎等于_mm_slli_epi16(_mm_mulhi_epi16(x, b), 1)，其中b = (1<<15) - (1<<(15-a))..几乎-是因为slli。虽然失去了一点精度，但第二种情况允许更精细的参数调整。