将整数范围映射到另一个范围

实际的运行时除法(FP和integer)非常慢，所以您肯定希望避免这种情况。您编写该表达式的方式可能会编译为包含除法，因为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/时钟负载吞吐量，因此将8xuint32_t重新加载到整数寄存器中的最佳选择可能是矢量存储/整数重载，而不是为ALU混洗的东西每个元素花费2个uops。这有一些延迟，所以我建议在循环通过标量之前，转换成一个可能为16或32 int(64或128字节)的tmp缓冲区，即2x或4x__m256i。

或者，可以交替转换和存储一个向量，然后在的8个元素上循环，另一个的元素是您之前转换的。即软件流水线。无序执行可以隐藏延迟，但无论您对内存做什么，都会扩展其针对缓存未命中的延迟隐藏功能。

根据您的CPU(例如Haswell或一些Skylake)，使用256位矢量指令可能会使您的最大turbo比其他情况下的略低。你可能会考虑一次只做4个向量，但你会在每个元素上花费更多的uop。

如果不是SIMD，那么对于vfmadd213sd，即使是标量C++fma()仍然是好的，但使用内部函数是从float->int(vcvtps2dq而不是vcvttps2dq)获得舍入(而不是截断)的一种非常方便的方法。

注意uint32_t<->float转换直到AVX512才直接可用。对于标量，您只需将int64_t转换为无符号下半部分的截断/零扩展即可。

非常方便的是(如注释中所述)您的输入是范围有限的，因此如果您将它们解释为有符号整数，则它们具有相同的值(有符号非负)。x和x-a(以及b-a)都已知是正的并且&ltINT32_MAX，即0x7FFFFFFF。(或者至少是非负的。零是好的。)

浮动舍入

对于SIMD，单精度float对SIMD吞吐量非常有利。到有符号int32_t的高效压缩转换。但并不是每个CCD_ 35都可以精确地表示为CCD_。越大的值四舍五入到最接近的偶数，最接近的2^2，2^3的倍数，或更大的2^24以上的值。

使用SIMDdouble是可能的，但需要一些混洗。

我不认为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<=2^²⁰-1我们知道float可以精确地表示该c.d范围内的每个输出整数值

如果您没有输入范围约束，您可以通过在输入上使用整数(x-a)+0x80000000U，在输出上使用...+c+0x80000000U(四舍五入到最接近的int32_t后)，将范围移位到缩放前的带符号。但对于小的uint32_t输入(接近0)，这将引入巨大的float舍入误差，这些输入的范围将移动到接近INT_MIN。

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

示例：

CCD_ 63向量应该在该内联之后由编译器从循环中提升出来，或者你可以手动完成。

这需要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__builtin_ctz()，用来计算尾随零。)

或者您能确保(b - a)总是2的幂吗？

或者更好的是，如果(b - a) / (d - c)是2的精确幂，那么整个事情可以是次/右移/加法

如果你不能总是确保你有时仍然需要一般情况，但也许可以有效地做到这一点。