为啥 SSE scalar sqrt(x) 比 rsqrt(x) * x 慢？

Posted 2023-02-19

【中文标题】为啥 SSE scalar sqrt(x) 比 rsqrt(x) * x 慢？

【英文标题】：Why is SSE scalar sqrt(x) slower than rsqrt(x) * x?为什么 SSE scalar sqrt(x) 比 rsqrt(x) * x 慢？

【发布时间】：2010-12-04 11:08:36

【问题描述】：

我一直在分析我们在英特尔酷睿双核上的一些核心数学，在研究平方根的各种方法时，我发现了一些奇怪的东西：使用 SSE 标量运算，取倒数平方根更快并乘以得到 sqrt，而不是使用本机 sqrt 操作码！

我正在使用类似以下的循环对其进行测试：

inline float TestSqrtFunction( float in );

void TestFunc()

  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );

我已经为 TestSqrtFunction 尝试了几个不同的主体，但有些时间确实让我摸不着头脑。到目前为止，最糟糕的是使用本机 sqrt() 函数并让“智能”编译器“优化”。在 24ns/float 时，使用 x87 FPU 这太糟糕了：

inline float TestSqrtFunction( float in )
  return sqrt(in); 

接下来我尝试使用内部函数强制编译器使用 SSE 的标量 sqrt 操作码：

inline void SSESqrt( float * restrict pOut, float * restrict pIn )

   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss

这更好，为 11.9ns/float。我还尝试了Carmack's wacky Newton-Raphson approximation technique，它在 4.3ns/float 上比硬件运行得更好，尽管误差为 1 in 2¹⁰（这对我的目的来说太多了）。

当我尝试对 reciprocal 平方根进行 SSE 运算，然后使用乘法得到平方根 ( x * 1/√x = √x ) 时，这真是太棒了。尽管这需要两个相关操作，但它是迄今为止最快的解决方案，1.24ns/float 并且精确到 2^-14：

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )

   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss

我的问题基本上是什么给了？ 为什么 SSE 的硬件内置平方根操作码比从其他两个数学运算中合成它要慢？

我确定这确实是操作本身的成本，因为我已经验证：

所有数据都适合缓存，并且 访问是顺序的 函数是内联的 展开循环没有区别 编译器标志设置为完全优化（我检查过，程序集很好）

(edit：stephentyrone 正确指出对长字符串的操作应该使用矢量化 SIMD 打包操作，例如 rsqrtps — 但这里的数组数据结构仅用于测试目的：什么我真的想衡量的是 scalar 在无法矢量化的代码中使用的性能。）

【问题讨论】：

x / sqrt(x) = sqrt(x)。或者，换一种说法：x^1 * x^(-1/2) = x^(1 - 1/2) = x^(1/2) = sqrt(x)

当然，inline float SSESqrt( float restrict fIn ) float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; 。但这是一个坏主意，因为如果 CPU 将浮点数写入堆栈然后立即将它们读回，它很容易导致加载命中存储停顿——特别是从向量寄存器到浮点寄存器以获取返回值是个坏消息。此外，SSE 内在函数表示的底层机器操作码无论如何都采用地址操作数。

LHS 的重要性取决于给定 x86 的特定 gen 和步进：我的经验是，在 i7 之前的任何设备上，在寄存器集之间移动数据（例如 FPU 到 SSE 到 eax）是非常糟糕，而 xmm0 和堆栈之间的往返不是，因为英特尔的存储转发。您可以自己确定时间来确定。通常，查看潜在 LHS 的最简单方法是查看发出的程序集并查看数据在寄存器组之间的位置；您的编译器可能会做聪明的事，也可能不会。至于标准化向量，我在这里写下了我的结果：bit.ly/9W5zoU

对于 PowerPC，是的：IBM 有一个 CPU 模拟器，可以通过静态分析预测 LHS 和许多其他管道气泡。一些 PPC 还具有您可以轮询的 LHS 硬件计数器。 x86 更难；好的分析工具越来越少（VTune 这些天有些坏了），重新排序的管道不太确定。您可以尝试通过测量每个周期的指令来凭经验测量它，这可以通过硬件性能计数器精确完成。可以使用例如 PAPI 或 PerfSuite (bit.ly/an6cMt) 读取“指令退休”和“总周期”寄存器。

你也可以简单地在一个函数上写一些排列，然后给它们计时，看看是否有任何特别受停顿的影响。英特尔没有公布关于他们的管道工作方式的许多细节（他们的 LHS 完全是一个肮脏的秘密），所以我学到的很多东西都是通过查看导致其他拱门停滞的场景（例如 PPC )，然后构建一个受控实验，看看 x86 是否也有。

