【问题标题】:Map integer range onto another range将整数范围映射到另一个范围
【发布时间】:2020-03-03 11:48:02
【问题描述】:

在运行时,我有 2 个由它们的 uint32_t 边界 a..bc..d 定义的范围。第一个范围往往比第二个大得多:8 < (b - a) / (d - c) < 64
确切限制:a >= 0b <= 2^31 - 1c >= 0d <= 2^20 - 1

我需要一个例程来执行从第一个范围到第二个范围的整数的线性映射:f(uint32_t x) -> round_to_uint32_t((float)(x - a) / (b - a) * (d - c) + c)
b - a >= d - c 时,保持比率尽可能接近理想是很重要的,否则如果来自[a; b] 的元素可以映射到来自[c; d] 的多个整数上,则可以返回这些整数中的任何一个。

听起来像是一个简单的比率问题,并且已经在许多问题中得到了解答,例如
Convert a number range to another range, maintaining ratio
但在这里我需要一个非常快速的解决方案。

此例程是专用排序算法的关键部分,对于已排序数组的每个元素至少会调用一次。

SIMD 解决方案在不降低整体性能的情况下也是可以接受的。

【问题讨论】:

  • “不 AVX”是什么意思? _mm256_fmadd_ps 看起来是显而易见的方法,为乘数 (d - c) / (b - a) 预先计算循环不变值。您可以对连续内存中的一组输入有效地执行此操作吗?如果不是,那么标量 fma() 仍然很好,并简化了 uint32_t 浮点转换。 (您可以只在标量 64 位整数之间进行转换,并将低半部分视为无符号)。
  • 精确限制:a >= 0b <= 2^31 - 1c >= 0d <= 2^20 - 1。它总是可能为零c
  • 完美,因此您可以将输入视为有符号正数或无符号数,并且x-a 始终是正数符号,而不是仍处于无符号范围的上半部分。并且可以四舍五入到单精度float?将x 舍入到float 而不是x-a 可以吗?这让您可以使用 cvt + vfmaddps + cvt 来完成(因为 -a * scale 可以在末尾与 +c 结合使用)。
  • IDK 你所说的“它们是 only 无符号整数”是什么意思。 uint32_t 的全范围是0 .. 2^32 - 1。如果x-a 可能 >= INT_MAX,您将无法有效地实现这一点。因此,如果您的数字被解释为int32_t(即它们的高位未设置),那么您的数字具有相同的值是一件好事。无论如何,我了解实际范围是多少,我只是在您的后续评论中对术语进行挑剔。
  • 如果部分 C++ 代码可以矢量化,那么您在 RAX 中就不需要它了。这就是为什么我问如何你使用这个值(特别是作为内存地址的一部分)。显然,这个重要问题的答案是肯定的,所以你在 RAX(或另一个整数 reg)中确实需要它,但是在 C++ 中多次使用它并不会自动遵循。你标记了这个程序集;我正在考虑我们希望编译器生成什么代码。 (因为这就是你编写高效 C++ 的方式)

标签: c++ assembly optimization x86-64 simd


【解决方案1】:

实际运行时除法(FP 和整数)非常慢,因此您绝对希望避免这种情况。您编写该表达式的方式可能编译为包含除法,因为 FP 数学不是关联的(没有-ffast-math);编译器无法为您将x / foo * bar 转换为x * (bar/foo),即使循环不变bar/foo 非常好。您确实需要浮点或 64 位整数来避免乘法溢出,但只有 FP 允许您重用非整数循环不变的除法结果。

_mm256_fmadd_ps 看起来是显而易见的方法,为乘数 (d - c) / (b - a) 提供了一个预先计算的循环不变值。如果float 舍入对于严格按顺序进行(乘然后除)不是问题,那么在循环之外首先进行这种不精确的除法可能是可以的。喜欢
_mm256_set1_ps((d - c) / (double)(b - a))。使用double 进行此计算可避免在转换为除法操作数的 FP 期间出现舍入错误。

您对许多 x 重复使用相同的 a、b、c、d,可能来自连续内存。您将结果用作内存地址的一部分,因此不幸的是,您最终确实需要将 SIMD 中的结果返回到整数寄存器中。 (可能使用 AVX512 分散存储可以避免这种情况。)

