在 C++ SIMD 中将带符号的 short 转换为 float

Question

我有一个带符号的短数组，我想除以 2048 并得到一个浮点数组。

我发现 SSE: convert short integer to float 允许将 unsigned 短裤转换为浮点数，但我也想处理签名短裤。

下面的代码有效，但仅适用于正面短裤。

// We want to divide some signed short by 2048 and get a float.
const auto floatScale = _mm256_set1_ps(2048);

short* shortsInput = /* values from somewhere */;
float* floatsOutput = /* initialized */;

__m128i* m128iInput = (__m128i*)&shortsInput[0];

// Converts the short vectors to 2 float vectors. This works, but only for positive shorts.
__m128i m128iLow = _mm_unpacklo_epi16(m128iInput[0], _mm_setzero_si128());
__m128i m128iHigh = _mm_unpackhi_epi16(m128iInput[0], _mm_setzero_si128());
__m128 m128Low = _mm_cvtepi32_ps(m128iLow);
__m128 m128High = _mm_cvtepi32_ps(m128iHigh);

// Puts the 2 __m128 vectors into 1 __m256.
__m256 singleComplete = _mm256_castps128_ps256(m128Low);
singleComplete = _mm256_insertf128_ps(singleComplete, m128High, 1);

// Finally do the math
__m256 scaledVect = _mm256_div_ps(singleComplete, floatScale);

// and puts the result where needed.
_mm256_storeu_ps(floatsOutput[0], scaledVect);

如何将我的签名短裤转换为花车？或者也许有更好的方法来解决这个问题？

---

编辑： 与非 SIMD 算法相比，我尝试了不同的答案，在 ~3.2GHz 的 AMD Ryzen 7 2700 上，在 2048 阵列上执行了 10M 次。我正在使用 Visual 15.7.3，主要是默认配置：

/permissive- /Yu"stdafx.h" /GS /GL /W3 /Gy /Zc:wchar_t /Zi /Gm- /O2 /sdl 
/Fd"x64\Release\vc141.pdb" /Zc:inline /fp:precise /D "NDEBUG" /D "_CONSOLE"
/D "_UNICODE" /D "UNICODE" /errorReport:prompt /WX- /Zc:forScope
/arch:AVX2 /Gd /Oi /MD /openmp /FC /Fa"x64\Release\" /EHsc /nologo
/Fo"x64\Release\" /Fp"x64\Release\test.pch" /diagnostics:classic 

请注意，我对 SIMD 非常陌生，并且很长时间没有使用 C++。这是我得到的（我分别重新运行每个测试，而不是一个接一个地重新运行，并获得了更好的结果）：

• 无 SIMD：7300 毫秒
• wim 的回答：2300 毫秒
• chtz 的 SSE2 答案：1650ms
• chtz 的 AVX2 答案：2100ms

因此，我通过使用 SIMD 获得了很好的加速，而 chtz 的 SSE2 答案虽然更冗长且难以理解，但速度更快。 （至少在启用 AVX 的情况下编译时，它避免了使用 3 操作数 VEX 编码指令来复制寄存器的额外指令。在 Intel CPU 上，AVX2 版本应该比 128 位版本快得多。）

这是我的测试代码：

const int size = 2048;
const int loopSize = (int)1e7;

float* noSimd(short* shortsInput) {
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j++) {
            floatsOutput[j] = shortsInput[j] / 2048.0f;
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld noSimd\n", totalTime);

    return floatsOutput;
}

float* wimMethod(short* shortsInput) {
    const auto floatScale = _mm256_set1_ps(1.0f / 2048.0f);
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j += 8) {
            __m128i short_vec = _mm_loadu_si128((__m128i*)&shortsInput[j]);
            __m256i int_vec = _mm256_cvtepi16_epi32(short_vec);
            __m256  singleComplete = _mm256_cvtepi32_ps(int_vec);

            // Finally do the math
            __m256 scaledVect = _mm256_mul_ps(singleComplete, floatScale);

            // and puts the result where needed.
            _mm256_storeu_ps(&floatsOutput[j], scaledVect);
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld wimMethod\n", totalTime);

    return floatsOutput;
}

