实际运行时除法(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 中所讨论的)您的输入是有范围限制的,因此如果您将它们解释为有符号整数,它们具有相同的值(有符号非负数)。已知x 和x-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 结合起来,那么
- 在要添加的值中,我们可能会因灾难性取消而损失精度。
- 我们需要对原始输入值执行
(float)x。如果a 很大而b-a 很小,即接近可能输入范围顶部的小范围,则舍入误差可以将所有可能的x 值映射到相同的浮点数。
- 为了充分利用 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-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 的精确幂,那么整个事情就可以是 sub/right shift/add。
如果您不能始终确保有时您仍然需要一般情况,但也许可以有效地做到这一点。