将整数向量转换为介于0和1之间的浮点数的最快准确方法

Question

考虑一个随机生成的__m256i向量。是否有一种更快、更精确的方法将它们转换为__m256浮点数向量，范围在0（包括）和1（不包括）之间，而不是通过float(1ull<<32)进行除法运算？

到目前为止，我尝试了以下内容，其中iRand是输入，ans是输出：

const __m256 fRand = _mm256_cvtepi32_ps(iRand);
const __m256 normalized = _mm256_div_ps(fRand, _mm256_set1_ps(float(1ull<<32)));
const __m256 ans = _mm256_add_ps(normalized, _mm256_set1_ps(0.5f));

- Serge Rogatch

---

乘以 0x1p-31f。通常情况下，除法比乘法更慢或需要更多资源。 - Eric Postpischil

由于您正在立即添加，建议查看SIMD融合乘加指令。 - Eric Postpischil

@EricPostpischil，非常感谢。我本来想用一些位操作的魔法，但是你简单的建议仍然比我的原始版本快得多。 - Serge Rogatch

当适用时，我使用位操作技巧。有时，当位数少于24位时，您可以将它们补丁到float的有效数字字段中，并使用浮点运算完成工作。但是，您显然支持31位（和一个符号），因此必须四舍五入。转换指令是为此设计的，因此您不太可能做得更好。 - Eric Postpischil

@chtz，你只需要计算 sqrt(-2log(1-x))。 - Severin Pappadeux

Answer 1

以下版本应该比您最初使用的 _mm256_div_ps 版本更快。 vdivps非常缓慢，例如，在我的Haswell Xeon上，它的延迟为18-21个周期，吞吐量为14个周期。顺便提一句，较新的CPU表现更好，Skylake为11/5，Ryzen为10/6。

如评论中所述，通过将除法替换为乘法，并进一步改进FMA，可以解决性能问题。 这种方法的问题在于分布的质量。 如果您尝试通过四舍五入或剪裁来获得输出区间中的这些数字，则会在输出数字的概率分布中引入峰值。

我的实现也不是完美的，它不能在输出区间内输出所有可能的值，跳过了许多可表示的浮点数，特别是接近0的数。 但至少分布非常均匀。

__m256 __vectorcall randomFloats( __m256i randomBits )
{
    // Convert to random float bits
    __m256 result = _mm256_castsi256_ps( randomBits );

    // Zero out exponent bits, leave random bits in mantissa.
    // BTW since the mask value is constexpr, we don't actually need AVX2 instructions for this, it's just easier to code with set1_epi32.
    const __m256 mantissaMask = _mm256_castsi256_ps( _mm256_set1_epi32( 0x007FFFFF ) );
    result = _mm256_and_ps( result, mantissaMask );

    // Set sign + exponent bits to that of 1.0, which is sign=0, exponent=2^0.
    const __m256 one = _mm256_set1_ps( 1.0f );
    result = _mm256_or_ps( result, one );

    // Subtract 1.0. The above algorithm generates floats in range [1..2).
    // Can't use bit tricks to generate floats in [0..1) because it would cause them to be distributed very unevenly.
    return _mm256_sub_ps( result, one );
}

更新：如果您希望获得更好的准确性，请使用以下版本。但它不再是“最快”的。

__m256 __vectorcall randomFloats_32( __m256i randomBits )
{
    // Convert to random float bits
    __m256 result = _mm256_castsi256_ps( randomBits );
    // Zero out exponent bits, leave random bits in mantissa.
    const __m256 mantissaMask = _mm256_castsi256_ps( _mm256_set1_epi32( 0x007FFFFF ) );
    result = _mm256_and_ps( result, mantissaMask );
    // Set sign + exponent bits to that of 1.0, which is sign=0, exponent = 2^0.
    const __m256 one = _mm256_set1_ps( 1.0f );
    result = _mm256_or_ps( result, one );
    // Subtract 1.0. The above algorithm generates floats in range [1..2).
    result = _mm256_sub_ps( result, one );

    // Use 9 unused random bits to add extra randomness to the lower bits of the values.
    // This increases precision to 2^-32, however most floats in the range can't store that many bits, fmadd will only add them for small enough values.

    // If you want uniformly distributed floats with 2^-24 precision, replace the second argument in the following line with _mm256_set1_epi32( 0x80000000 ).
    // In this case you don't need to set rounding mode bits in MXCSR.
    __m256i extraBits = _mm256_and_si256( randomBits, _mm256_castps_si256( mantissaMask ) );
    extraBits = _mm256_srli_epi32( extraBits, 9 );
    __m256 extra = _mm256_castsi256_ps( extraBits );
    extra = _mm256_or_ps( extra, one );
    extra = _mm256_sub_ps( extra, one );
    _MM_SET_ROUNDING_MODE( _MM_ROUND_DOWN );
    constexpr float mul = 0x1p-23f; // The initial part of the algorithm has generated uniform distribution with the step 2^-23.
    return _mm256_fmadd_ps( extra, _mm256_set1_ps( mul ), result );
}