float* chtzMethodSSE2(short* shortsInput) {
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j += 8) {
            // get input:
            __m128i val = _mm_loadu_si128((__m128i*)&shortsInput[j]);
            // add 0x8000 to wrap to unsigned short domain:
            val = _mm_add_epi16(val, const0x8000);
            // interleave with upper part of float(1<<23)/2048.f:
            __m128i lo = _mm_unpacklo_epi16(val, const0x4580);
            __m128i hi = _mm_unpackhi_epi16(val, const0x4580);
            // interpret as float and subtract float((1<<23) + (0x8000))/2048.f
            __m128 lo_f = _mm_sub_ps(_mm_castsi128_ps(lo), constFloat);
            __m128 hi_f = _mm_sub_ps(_mm_castsi128_ps(hi), constFloat);
            // store:
            _mm_storeu_ps(&floatsOutput[j], lo_f);
            _mm_storeu_ps(&floatsOutput[j] + 4, hi_f);
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld chtzMethod\n", totalTime);

    return floatsOutput;
}

float* chtzMethodAVX2(short* shortsInput) {
    const auto floatScale = _mm256_set1_ps(1.0f / 2048.0f);
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j += 8) {

            // get input:
            __m128i val = _mm_loadu_si128((__m128i*)&shortsInput[j]);
            // interleave with 0x0000
            __m256i val_unpacked = _mm256_cvtepu16_epi32(val);

            // 0x4580'8000
            const __m256 magic = _mm256_set1_ps(float((1 << 23) + (1 << 15)) / 2048.f);
            const __m256i magic_i = _mm256_castps_si256(magic);

            /// convert by xor-ing and subtracting magic value:
            // VPXOR avoids port5 bottlenecks on Intel CPUs before SKL
            __m256 val_f = _mm256_castsi256_ps(_mm256_xor_si256(val_unpacked, magic_i));
            __m256 converted = _mm256_sub_ps(val_f, magic);
            // store:
            _mm256_storeu_ps(&floatsOutput[j], converted);
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld chtzMethod2\n", totalTime);

    return floatsOutput;
}

Answer 1

使用 AVX2 无需单独转换高低片段：

const auto floatScale = _mm256_set1_ps(1.0f/2048.0f);

short* shortsInput = /* values from somewhere */;
float* floatsOutput = /* initialized */;

__m128i short_vec = _mm_loadu_si128((__m128i*)shortsInput);
__m256i int_vec =  _mm256_cvtepi16_epi32 (short_vec);
__m256  singleComplete = _mm256_cvtepi32_ps (int_vec);

// Finally do the math
__m256 scaledVect = _mm256_mul_ps(singleComplete, floatScale);

// and puts the result where needed.
_mm256_storeu_ps(floatsOutput, scaledVect);

这很好地编译了on the Godbolt compiler explorer，并且在 L1d 缓存中输入/输出热和对齐的输入/输出数组，在 Skylake i7-6700k 上以约 360 个时钟周期转换了一个包含 2048 个元素的数组（在重复循环中测试）。每个元素约 0.18 个周期，或每个时钟周期约 5.7 个转换。或者每个向量约 1.4 个周期，包括存储。它主要是前端吞吐量的瓶颈（每个时钟 3.75 个融合域微指令），即使在 clang 的循环展开时也是如此，因为转换是 5 微指令。

请注意，即使在 Haswell/Skylake 上使用简单的寻址模式，vpmovsxwd ymm, [mem] 也无法微融合到单个 uop，因此在这种情况下，最好使用最近的 gcc/clang 将指针增量转换为索引寻址循环计数器。对于大多数内存源向量指令（如vpmovsxwd xmm, [mem]），这将花费额外的微指令：Micro fusion and addressing modes。

一加载一存储，存储不能在Haswell/Skylake的port7存储AGU上运行没关系，它只处理非索引寻址模式。

英特尔 CPU 上的最大吞吐量需要循环展开（如果没有内存瓶颈），因为加载 + 转换 + 存储已经是 4 微指令。与@chtz 的回答相同。

如果您只需要读取几次浮点值，最好立即使用向量结果进行进一步计算。它只有 3 条指令（但对于乱序的 exec 隐藏确实有一些延迟）。在需要时重做转换可能比使用更大的缓存占用空间来将两倍大的float[] 结果存储在内存中要好；这取决于您的用例和硬件。

Answer 2

您可以通过手动组合一个浮点数来替换执行转换 epi16->epi32->float 并乘以 1.f/2048.f 的标准方法。

之所以可行，是因为除数是 2 的幂，因此手动组合浮点数只是意味着不同的指数。

感谢@PeterCordes，这是这个想法的优化 AVX2 版本，使用 XOR 设置 32 位浮点数的高字节，同时翻转整数值的符号位。 FP SUB 将尾数的那些低位转换为正确的 FP 值：

