为什么GCC或Clang在使用快速数学时不能优化1指令的倒数

Question

有谁知道为什么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

Answer 1

由于精度RCPPS是很多低于float分工.

启用该优化的选项不适合作为其中一部分-ffast-math.

gcc手册的x86目标选项说实际上有一个选项(with -ffast-math)确实让gcc使用它们(使用Newton-Raphson迭代):

• -mrecip 此选项允许使用RCPSS和RSQRTSS指令(及其矢量化变量RCPPS和RSQRTPS)以及额外的Newton-Raphson步骤来提高精度,而不是DIVSS和SQRTSS(及其矢量化变量)用于单精度浮点参数.仅当-funsafe-math-optimizations与-finite-math-only和-fno-trapping-math一起启用时,才会生成这些指令.请注意,虽然序列的吞吐量高于非互易指令的吞吐量,但序列的精度可降低最多2 ulp(即1.0的倒数等于0.99999994).

  注意,GCC在已经使用-ffast-math(或上述选项组合)的RSQRTSS(或RSQRTPS)方面实现了1.0f/sqrtf(x),并且不需要-mrecip.

  另请注意,GCC使用额外的Newton-Raphson步骤为向量化单浮点除法和向量化sqrtf(x)已经使用-ffast-math(或上述选项组合)发出上述序列,并且不需要-mrecip.

• -mrecip=opt

该选项控制可以使用哪些倒数估计指令.opt是以逗号分隔的选项列表,可以在前面加上'!' 反转选项:

’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. 

因此,例如,-mrecip = all,!sqrt启用除平方根之外的所有倒数近似.

请注意,英特尔新的Skylake设计进一步提高了FP部门的性能,延迟为8-11c,吞吐量为1/3c.(或者对于256b矢量,每5c吞吐量一个,但相同的延迟vdivps).它们扩大了分频器,因此AVX vdivps ymm现在与128b矢量的延迟相同.

(SnB到Haswell的256b div和sqrt的延迟/接收吞吐量大约是其两倍,所以他们显然只有128b的分频器.)Skylake还管理了两个操作,因此大约有4个div操作可以在飞行中.sqrt也更快.

所以在几年之后,一旦Skylake普及,只有rcpps你需要多次划分相同的东西才值得做. rcpps并且一对fma可能具有略高的吞吐量但延迟更差.而且,vdivps只有一个uop; 因此,更多的执行资源可用于与分部同时发生的事情.

AVX512的初始实现方式还有待观察.rcpps如果FP部门的表现成为瓶颈,那么推测和Newton-Raphson迭代的几个FMA将是一场胜利.如果uop吞吐量是一个瓶颈,并且在划分飞行时还有很多其他工作要做,那么vdivps zmm可能仍然很好(当然,除非重复使用相同的除数).

Answer 2

我正在我的一个应用程序中试验浮点数学密集型热路径，并发现了类似的东西。我通常不会查看编译器发出的指令，所以我有点惊讶并深入研究了数学细节。

这是 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)

的相对误差rcpps 小于 1.5 * 2^-12, so eps <= 1.5 * 2^-12, so：

eps^2 <= 2.25 * 2^-24
      <  1.5  * 2^-23

因此，通过执行这些额外的指令，我们将精度从 12 位提高到了 23 位。请注意，单精度浮点数有 24 位精度，因此我们在这里几乎获得全精度。

那么，这只是一些神奇的指令序列，恰好可以让我们获得额外的精度吗？不完全的。它基于牛顿方法（我认为经常使用汇编的人将其称为牛顿拉夫森方法）。

牛顿法是一种求根法。给定某个函数，它通过从近似解开始并迭代改进它来f(x)找到 的近似解。牛顿迭代由下式给出：f(x) = 0x_0

x_n+1 = x_n - f(x_n) / f'(x_n)

在我们的例子中，我们可以将求1/a的倒数重新表述a为求函数 的根f(x) = a*x - 1，导数f'(x) = a。将其代入牛顿迭代方程，我们得到：

x_n+1 = x_n - (a*x_n - 1) / a

两个观察结果：

1. 在这种情况下，牛顿迭代实际上给了我们一个精确的结果，而不仅仅是一个更好的近似值。f这是有道理的，因为牛顿方法的工作原理是对周围进行线性近似x_n。在这种情况下f是完全线性的，因此近似是完美的。然而...

2. 计算牛顿迭代需要我们除以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