现代 x86 CPU 具有 2 个时钟的负载吞吐量,因此将 8x uint32_t 返回整数寄存器的最佳选择可能是向量存储/整数重新加载,而不是为 ALU shuffle 的每个元素花费 2 微秒。这有一些延迟,所以我建议在循环通过该标量之前转换为 16 或 32 个整数(64 或 128 字节)的 tmp 缓冲区,即 2x 或 4x __m256i

或者可能交替转换和存储一个向量,然后循环遍历您之前转换的 another 的 8 个元素。即软件流水线。乱序执行可以隐藏延迟,但您已经将其延迟隐藏功能扩展到缓存未命中,无论您使用内存做什么。

根据您的 CPU(例如 Haswell 或某些 Skylake),使用 256 位向量指令可能会限制您的最大涡轮增压比其他情况略低。您可能会考虑一次只做 4 个向量,但是每个元素会花费更多的微指令。

如果不是 SIMD,那么即使是标量 C++ fma() 对于 vfmadd213sd 来说仍然很好,但是使用内部函数是从 float -> int (vcvtps2dq 而不是vcvttps2dq)。


请注意,uint32_t float 转换直到 AVX512 才直接可用。对于标量,您可以在 int64_t 之间进行转换,并使用无符号低半部分的截断/零扩展。

非常方便(如 cmets 中所讨论的)您的输入是有范围限制的,因此如果您将它们解释为有符号整数,它们具有相同的值(有符号非负数)。已知xx-a(和b-a)都是正数并且0x7FFFFFFF。 (或者至少是非负数。零很好。)

浮点舍入

对于 SIMD,单精度 float非常有利于 SIMD 吞吐量。与签名 int32_t 的高效打包转换。但是not every int32_t can be exactly represented as a float。较大的值会四舍五入到最接近的偶数,最接近 2^2、2^3 的倍数,或者越远高于 2^24 的值。

使用 SIMD double 是可能的,但需要一些改组。

我不认为float 对于用(float)(x-a) 编写的公式通常不是问题。如果b-a 输入范围很大,这意味着两个范围都很大,并且舍入误差不会将所有可能的x 值映射到相同的输出中。根据乘数的不同,输入舍入误差可能比输出舍入误差更差,对于更高的x-a 值,可能会留下一些可表示的输出浮点数。

但是,如果我们想将-a * (d - c) / (b - a) 部分分解出来并在最后将其与+c 结合起来,那么

  1. 在要添加的值中,我们可能会因灾难性取消而损失精度。
  2. 我们需要对原始输入值执行(float)x。如果a 很大而b-a 很小,即接近可能输入范围顶部的小范围,则舍入误差可以将所有可能的x 值映射到相同的浮点数。
  3. 为了充分利用 FMA,我们希望在转换回整数之前执行 +c,如果 d-c 是一个小输出范围但 c 是巨大的,这又会带来输出舍入错误的风险。在你的情况下不是问题;使用d 20 - 1 我们知道float 可以准确地表示该 c..d 范围内的每个输出整数值。

如果您没有输入范围限制,您可以在缩放之前通过在输入上使用整数 (x-a)+0x80000000U 和在输出上使用 ...+c+0x80000000U(四舍五入到最接近的 int32_t 之后)进行范围转换。但这会为小的 uint32_t 输入(接近 0)引入巨大的 float 舍入误差,这些输入将范围移动到接近 INT_MIN

我们不需要对 b-ad-c 进行范围移位,因为 + 或 - 或与 0x80000000U 的 XOR 会在减法中抵消。


示例:

在这个内联之后,const 向量应该被编译器从循环中提升出来, 或者您可以手动执行此操作。

这需要 AVX1 + FMA(例如 AMD Piledriver 或 Intel Haswell 或更高版本)。未经测试,抱歉我什至没有把它扔到 Godbolt 上看看它是否编译。

// fastest but not safe if b-a is small and  a > 2^24
static inline
__m256i range_scale_fast_fma(__m256i data, uint32_t a, uint32_t b, uint32_t c, uint32_t d)
{
     // avoid rounding errors when computing the scale factor, but convert double->float on the final result
    double scale_scalar = (d - c) / (double)(b - a);
    const __m256 scale = _mm256_set1_ps(scale_scalar);
    const __m256 add = _m256_set1_ps(-a*scale_scalar + c);
    //    (x-a) * scale + c
    // =  x * scale + (-a*scale + c)   but with different rounding error from doing -a*scale + c

    __m256  in = _mm256_cvtepi32_ps(data);
    __m256  out = _mm256_fmadd_ps(in, scale, add);
    return _mm256_cvtps_epi32(out);   // convert back with round to nearest-even
                                   // _mm256_cvttps_epi32 truncates, matching C rounding; maybe good for scalar testing
}