【参考方案1】：

sqrtss 给出正确舍入的结果。 rsqrtss 给出倒数的近似值，精确到大约 11 位。

sqrtss 正在生成更准确的结果，用于需要准确度的情况。 rsqrtss 存在于近似值足够但需要速度的情况。如果您阅读 Intel 的文档，您还会发现一个指令序列（倒数平方根近似，后跟一个 Newton-Raphson 步骤），它提供了几乎完全的精度（大约 23 位的精度，如果我没记错的话），并且仍然有点比sqrtss 快。

编辑：如果速度很关键，并且您确实在循环中为许多值调用它，您应该使用这些指令的矢量化版本，rsqrtps 或 sqrtps，每条指令都处理四个浮点数。

【讨论】：

n/r 步骤为您提供 22 位的准确度（它加倍）； 23 位将完全准确。

@Jasper Bekkers：不，不会。首先，float 具有 24 位精度。其次，sqrtss 被正确舍入，在舍入之前需要约 50 位，并且无法使用单精度的简单 N/R 迭代来实现。

肯定是这个原因。为了扩展这个结果：英特尔的 Embree 项目 (software.intel.com/en-us/articles/…) 将矢量化用于其数学运算。您可以在该链接下载源代码并查看他们如何制作 3/4 D 矢量。他们的向量归一化使用 rsqrt 后跟 newton-raphson 的迭代，这样就非常准确并且仍然比 1/ssqrt 快！

一个小警告：如果 x 为零或无穷大，xrsqrt(x) 将导致 NaN。 0*rsqrt(0) = 0 * INF = NaN。 INFrsqrt(INF) = INF * 0 = NaN。出于这个原因，NVIDIA GPU 上的 CUDA 将近似的单精度平方根计算为 recip(rsqrt(x))，而硬件提供了对倒数和倒数平方根的快速逼近。显然，处理这两种特殊情况的显式检查也是可能的（但在 GPU 上会更慢）。

@BrandonPelfrey 您在哪个文件中找到了 Newton Rhapson 步骤？

【参考方案2】：

除法也是如此。 MULSS(a,RCPSS(b)) 比 DIVSS(a,b) 快得多。事实上，即使您通过 Newton-Raphson 迭代提高其精度，它仍然更快。

英特尔和 AMD 都在其优化手册中推荐了这种技术。在不需要符合 IEEE-754 的应用程序中，使用 div/sqrt 的唯一原因是代码可读性。

【讨论】：

Broadwell 和更高版本具有更好的 FP 除法性能，因此像 clang 这样的编译器选择不在最近的 CPU 上使用倒数 + 牛顿作为标量，因为它通常不更快。在大多数循环中，div 并不是唯一的操作，因此即使存在 divps 或 divss，总的 uop 吞吐量通常也是瓶颈。请参阅Floating point division vs floating point multiplication，我的答案中有一节说明为什么rcpps 不再是吞吐量胜利。 （或延迟获胜），以及划分吞吐量/延迟的数字。

如果您的精度要求太低以至于您可以跳过牛顿迭代，那么a * rcpss(b) 可以更快，但它仍然比a/b 更多微秒！

【参考方案3】：

实际上可能是不正确的，而不是提供答案（我也不会检查或争论缓存和其他东西，假设它们是相同的）我会尝试将您指向可以回答的来源你的问题。 区别可能在于 sqrt 和 rsqrt 的计算方式。你可以在这里阅读更多内容http://www.intel.com/products/processor/manuals/。我建议从阅读您正在使用的处理器功能开始，有一些信息，尤其是关于 rsqrt 的信息（cpu 正在使用具有巨大近似值的内部查找表，这使得获得结果变得更加简单）。看起来，rsqrt 比 sqrt 快得多，1 个额外的 mul 操作（成本不高）可能不会改变这里的情况。

编辑：可能值得一提的事实很少： 1. 一旦我对我的图形库进行了一些微优化，我就使用 rsqrt 来计算向量的长度。 （而不是 sqrt，我将平方和乘以它的 rsqrt，这正是您在测试中所做的），并且它表现更好。 2. 使用简单的查找表计算 rsqrt 可能更容易，对于 rsqrt，当 x 趋于无穷大时，1/sqrt(x) 趋于 0，因此对于小的 x 函数值不会改变（很多），而对于sqrt - 它趋于无穷大，所以就是这么简单的情况;)。

