为什么 GCC 或 Clang 在使用快速数学时不优化 1 条指令的倒数
Why does GCC or Clang not optimise reciprocal to 1 instruction when using fast-math
有谁知道为什么 GCC/Clang 不会在下面的代码示例中使用乐观函数 test1 在使用 fast-math 选项时只使用 RCPPS 指令？是否有另一个编译器标志可以生成此代码？
typedef float float4 __attribute__((vector_size(16)));
float4 test1(float4 v)
{
return 1.0f / v;
}
你可以在这里看到编译后的输出：https://goo.gl/jXsqat
因为RCPPS的精度比float低很多 division.
启用该优化的选项不适合作为 -ffast-math 的一部分。
x86 target options of the gcc manual says there in fact is an option that (with -ffast-math) does get gcc to use them (with a Newton-Raphson iteration - / Newton Raphson with SSE2 - can someone explain me these 3 lines - SIMD 和标量每条指令的性能基本相同，牛顿迭代数学相同):
-mrecip
This option enables use of RCPSS and RSQRTSS instructions (and their vectorized variants RCPPS and RSQRTPS) with an additional
Newton-Raphson step to increase precision instead of DIVSS and SQRTSS
(and their vectorized variants) for single-precision floating-point
arguments. These instructions are generated only when
-funsafe-math-optimizations is enabled together with -finite-math-only and -fno-trapping-math. Note that while the throughput of the sequence
is higher than the throughput of the non-reciprocal instruction, the
precision of the sequence can be decreased by up to 2 ulp (i.e. the
inverse of 1.0 equals 0.99999994).
Note that GCC implements 1.0f/sqrtf(x) in terms of RSQRTSS (or RSQRTPS) already with -ffast-math (or the above option combination),
and doesn't need -mrecip.
Also note that GCC emits the above sequence with additional Newton-Raphson step for vectorized single-float division and
vectorized sqrtf(x) already with -ffast-math (or the above option
combination), and doesn't need -mrecip.
-mrecip=opt
This option controls which reciprocal estimate instructions may be used. opt is a comma-separated list of options, which may be preceded
by a ‘!’ to invert the option:
’all’
Enable all estimate instructions.
‘default’
Enable the default instructions, equivalent to -mrecip.
‘none’
Disable all estimate instructions, equivalent to -mno-recip.
‘div’
Enable the approximation for scalar division.
‘vec-div’
Enable the approximation for vectorized division.
‘sqrt’
Enable the approximation for scalar square root.
‘vec-sqrt’
Enable the approximation for vectorized square root.
So, for example, -mrecip=all,!sqrt enables all of the reciprocal approximations, except for square root.
请注意，Intel 的新 Skylake 设计 further improves FP division performance，延迟为 8-11c，吞吐量为 1/3c。 （或者 256b 向量每 5c 吞吐量一个，但 vdivps 的延迟相同）。他们 widened dividers，因此 AVX vdivps ymm 现在与 128b 向量的延迟相同。
(Haswell 的 SnB 做了 256b div 和 sqrt，延迟/接收吞吐量大约是两倍，所以他们显然只有 128b-wide dividers.) Skylake 还对这两个操作进行了更多的流水线处理，因此大约有 4 div 个操作可以进行。 sqrt 也更快。
所以几年后，一旦 Skylake 得到ide传播，只有当你需要 divide 时，才值得做 rcpps事情多次。 rcpps 和一对 fma 的吞吐量可能略高但延迟更差。此外，vdivps 只是一个 uop；因此更多的执行资源将可用于与 division.
同时发生的事情
AVX512 的初始实现会是什么样子还有待观察。据推测，如果 FP division 性能成为瓶颈，rcpps 和几个用于 Newton-Raphson 迭代的 FMA 将是一个胜利。如果 uop 吞吐量是一个瓶颈，并且在 divisions 运行期间还有很多其他工作要做，那么 vdivps zmm 可能仍然很好（除非重复使用相同的 divisor，当然）。
Floating point division vs floating point multiplication 详细介绍了 FP 吞吐量与延迟。
我在我的一个应用程序中试验了一个浮点数学密集型热路径，并发现了类似的东西。我通常不看我的编译器发出的指令，所以我有点惊讶并深入研究了数学细节。
这是 gcc 生成的指令集，由我用携带的计算注释：
test1(float __vector): ; xmm0 = a
rcpps xmm1, xmm0 ; xmm1 = 1 / xmm0 = 1/a
mulps xmm0, xmm1 ; xmm0 = xmm0 * xmm1 = a * 1/a
mulps xmm0, xmm1 ; xmm0 = xmm0 * xmm1 = a * (1/a)^2
addps xmm1, xmm1 ; xmm1 = xmm1 + xmm1 = 2 * (1/a)
subps xmm1, xmm0 ; xmm1 = xmm1 - xmm0 = 2 * (1/a) - a * (1/a)^2
movaps xmm0, xmm1 ; xmm0 = xmm1 = 2 * (1/a) - a * (1/a)^2
ret
这是怎么回事？为什么要浪费额外的 4 条指令来计算数学上只等价于倒数的表达式？
嗯，rcpps 指令只是一个 近似值 倒数。其他算术指令（mulps、addps、subps）精确到单精度。近似倒数函数写成r(x)。然后最后变成y = 2*r(a) - a*r(a)^2。如果我们将 r(x) = (1 + eps) * (1/x) 替换为 eps 是相对误差，我们将得到：
y = 2 * (1 + eps) * (1/a) - a * (1 + eps)^2 * (1/a)^2
= (2 + 2*eps - (1 + eps)^2) * (1/a)
= (2 + 2*eps - (1 + 2*eps + eps^2)) * (1/a)
= (1 - eps^2) * (1/a)
rcppsis less than1.5 * 2^-12的相对误差，所以eps <= 1.5 * 2^-12，所以：
eps^2 <= 2.25 * 2^-24
< 1.5 * 2^-23
因此，通过执行这些额外的指令，我们将精度从 12 位提高到 23 位。请注意，单精度浮点数具有 24 位精度，因此我们 almost 在这里获得全精度。
那么这只是一些神奇的指令序列，恰好让我们获得了额外的精度吗？不完全的。它基于 Newton's method（据我所知，经常使用汇编的人将其称为 Newton-Raphson）。
牛顿法是一种求根法。给定一些函数 f(x)，它通过从近似解 x_0 开始并迭代改进它来找到 f(x) = 0 的近似解。牛顿迭代由下式给出：
x_n+1 = x_n - f(x_n) / f'(x_n)
在我们的例子中，我们可以将寻找 a 的倒数 1/a 重新表述为寻找函数 f(x) = a*x - 1 的根，以及导数 f'(x) = a。将其代入牛顿迭代方程，我们得到：
x_n+1 = x_n - (a*x_n - 1) / a
两个观察：
在这种情况下，牛顿迭代实际上给了我们一个精确的结果，而不仅仅是一个更好的近似值。这是有道理的，因为牛顿法的工作原理是在 x_n 附近对 f 进行线性近似。在这种情况下 f 是完全线性的，因此近似是完美的。然而...
计算牛顿迭代需要我们除以 a，这是我们试图近似的精确计算。这就产生了一个循环问题。我们通过修改牛顿迭代以使用我们的近似倒数 x_n 来除以 a:
来打破循环
x_n+1 = x_n - (a*x_n - 1) * x_n
~= x_n - (a*x_n - 1) / a
这个迭代会工作得很好，但从矢量数学的角度来看并不好：它需要减去 1。要使用向量数学来做到这一点，需要准备一个包含 1 序列的向量寄存器。这需要额外的指令和额外的寄存器。
我们可以重写迭代来避免这种情况：
x_n+1 = x_n - (a*x_n - 1) * x_n
= x_n - (a*x_n^2 - x_n)
= 2*x_n - a*x_n^2
现在替换x_0 = r(a)，我们从上面恢复我们的表达式：
y = x_1 = 2*r(a) - a*r(a)^2
有谁知道为什么 GCC/Clang 不会在下面的代码示例中使用乐观函数 test1 在使用 fast-math 选项时只使用 RCPPS 指令？是否有另一个编译器标志可以生成此代码？
typedef float float4 __attribute__((vector_size(16)));
float4 test1(float4 v)
{
return 1.0f / v;
}
你可以在这里看到编译后的输出：https://goo.gl/jXsqat
因为RCPPS的精度比float低很多 division.
启用该优化的选项不适合作为 -ffast-math 的一部分。
x86 target options of the gcc manual says there in fact is an option that (with -ffast-math) does get gcc to use them (with a Newton-Raphson iteration -
-mrecipThis option enables use of RCPSS and RSQRTSS instructions (and their vectorized variants RCPPS and RSQRTPS) with an additional Newton-Raphson step to increase precision instead of DIVSS and SQRTSS (and their vectorized variants) for single-precision floating-point arguments. These instructions are generated only when -funsafe-math-optimizations is enabled together with -finite-math-only and -fno-trapping-math. Note that while the throughput of the sequence is higher than the throughput of the non-reciprocal instruction, the precision of the sequence can be decreased by up to 2 ulp (i.e. the inverse of 1.0 equals 0.99999994).Note that GCC implements 1.0f/sqrtf(x) in terms of RSQRTSS (or RSQRTPS) already with -ffast-math (or the above option combination), and doesn't need -mrecip.
Also note that GCC emits the above sequence with additional Newton-Raphson step for vectorized single-float division and vectorized sqrtf(x) already with -ffast-math (or the above option combination), and doesn't need -mrecip.
-mrecip=optThis option controls which reciprocal estimate instructions may be used. opt is a comma-separated list of options, which may be preceded by a ‘!’ to invert the option:
’all’ Enable all estimate instructions. ‘default’ Enable the default instructions, equivalent to -mrecip. ‘none’ Disable all estimate instructions, equivalent to -mno-recip. ‘div’ Enable the approximation for scalar division. ‘vec-div’ Enable the approximation for vectorized division. ‘sqrt’ Enable the approximation for scalar square root. ‘vec-sqrt’ Enable the approximation for vectorized square root.So, for example, -mrecip=all,!sqrt enables all of the reciprocal approximations, except for square root.
请注意，Intel 的新 Skylake 设计 further improves FP division performance，延迟为 8-11c，吞吐量为 1/3c。 （或者 256b 向量每 5c 吞吐量一个，但 vdivps 的延迟相同）。他们 widened dividers，因此 AVX vdivps ymm 现在与 128b 向量的延迟相同。
(Haswell 的 SnB 做了 256b div 和 sqrt，延迟/接收吞吐量大约是两倍，所以他们显然只有 128b-wide dividers.) Skylake 还对这两个操作进行了更多的流水线处理，因此大约有 4 div 个操作可以进行。 sqrt 也更快。
所以几年后，一旦 Skylake 得到ide传播，只有当你需要 divide 时，才值得做 rcpps事情多次。 rcpps 和一对 fma 的吞吐量可能略高但延迟更差。此外，vdivps 只是一个 uop；因此更多的执行资源将可用于与 division.
AVX512 的初始实现会是什么样子还有待观察。据推测，如果 FP division 性能成为瓶颈，rcpps 和几个用于 Newton-Raphson 迭代的 FMA 将是一个胜利。如果 uop 吞吐量是一个瓶颈，并且在 divisions 运行期间还有很多其他工作要做，那么 vdivps zmm 可能仍然很好（除非重复使用相同的 divisor，当然）。
Floating point division vs floating point multiplication 详细介绍了 FP 吞吐量与延迟。
我在我的一个应用程序中试验了一个浮点数学密集型热路径，并发现了类似的东西。我通常不看我的编译器发出的指令，所以我有点惊讶并深入研究了数学细节。
这是 gcc 生成的指令集，由我用携带的计算注释：
test1(float __vector): ; xmm0 = a
rcpps xmm1, xmm0 ; xmm1 = 1 / xmm0 = 1/a
mulps xmm0, xmm1 ; xmm0 = xmm0 * xmm1 = a * 1/a
mulps xmm0, xmm1 ; xmm0 = xmm0 * xmm1 = a * (1/a)^2
addps xmm1, xmm1 ; xmm1 = xmm1 + xmm1 = 2 * (1/a)
subps xmm1, xmm0 ; xmm1 = xmm1 - xmm0 = 2 * (1/a) - a * (1/a)^2
movaps xmm0, xmm1 ; xmm0 = xmm1 = 2 * (1/a) - a * (1/a)^2
ret
这是怎么回事？为什么要浪费额外的 4 条指令来计算数学上只等价于倒数的表达式？
嗯，rcpps 指令只是一个 近似值 倒数。其他算术指令（mulps、addps、subps）精确到单精度。近似倒数函数写成r(x)。然后最后变成y = 2*r(a) - a*r(a)^2。如果我们将 r(x) = (1 + eps) * (1/x) 替换为 eps 是相对误差，我们将得到：
y = 2 * (1 + eps) * (1/a) - a * (1 + eps)^2 * (1/a)^2
= (2 + 2*eps - (1 + eps)^2) * (1/a)
= (2 + 2*eps - (1 + 2*eps + eps^2)) * (1/a)
= (1 - eps^2) * (1/a)
rcppsis less than1.5 * 2^-12的相对误差，所以eps <= 1.5 * 2^-12，所以：
eps^2 <= 2.25 * 2^-24
< 1.5 * 2^-23
因此，通过执行这些额外的指令，我们将精度从 12 位提高到 23 位。请注意，单精度浮点数具有 24 位精度，因此我们 almost 在这里获得全精度。
那么这只是一些神奇的指令序列，恰好让我们获得了额外的精度吗？不完全的。它基于 Newton's method（据我所知，经常使用汇编的人将其称为 Newton-Raphson）。
牛顿法是一种求根法。给定一些函数 f(x)，它通过从近似解 x_0 开始并迭代改进它来找到 f(x) = 0 的近似解。牛顿迭代由下式给出：
x_n+1 = x_n - f(x_n) / f'(x_n)
在我们的例子中，我们可以将寻找 a 的倒数 1/a 重新表述为寻找函数 f(x) = a*x - 1 的根，以及导数 f'(x) = a。将其代入牛顿迭代方程，我们得到：
x_n+1 = x_n - (a*x_n - 1) / a
两个观察：
在这种情况下，牛顿迭代实际上给了我们一个精确的结果，而不仅仅是一个更好的近似值。这是有道理的，因为牛顿法的工作原理是在
x_n附近对f进行线性近似。在这种情况下f是完全线性的，因此近似是完美的。然而...计算牛顿迭代需要我们除以
来打破循环a，这是我们试图近似的精确计算。这就产生了一个循环问题。我们通过修改牛顿迭代以使用我们的近似倒数x_n来除以a:x_n+1 = x_n - (a*x_n - 1) * x_n ~= x_n - (a*x_n - 1) / a
这个迭代会工作得很好，但从矢量数学的角度来看并不好：它需要减去 1。要使用向量数学来做到这一点，需要准备一个包含 1 序列的向量寄存器。这需要额外的指令和额外的寄存器。
我们可以重写迭代来避免这种情况：
x_n+1 = x_n - (a*x_n - 1) * x_n
= x_n - (a*x_n^2 - x_n)
= 2*x_n - a*x_n^2
现在替换x_0 = r(a)，我们从上面恢复我们的表达式：
y = x_1 = 2*r(a) - a*r(a)^2