使用fast-math时，为什么GCC或Clang不会将倒数优化为一条指令？

Question

使用fast-math时，为什么GCC或Clang不会将倒数优化为一条指令？

c++ssecompiler-optimizationsimdfast-math

17

有没有人知道为什么使用fast-math选项时，GCC/Clang不会将下面代码示例中的函数test1优化为只使用RCPPS指令？是否有另一个编译器标志可以生成这个代码？

typedef float float4 __attribute__((vector_size(16)));

float4 test1(float4 v)
{
    return 1.0f / v;
}

您可以在此处查看编译后的输出：https://goo.gl/jXsqat

- Chris_F

2个回答

5

我正在尝试在我的一个应用程序中进行浮点数计算密集型热路径的实验，并发现了类似的情况。我通常不会查看编译器生成的指令，所以我有些惊讶并深入研究了数学细节。

这是由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，因此eps <= 1.5 * 2^-12，所以：

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

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

那么这只是一些神奇的指令序列，碰巧让我们获得了额外的精度吗？并不完全是这样。它基于牛顿迭代法（我了解到在与汇编工作密切相关的人中被称为牛顿-拉夫逊方法）。

牛顿迭代法是一种寻找根的方法。给定一些函数f（x），它通过从近似解x_0开始，并逐步改进来查找f（x）= 0的近似解。牛顿迭代公式如下：

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

在我们的情况下，我们可以将找到倒数 1/a 转化为找到函数 f(x) = a*x - 1 的根，其中导数为 f'(x) = a。将其代入牛顿迭代的方程中，我们得到：

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

两个观察结果：

In this case the Newton iteration actually gives us an exact result, rather than just a better approximation. This makes sense, because Newton's method works by making a linear approximation of f around x_n. In this case f is completely linear, so the approximation is perfect. However...
Computing the Newton iteration requires us to divide by a, which is the exact computation we're trying to approximate. This creates a circular problem. We break the cycle by modifying the Newton iteration to use our approximate reciprocal x_n for the division by 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

- JSQuareD

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Peter Cordes · Accepted Answer

因为RCPPS的精度远低于float除法，所以启用这种优化选项不适合作为-ffast-math的一部分。 gcc手册的x86目标选项中确实有一个选项（使用牛顿-拉弗森迭代），可以让gcc在使用-ffast-math时使用它们。SIMD和标量的每个指令基本上具有相同的性能，而牛顿迭代数学是相同的：SSE/AVX快速矢量rsqrt和倒数，视精度而定 / 使用SSE2的牛顿-拉弗森方法 - 能否有人解释这三行代码。

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

请注意，英特尔的新款Skylake设计进一步提高了FP除法性能，将延迟降至8-11个时钟周期，吞吐量为1/3个时钟周期（对于256位向量，吞吐量为每5个时钟周期一个，但vdivps的延迟相同）。他们扩大了除数，因此AVX vdivps ymm现在与128位向量相同的延迟。

（从Snb到Haswell使用256b的div和sqrt，延迟/倒数吞吐量大约是128b的两倍，因此它们显然只有128b宽的除法器。）Skylake还将两个操作都进行了更多的流水线处理，因此大约可以同时进行4个除法操作。sqrt也更快。

在未来几年中，一旦Skylake普及，只有需要多次除以相同的东西时才值得使用rcpps。rcpps和一些fma可能具有稍高的吞吐量但延迟更差。此外，vdivps只是单个uop；因此，在除法执行期间将有更多的执行资源可用于同时发生的其他操作。

AVX512的初始实现如何尚不确定。假设如果FP除法性能是瓶颈，则通过Newton-Raphson迭代进行rcpps和一对FMA可能会获胜。如果uop吞吐量是瓶颈，并且有大量其他工作可以在除法执行期间完成，则vdivps zmm仍然很好（当然，除数重复使用时除外）。浮点除法与浮点乘法有关FP吞吐量与延迟的更多信息。