另外，澄清一下：我不确定我在我链接的书中的哪里找到它，但我很确定我已经读过 rsqrt 正在使用一些查找表，它应该只使用，当结果不需要精确时，虽然 - 我也可能错了，就像前一段时间一样:)。

【参考方案4】：

几年前已经有许多其他答案了。以下是共识正确的地方：

rsqrt* 指令计算平方根倒数的近似值，大约 11-12 位。 它是用尾数索引的查找表（即 ROM）实现的。 （实际上，它是一个压缩查找表，类似于旧的数学表，使用调整低位来节省晶体管。） 之所以可用，是因为它是 FPU 用于“真实”平方根算法的初始估计值。 还有一个近似的倒数指令，rcp。这两条指令都是 FPU 如何实现平方根和除法的线索。

以下是共识的错误：

SSE 时代的 FPU 不使用 Newton-Raphson 来计算平方根。这在软件中是一种很好的方法，但在硬件中以这种方式实现是错误的。

计算平方根倒数的 N-R 算法有这个更新步骤，正如其他人所指出的：

x' = 0.5 * x * (3 - n*x*x);

这是很多数据相关的乘法和一个减法。

以下是现代 FPU 实际使用的算法。

给定b[0] = n，假设我们可以找到一系列数字Y[i]使得b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2接近1。然后考虑：

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

显然x[n] 接近sqrt(n) 和y[n] 接近1/sqrt(n)。

我们可以对倒数平方根使用 Newton-Raphson 更新步骤来得到一个好的Y[i]：

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

然后：

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

和：

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

下一个关键观察是b[i] = x[i-1] * y[i-1]。所以：

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

然后：

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

也就是说，给定初始 x 和 y，我们可以使用以下更新步骤：

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

或者，更高级的，我们可以设置h = 0.5 * y。这是初始化：

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

这是更新步骤：

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

这是 Goldschmidt 的算法，如果你在硬件中实现它，它有一个巨大的优势：“内循环”是三个乘加，没有别的，其中两个是独立的，可以流水线化。

在 1999 年，FPU 已经需要一个流水线加/减电路和一个流水线乘法电路，否则 SSE 不会很“流”。 1999 年，每个电路只需要一个，就可以以完全流水线的方式实现这个内部循环，而不会浪费大量的硬件来计算平方根。

今天，当然，我们已经向程序员公开了乘加法。同样，内部循环是三个流水线 FMA，即使您不计算平方根，它们（再次）通常也很有用。

【讨论】：

相关：How sqrt() of GCC works after compiled? Which method of root is used? Newton-Raphson? 有一些硬件 div/sqrt 执行单元设计的链接。 Fast vectorized rsqrt and reciprocal with SSE/AVX depending on precision - 软件中的一次牛顿迭代，带或不带 FMA，与 _mm256_rsqrt_ps 一起使用，带有 Haswell 性能分析。如果您在循环中没有其他工作并且会严重限制分频器吞吐量，通常只有一个好主意。硬件 sqrt 是单 uop，因此可以与其他工作混合使用。

【参考方案5】：

Newton-Raphson 使用等于-f/f' 的增量收敛到f(x) 的零，其中f' 是导数。

对于x=sqrt(y)，可以尝试使用f(x) = x^2 - y解决f(x) = 0对于x；

那么增量为：dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x 其中有一个缓慢的划分。

您可以尝试其他功能（例如f(x) = 1/y - 1/x^2），但它们会同样复杂。

现在让我们看看1/sqrt(y)。您可以尝试f(x) = x^2 - 1/y，但同样复杂：例如dx = 2xy / (y*x^2 - 1)。 f(x) 的一个不明显的替代选择是：f(x) = y - 1/x^2

然后：dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

啊！这不是一个微不足道的表达式，但你只有乘法，没有除法。 => 更快！

并且：完整的更新步骤new_x = x + dx 则为：

x *= 3/2 - y/2 * x * x 这也很简单。

【参考方案6】：

因为这些指令忽略舍入模式，并且不处理浮点异常或非规范化数字，所以速度更快。由于这些原因，流水线、推测和乱序执行其他 fp 指令要容易得多。

【讨论】：

显然错了。 FMA 取决于当前的舍入模式，但在 Haswell 及更高版本上每个时钟的吞吐量为两个。凭借两个全流水线 FMA 装置，Haswell 可以同时运行多达 10 个 FMA。正确答案是rsqrt 的准确率多 低，这意味着在查找表格以获得初步猜测之后要做的工作要少得多（或根本没有？）。