- Soonts

---

可能不会比Eric在Intel CPU上的cvt + FMA建议更快，但我认为可以避免Antti指出的convert中大整数四舍五入超过+-2^24范围的正确性问题。在只有一个FMA单元的AMD CPU上，这也可能更快。或者在Intel上，如果周围的代码在FMA / mul / add上瓶颈并使端口5空闲，则甚至可能更快。 - Peter Cordes

@PeterCordes 对，我所说的“更快”是指与原始OP使用div_ps相比。虽然它不会在输出间隔中输出所有可能的浮点数，但是这种方法的分布很好。如果您使用_MM_ROUND_TOWARD_ZERO在MXCSR中解决cvtepi32_ps问题，则分布将略微偏向0.5。 - Soonts

哦，关于分布的观点很好。这是在实数范围内均匀分布的，而不是在可表示的浮点位模式上分布。在答案中说一些这方面的内容是个好主意，还要说明为什么它更快（因为div非常慢），并且额外加分的话可以与FMA进行比较。因为除了速度之外，这里还有许多有趣的问题。 - Peter Cordes

从OP在评论中澄清的内容来看，我猜测使用convert+FMA会更好，因为它生成的值更接近于0.0f（并且浪费较少的随机整数位）。对于Box-Muller来说，生成1.0f不应该是一个问题（但确切地生成0.0f则有可能）。 - chtz

@PeterCordes 是的，有很多有趣的问题。我发布了另一个转换程序的版本，欢迎您的评论。 - Severin Pappadeux

Answer 2

首先，没有除法，用乘法代替。虽然@Soonts对你来说可能足够好，但由于使用映射到[1...2)区间，它会产生形式为k/2⁻²³的均匀二进制有理数，这只有可能生成一半。我更喜欢S.Vigna（底部）的方法，其中所有形式为k/2⁻²⁴的二进制有理数都是等可能的。

代码，VC++2019，x64，Win10，Intel i7 Skylake

#include <random>

#include "immintrin.h"

auto p256_dec_u32(__m256i in) -> void {
    alignas(alignof(__m256i)) uint32_t v[8];
    _mm256_store_si256((__m256i*)v, in);
    printf("v8_u32: %u %u %u %u %u %u %u %u\n", v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
}

auto p256_dec_f32(__m256 in) -> void {
    alignas(alignof(__m256)) float v[8];
    _mm256_store_ps(v, in);
    printf("v8_float: %e %e %e %e %e %e %e %e\n", v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
}

auto main() -> int {
    const float c = 0x1.0p-24f; // or (1.0f / (uint32_t(1) << 24));

    const int N = 1000000;

    std::mt19937 rng{ 987654321ULL };

    __m256 sum = _mm256_set1_ps(0.0f);

    for (int k = 0; k != N; ++k) {
        alignas(alignof(__m256i)) uint32_t rnd[8] = { rng(), rng(), rng(), rng(), rng(), rng(), rng(), rng() };

        __m256i r = _mm256_load_si256((__m256i*)rnd);
        __m256  q = _mm256_mul_ps(_mm256_cvtepi32_ps(_mm256_srli_epi32(r, 8)), _mm256_set1_ps(c));

        sum = _mm256_add_ps(sum, q);
    }

    sum = _mm256_div_ps(sum, _mm256_set1_ps((float)N)); // computing average

    p256_dec_f32(sum);

    return 0;
}

带有输出

5.002970e-01 4.997833e-01 4.996118e-01 5.004955e-01 5.002163e-01 4.997193e-01 4.996586e-01 5.001499e-01

- Severin Pappadeux

---

这是单精度。k/2^-52将用于__m256d。 - Peter Cordes

为什么向右移动8位？我能理解移动一位来去掉符号位，但你不必要地去掉7个随机位（这些位可以用于生成0到2^-24之间的数字）。 - chtz

@chtz 是的，这种方法只产生形式为k/2^−24的值，但是它们是均匀分布的，一些随机输入会丢失。为了获得所有可能的浮点数，您需要像http://mumble.net/~campbell/2014/04/28/uniform-random-float那样做，并查看S.Vigna页面的底部。 - Severin Pappadeux