// get input:
__m128i val = _mm_loadu_si128((__m128i*)input);
// interleave with 0x0000
__m256i val_unpacked = _mm256_cvtepu16_epi32(val);

// 0x4580'8000
const __m256 magic = _mm256_set1_ps(float((1<<23) + (1<<15))/2048.f);
const __m256i magic_i = _mm256_castps_si256(magic);

/// convert by xor-ing and subtracting magic value:
// VPXOR avoids port5 bottlenecks on Intel CPUs before SKL
__m256 val_f = _mm256_castsi256_ps(_mm256_xor_si256(val_unpacked, magic_i));
__m256 converted = _mm256_sub_ps(val_f, magic);
// store:
_mm256_storeu_ps(output, converted);

看on the Godbolt compiler explorer with gcc and clang;在 Skylake i7-6700k 上，高速缓存中热的 2048 元素循环需要约 360 个时钟周期，与 @wim 执行标准符号扩展/转换/乘法的版本相同的速度（在测量误差范围内）（具有相似数量的循环展开）。由@PeterCordes 使用 Linux perf 测试。但在 Ryzen 上，这可能会明显更快，因为我们避免了 _mm256_cvtepi32_ps（Ryzen 对 vcvtdq2ps ymm 有 1 per 2 时钟吞吐量：http://agner.org/optimize/。）

0x8000 与下半部分的异或相当于加/减0x8000，因为溢出/进位被忽略。巧合的是，这允许在异或和减法中使用相同的魔法常数。

奇怪的是，gcc 和 clang 更喜欢用加法 -magic 代替减法，这不会重复使用常量......他们更喜欢使用 add，因为它是可交换的，但在这种情况下没有任何好处，因为他们没有将它与内存操作数一起使用。

---

这是一个 SSE2 版本，它与设置 32 位 FP 位模式的高 2 个字节分开进行有符号/无符号翻转。

我们使用一个_mm_add_epi16、两个_mm_unpackXX_epi16 和两个_mm_sub_ps 来表示8 个值（_mm_castsi128_ps 是无操作的，_mm_set 将缓存在寄存器中）：

// get input:
__m128i val = _mm_loadu_si128((__m128i*)input);
// add 0x8000 to wrap to unsigned short domain:
// val = _mm_add_epi16(val, _mm_set1_epi16(0x8000));
val = _mm_xor_si128(val, _mm_set1_epi16(0x8000));  // PXOR runs on more ports, avoids competing with FP add/sub or unpack on Sandybridge/Haswell.

// interleave with upper part of float(1<<23)/2048.f:
__m128i lo = _mm_unpacklo_epi16(val, _mm_set1_epi16(0x4580));
__m128i hi = _mm_unpackhi_epi16(val, _mm_set1_epi16(0x4580));
// interpret as float and subtract float((1<<23) + (0x8000))/2048.f
__m128 lo_f = _mm_sub_ps(_mm_castsi128_ps(lo), _mm_set_ps1(float((1<<23) + (1<<15))/2048.f));
__m128 hi_f = _mm_sub_ps(_mm_castsi128_ps(hi), _mm_set_ps1(float((1<<23) + (1<<15))/2048.f));
// store:
_mm_storeu_ps(output, lo_f);
_mm_storeu_ps(output+4, hi_f);

使用演示： https://ideone.com/b8BfJd

如果您的输入是 unsigned short，则不需要 _mm_add_epi16（当然，需要删除 _mm_sub_ps 中的 1<<15）。然后你会在SSE: convert short integer to float 上看到 Marat 的回答。

这可以轻松移植到 AVX2，每次迭代的转换次数是原来的两倍，但必须注意输出元素的顺序（感谢 @wim 指出这一点）。

---

另外，对于纯 SSE 解决方案，可以简单地使用 _mm_cvtpi16_ps，但这是英特尔库函数。没有一条指令可以做到这一点。

// cast input pointer:
__m64* input64 = (__m64*)input;
// convert and scale:
__m128 lo_f = _mm_mul_ps(_mm_cvtpi16_ps(input64[0]), _mm_set_ps1(1.f/2048.f));
__m128 hi_f = _mm_mul_ps(_mm_cvtpi16_ps(input64[1]), _mm_set_ps1(1.f/2048.f));

我没有对任何解决方案进行基准测试（也没有检查理论吞吐量或延迟）