将整数范围映射到另一个范围
Map integer range onto another range
在运行时，我有 2 个范围由它们的 uint32_t 边界 a..b 和 c..d 定义。第一个范围往往比第二个范围大得多：8 < (b - a) / (d - c) < 64.
确切限制：a >= 0、b <= 2^31 - 1、c >= 0、d <= 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] 中的多个整数的情况下，可以 return 这些整数中的任何一个。
听起来像是一个简单的比率问题，并且已经在许多问题中得到了解答，例如
Convert a number range to another range, maintaining ratio
但在这里我需要一个非常非常快速的解决方案。
此例程是专门排序算法的关键部分，将对已排序数组的每个元素至少调用一次。
如果不降低整体性能，SIMD 解决方案也是可以接受的。
实际的运行时除法（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 CPUs 具有 2/clock 负载吞吐量，因此将 8x uint32_t 返回到整数寄存器的最佳选择可能是矢量存储/整数重新加载，而不是每次花费 2 微指令ALU 洗牌元素。这有一些延迟，所以我建议在循环通过该标量之前转换为可能为 16 或 32 个整数（64 或 128 字节）的 tmp 缓冲区，即 2x 或 4x __m256i。
或者交替转换和存储一个向量，然后循环遍历您之前转换的 另一个 的 8 个元素。即软件流水线。乱序执行可以隐藏延迟，但无论您对内存做什么，您都已经在扩展其延迟隐藏功能以防止缓存未命中。
根据您的 CPU（例如 Haswell 或某些 Skylake），使用 256 位矢量指令可能会限制您的最大涡轮增压比其他方式略低。您可能会考虑一次只处理 4 个向量，但您会在每个元素上花费更多的微指令。
如果不是 SIMD，那么即使是标量 C++ fma() 仍然很好，对于 vfmadd213sd，但是使用内部函数是从 float -> int 进行舍入（而不是截断）的一种非常方便的方法（vcvtps2dq 而不是 vcvttps2dq）。
请注意，uint32_t <-> float 转换在 AVX512 之前不能直接使用。对于标量，您只需将 to/from int64_t 转换为无符号低半部分的截断/零扩展。
非常方便（如评论中所述）您的输入是有范围限制的，因此如果您将它们解释为带符号的整数，它们具有相同的值（带符号的非负数）。已知 x 和 x-a（以及 b-a）都是正数并且 <= INT32_MAX 即 0x7FFFFFFF。 （或者至少是非负数。零也可以。）
浮点数舍入
对于 SIMD，单精度 float 非常 有利于 SIMD 吞吐量。高效打包转换 to/from 签名 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结合起来，那么
- 在要添加的值中，我们可能会因灾难性取消而导致精度损失。
- 我们需要对原始输入值执行
(float)x。如果 a 很大而 b-a 很小，即接近可能输入范围顶部的小范围，舍入误差可以将所有可能的 x 值映射到同一个浮点数。
- 为了充分利用 FMA，我们希望在转换回整数之前执行
+c，如果 d-c 是一个小输出范围但 [=51]，这将再次冒输出舍入错误的风险=] 是巨大的。在你的情况下没问题；使用 d <= 2^20 - 1 我们知道 float 可以准确地表示 c..d 范围内的每个输出整数值。
如果您没有输入范围限制，您可以在缩放前通过在输入上使用整数 (x-a)+0x80000000U 和在输出上使用 ...+c+0x80000000U 来对范围偏移 to/from 进行签名（之后四舍五入到最接近的 int32_t)。但这会为小的 uint32_t 输入（接近 0）引入巨大的 float 舍入误差，这些输入的范围会偏移到接近 INT_MIN.
我们不需要对 b-a 或 d-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 的幂，您可以使用 tzcnt 或 bsf 进行位扫描，然后相乘。 （这些有内在函数，比如 GNU C __builtin_ctz() 来计算尾随零。）
或者你能确保 (b - a) 总是 2 的幂吗？
或者更好的是，如果 (b - a) / (d - c) 是 2 的精确幂，那么整个事情就可以减/右移/加。
如果你不能总是确保你有时仍然需要一般情况，但也许可以有效地做到这一点。
在运行时，我有 2 个范围由它们的 uint32_t 边界 a..b 和 c..d 定义。第一个范围往往比第二个范围大得多：8 < (b - a) / (d - c) < 64.
确切限制：a >= 0、b <= 2^31 - 1、c >= 0、d <= 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] 中的多个整数的情况下，可以 return 这些整数中的任何一个。
听起来像是一个简单的比率问题，并且已经在许多问题中得到了解答，例如
Convert a number range to another range, maintaining ratio
但在这里我需要一个非常非常快速的解决方案。
此例程是专门排序算法的关键部分，将对已排序数组的每个元素至少调用一次。
如果不降低整体性能，SIMD 解决方案也是可以接受的。
实际的运行时除法（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 CPUs 具有 2/clock 负载吞吐量，因此将 8x uint32_t 返回到整数寄存器的最佳选择可能是矢量存储/整数重新加载，而不是每次花费 2 微指令ALU 洗牌元素。这有一些延迟，所以我建议在循环通过该标量之前转换为可能为 16 或 32 个整数（64 或 128 字节）的 tmp 缓冲区，即 2x 或 4x __m256i。
或者交替转换和存储一个向量，然后循环遍历您之前转换的 另一个 的 8 个元素。即软件流水线。乱序执行可以隐藏延迟，但无论您对内存做什么，您都已经在扩展其延迟隐藏功能以防止缓存未命中。
根据您的 CPU（例如 Haswell 或某些 Skylake），使用 256 位矢量指令可能会限制您的最大涡轮增压比其他方式略低。您可能会考虑一次只处理 4 个向量，但您会在每个元素上花费更多的微指令。
如果不是 SIMD，那么即使是标量 C++ fma() 仍然很好，对于 vfmadd213sd，但是使用内部函数是从 float -> int 进行舍入（而不是截断）的一种非常方便的方法（vcvtps2dq 而不是 vcvttps2dq）。
请注意，uint32_t <-> float 转换在 AVX512 之前不能直接使用。对于标量，您只需将 to/from int64_t 转换为无符号低半部分的截断/零扩展。
非常方便（如评论中所述）您的输入是有范围限制的，因此如果您将它们解释为带符号的整数，它们具有相同的值（带符号的非负数）。已知 x 和 x-a（以及 b-a）都是正数并且 <= INT32_MAX 即 0x7FFFFFFF。 （或者至少是非负数。零也可以。）
浮点数舍入
对于 SIMD，单精度 float 非常 有利于 SIMD 吞吐量。高效打包转换 to/from 签名 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结合起来，那么
- 在要添加的值中，我们可能会因灾难性取消而导致精度损失。
- 我们需要对原始输入值执行
(float)x。如果a很大而b-a很小，即接近可能输入范围顶部的小范围，舍入误差可以将所有可能的x值映射到同一个浮点数。 - 为了充分利用 FMA，我们希望在转换回整数之前执行
+c，如果d-c是一个小输出范围但 [=51]，这将再次冒输出舍入错误的风险=] 是巨大的。在你的情况下没问题；使用d<= 2^20 - 1 我们知道float可以准确地表示 c..d 范围内的每个输出整数值。
如果您没有输入范围限制，您可以在缩放前通过在输入上使用整数 (x-a)+0x80000000U 和在输出上使用 ...+c+0x80000000U 来对范围偏移 to/from 进行签名（之后四舍五入到最接近的 int32_t)。但这会为小的 uint32_t 输入（接近 0）引入巨大的 float 舍入误差，这些输入的范围会偏移到接近 INT_MIN.
我们不需要对 b-a 或 d-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 的幂，您可以使用 tzcnt 或 bsf 进行位扫描，然后相乘。 （这些有内在函数，比如 GNU C __builtin_ctz() 来计算尾随零。）
或者你能确保 (b - a) 总是 2 的幂吗？
或者更好的是，如果 (b - a) / (d - c) 是 2 的精确幂，那么整个事情就可以减/右移/加。
如果你不能总是确保你有时仍然需要一般情况，但也许可以有效地做到这一点。