c++ - 如何使用 SSE/AVX 高效地执行 double/int64 转换？

Question

SSE2 有在单精度 float 和 32 位整数之间转换 vector 的指令。

• _mm_cvtps_epi32()
• _mm_cvtepi32_ps()

但是没有 double 和 64 位整数的等价物。换句话说，它们不见了:

• _mm_cvtpd_epi64()
• _mm_cvtepi64_pd()

AVX 好像也没有。

模拟这些内在函数的最有效方法是什么？

Answer 1

在 AVX512 之前没有单一指令，它添加了与 64 位整数(有符号或无符号)的转换。 (还支持与 32 位无符号的转换)。查看像 _mm512_cvtpd_epi64 这样的内在函数和更窄的 AVX512VL 版本，如 _mm256_cvtpd_epi64 .

如果您只有 AVX2 或更少，则需要以下技巧来进行打包转换。 (对于标量，x86-64 具有来自 SSE2 的标量 int64_t <-> double 或 float，但标量 uint64_t <-> FP 需要技巧，直到 AVX512 添加无符号转换。标量 32 位无符号可以通过零扩展到 64 位来完成签名。)

---

如果你愿意偷工减料，double <-> int64只需两条指令即可完成转换:

• 如果您不关心无穷大或NaN .
• 对于 double <-> int64_t ，您只关心 [-2^51, 2^51] 范围内的值.
• 对于 double <-> uint64_t ，您只关心 [0, 2^52) 范围内的值.

double -> uint64_t

//  Only works for inputs in the range: [0, 2^52)
__m128i double_to_uint64(__m128d x){
    x = _mm_add_pd(x, _mm_set1_pd(0x0010000000000000));
    return _mm_xor_si128(
        _mm_castpd_si128(x),
        _mm_castpd_si128(_mm_set1_pd(0x0010000000000000))
    );
}

double -> int64_t

//  Only works for inputs in the range: [-2^51, 2^51]
__m128i double_to_int64(__m128d x){
    x = _mm_add_pd(x, _mm_set1_pd(0x0018000000000000));
    return _mm_sub_epi64(
        _mm_castpd_si128(x),
        _mm_castpd_si128(_mm_set1_pd(0x0018000000000000))
    );
}

uint64_t -> 双

//  Only works for inputs in the range: [0, 2^52)
__m128d uint64_to_double(__m128i x){
    x = _mm_or_si128(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)));
    return _mm_sub_pd(_mm_castsi128_pd(x), _mm_set1_pd(0x0010000000000000));
}

int64_t -> 双倍

//  Only works for inputs in the range: [-2^51, 2^51]
__m128d int64_to_double(__m128i x){
    x = _mm_add_epi64(x, _mm_castpd_si128(_mm_set1_pd(0x0018000000000000)));
    return _mm_sub_pd(_mm_castsi128_pd(x), _mm_set1_pd(0x0018000000000000));
}

---

舍入行为:

• 对于double -> uint64_t转换，舍入按照当前舍入模式正常工作。 (通常是四舍五入)
• 对于double -> int64_t转换，舍入将遵循除截断之外的所有模式的当前舍入模式。如果当前的舍入模式是截断(向零舍入)，它实际上会向负无穷舍入。

---

它是如何工作的？

尽管这个技巧只有 2 条指令，但并不完全不言自明。

关键是要认识到，对于 double 浮点，[2^52, 2^53) 范围内的值在尾数的最低位下方有“二进制位置”。换句话说，如果将指数位和符号位清零，尾数就变成了整数表示。

转换 x来自 double -> uint64_t ，您添加魔数(Magic Number)M这是 2^52 的浮点值.这把 x进入 [2^52, 2^53) 的“标准化”范围并方便地舍入小数部分。

现在剩下的就是删除高 12 位。这很容易通过掩盖它来完成。最快的方法是识别那些高 12 位与 M 相同。 .因此，我们可以简单地通过 M 减去或异或，而不是引入额外的掩码常数。 . XOR 具有更高的吞吐量。

从 uint64_t -> double 转换只是这个过程的逆过程。您将 M 的指数位加回.然后通过减去 M 来取消归一化数字。 float 。

有符号整数转换稍微复杂一些，因为您需要处理 2 的补码符号扩展。我会把这些留给读者作为练习。

相关: A fast method to round a double to a 32-bit int explained

---

全范围 int64 -> double :

多年后，我终于需要这个了。

• uint64_t -> double 的 5 条说明
• int64_t -> double 的 6 条说明

uint64_t -> 双

__m128d uint64_to_double_full(__m128i x){
    __m128i xH = _mm_srli_epi64(x, 32);
    xH = _mm_or_si128(xH, _mm_castpd_si128(_mm_set1_pd(19342813113834066795298816.)));          //  2^84
    __m128i xL = _mm_blend_epi16(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)), 0xcc);   //  2^52
    __m128d f = _mm_sub_pd(_mm_castsi128_pd(xH), _mm_set1_pd(19342813118337666422669312.));     //  2^84 + 2^52
    return _mm_add_pd(f, _mm_castsi128_pd(xL));
}

int64_t -> 双倍

__m128d int64_to_double_full(__m128i x){
    __m128i xH = _mm_srai_epi32(x, 16);
    xH = _mm_blend_epi16(xH, _mm_setzero_si128(), 0x33);
    xH = _mm_add_epi64(xH, _mm_castpd_si128(_mm_set1_pd(442721857769029238784.)));              //  3*2^67
    __m128i xL = _mm_blend_epi16(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)), 0x88);   //  2^52
    __m128d f = _mm_sub_pd(_mm_castsi128_pd(xH), _mm_set1_pd(442726361368656609280.));          //  3*2^67 + 2^52
    return _mm_add_pd(f, _mm_castsi128_pd(xL));
}

这些适用于整个 64 位范围，并正确舍入为当前舍入行为。

这些是下面类似的 wim 的答案 - 但有更多的滥用优化。因此，解读这些内容也将留给读者作为练习。