或者更安全的版本,使用整数进行输入范围移位:如果需要可移植性(仅 AVX1),您可以在此处轻松避免 FMA,并且也可以使用整数加法作为输出。但是我们知道输出范围足够小,它总是可以精确地表示任何整数

static inline
__m256i range_scale_safe_fma(__m256i data, uint32_t a, uint32_t b, uint32_t c, uint32_t d)
{
     // avoid rounding errors when computing the scale factor, but convert double->float on the final result
    const __m256 scale = _mm256_set1_ps((d - c) / (double)(b - a));
    const __m256 cvec = _m256_set1_ps(c);

    __m256i in_offset = _mm256_add_epi32(data, _mm256_set1_epi32(-a));  // add can more easily fold a load of a memory operand than sub because it's commutative.  Only some compilers will do this for you.
    __m256  in_fp = _mm256_cvtepi32_ps(in_offset);
    __m256  out = _mm256_fmadd_ps(in_fp, scale, _mm256_set1_ps(c));  // in*scale + c
    return _mm256_cvtps_epi32(out);
}

没有 FMA,您仍然可以使用 vmulps。如果您这样做,最好在添加 c 之前转换回整数,尽管 vaddps 会是安全的。

你可以在循环中使用它

void foo(uint32_t *arr, ptrdiff_t len)
{
    if (len < 24) special case;

    alignas(32) uint32_t tmpbuf[16];

    // peel half of first iteration for software pipelining / loop rotation
    __m256i arrdata = _mm256_loadu_si256((const __m256i*)&arr[0]);
    __m256i outrange = range_scale_safe_fma(arrdata);
    _mm256_store_si256((__m256i*)tmpbuf, outrange);

    // could have used an unsigned loop counter
    // since we probably just need an if() special case handler anyway for small len which could give len-23 < 0
    for (ptrdiff_t i = 0 ; i < len-(15+8) ; i+=16 ) {

        // prep next 8 elements
        arrdata = _mm256_loadu_si256((const __m256i*)&arr[i+8]);
        outrange = range_scale_safe_fma(arrdata);
        _mm256_store_si256((__m256i*)&tmpbuf[8], outrange);

        // use first 8 elements
        for (int j=0 ; j<8 ; j++) {
            use tmpbuf[j] which corresponds to arr[i+j]
        }

        // prep 8 more for next iteration
        arrdata = _mm256_loadu_si256((const __m256i*)&arr[i+16]);
        outrange = range_scale_safe_fma(arrdata);
        _mm256_store_si256((__m256i*)&tmpbuf[0], outrange);

        // use 2nd 8 elements
        for (int j=8 ; j<16 ; j++) {
            use tmpbuf[j] which corresponds to arr[i+j]
        }
    }

    // use tmpbuf[0..7]
    // then cleanup: one vector at a time until < 8 or < 4 with 128-bit vectors, then scalar
}

这些变量名听起来很愚蠢,但我想不出更好的办法。

此软件流水线是一种优化;你可以让它工作/尝试一次使用一个向量。 (如果需要,可以使用 _mm_cvtsi128_si32(_mm256_castsi256_si128(outrange)) 优化从重新加载到 vmovd 的第一个元素的重新加载。)


特殊情况

如果您知道(b - a) 是2 的幂,您可以使用tzcntbsf 进行位扫描,然后相乘。 (它们有内在函数,比如 GNU C __builtin_ctz() 来计算尾随零。)

或者您能否确保(b - a) 始终是 2 的幂?

或者更好,如果 (b - a) / (d - c) 是 2 的精确幂,那么整个事情就可以是 sub/right shift/add。

如果您不能始终确保有时您仍然需要一般情况,但也许可以有效地做到这一点。

【讨论】:

  • 谢谢。包含该子任务的 MR 最终被合并。如果没有您的帮助,实施的效果会大大降低。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2023-03-03
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多