为什么SSE标量sqrt(x)比rsqrt(x)*x慢？

Question

我正在对我们核心数学代码进行性能分析，使用Intel Core Duo处理器。在尝试使用不同的计算平方根的方法时，我注意到了一个奇怪的现象：使用SSE标量运算，先取倒数平方根再乘以它比使用本机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()函数并让“智能”编译器进行优化。使用x87 FPU每秒24纳秒/浮点数，这个速度非常差：

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

接下来我尝试的是使用内置函数来强制编译器使用SSE的标量平方根指令：

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。我还试过卡马克疯狂的牛顿拉弗森逼近技术，速度比硬件还要快，只需要4.3ns/float，但误差为2的10次方之一（对于我的目的来说太大了）。

最棒的是，当我使用SSE操作进行倒数平方根运算，然后使用乘法得到平方根（x * 1/√x = √x）。虽然这需要两个相依的操作，但它是迄今为止速度最快的解决方案，每个浮点数只需1.24ns，且准确精度可达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内置的硬件平方根操作码比通过合成其他两个数学运算慢？

我确定这主要是由于操作本身的成本，因为我已经验证了：

• 所有数据都适合缓存，并且访问是顺序的
• 函数被内联
• 展开循环没有任何影响
• 编译器标志设置为完全优化（汇编代码很好，我检查过）

(编辑：stephentyrone正确指出，长数字字符串上的操作应使用矢量化SIMD打包操作，如rsqrtps，但此处的数组数据结构仅用于测试目的：我真正尝试衡量的是无法进行矢量化的代码中的标量性能。)

- Crashworks

---

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

当然，inline float SSESqrt(float restrict fIn)的意思是对输入参数进行平方根计算，并返回结果。但这种实现方式不好，因为如果CPU将浮点数写入堆栈，然后立即读取它们，就很容易引起加载阻塞存储。特别是将向量寄存器中的值转移到浮点寄存器以返回值，则更加不可取。此外，SSE intrinsic所表示的底层机器指令也需要使用地址操作数。 - Crashworks

LHS的重要性取决于具体的x86处理器的型号和版本：我的经验是，在i7之前的处理器上，将数据移动到不同寄存器集（例如从FPU移到SSE再到eax）非常糟糕，而在xmm0和堆栈之间来回传递则没有问题，因为Intel的存储转发技术。您可以自行测试以确保其准确性。通常，查看汇编代码并查看数据在哪些寄存器集之间操作是发现潜在LHS的最简单方法；您的编译器可能会做出聪明的选择，也可能不会。至于向量归一化，我的结果在这里：http://bit.ly/9W5zoU。 - Crashworks

对于PowerPC，是的：IBM有一个CPU模拟器，可以通过静态分析预测LHS和许多其他流水线气泡。一些PPC还具有用于轮询LHS的硬件计数器。对于x86来说更难；良好的分析工具较少（VTune现在有点问题），而且重新排序的流水线不太确定。您可以尝试通过测量每个周期的指令数来从经验上衡量它，这可以通过硬件性能计数器精确完成。可以使用PAPI或PerfSuite（http://bit.ly/an6cMt）读取“已退休指令”和“总周期”寄存器。 - Crashworks

你也可以简单地在一个函数上写几个排列，并计时，看是否有任何一个特别受到停顿的影响。英特尔并没有公布关于他们流水线工作方式的许多细节（他们左手边的一切都是一种肮脏的秘密），所以我学到的很多东西都是通过观察其他架构（例如 PPC）上引起停顿的情况，然后构建一个受控实验来查看 x86 是否也有这个问题。 - Crashworks

Answer 1

sqrtss可以给出正确舍入的结果。rsqrtss提供了一个约11位精度的倒数估计值。

sqrtss在需要精确度时提供更为准确的结果。rsqrtss用于当近似值足够且需要加快速度的情况。如果您查阅英特尔的文档，您将发现一种指令序列（倒数平方根估计值加上单个牛顿-拉弗森步骤），它可以提供接近全精度（如果我没记错，大约有23位精度），并且仍然比sqrtss稍快。

编辑：如果速度非常关键，并且您真的要在许多值的循环中调用此函数，那么您应该使用这些指令的矢量化版本rsqrtps或sqrtps，每个指令可以处理四个浮点数。

- Stephen Canon

---

n/r 步骤可以提供22位精度（将精度增加一倍）；23位精度即为完全精度。 - Jasper Bekkers

不会。首先，浮点数具有24位精度。其次，"sqrtss"是正确舍入的，这需要在舍入之前约50位，并且不能使用单精度的简单N/R迭代实现。 - Stephen Canon

这绝对是原因。为了扩展这个结果：英特尔的Embree项目（http://software.intel.com/en-us/articles/embree-photo-realistic-ray-tracing-kernels/）在其数学计算中使用矢量化。您可以在该链接上下载源代码并查看他们如何处理3/4D向量。他们的矢量归一化使用rsqrt，然后进行牛顿-拉弗森迭代，这样非常准确，而且仍然比1/ssqrt更快！ - Brandon Pelfrey

小小的注意事项：如果x为零或无穷大，则xrsqrt(x)的结果为NaN。0rsqrt(0) = 0 * INF = NaN. INF*rsqrt(INF) = INF * 0 = NaN。因此，NVIDIA GPU上的CUDA使用recip(rsqrt(x))计算近似单精度平方根，硬件提供快速逆和倒数平方根的近似值。显然，可以使用明确检查处理这两种特殊情况（但在GPU上速度会较慢）。 - njuffa

@BrandonPelfrey 你在哪个文件中找到了牛顿-拉普森步骤？ - fredoverflow

Answer 2

几年前已经有关于此问题的其他答案了。以下是达成一致的正确观点：

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

以下是达成一致的错误观点：

- SSE时代的FPU不使用牛顿-拉夫逊方法来计算平方根。在软件中，这是一种很好的方法，但在硬件中实现会是一个错误。 - 计算倒数平方根的N-R算法具有这样的更新步骤，正如其他人所指出的：

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

这需要进行很多依赖数据的乘法和一次减法。

接下来是现代浮点运算器实际使用的算法。

假设给定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)。

我们可以使用牛顿-拉弗森反平方根更新步骤来获得一个好的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年，浮点数单元已经需要一个流水线加/减电路和一个流水线乘法电路，否则SSE将不会非常适用于数据流。只需要每个电路中的一个就可以在1999年以全流水线方式实现这个内部循环，而不会浪费大量硬件资源在平方根上。

当然，今天我们还有融合乘加（FMA）暴露给程序员。同样，内部循环是三个流水线FMA，即使你不计算平方根，它们（再次）通常也很有用。

- Pseudonym

---

相关：GCC编译后如何运行sqrt()函数？使用哪种根号算法？牛顿迭代法？中有一些关于硬件div/sqrt执行单元设计的链接。使用SSE/AVX进行快速向量化的rsqrt和倒数，具体取决于精度——在软件中进行一个牛顿迭代，可以带或不带FMA，用于与_mm256_rsqrt_ps配合使用，带有Haswell性能分析。通常只有当您在循环中没有其他工作并且会在除法器吞吐量上遇到瓶颈时才是一个好主意。硬件sqrt是单个uop，因此可以与其他工作混合使用。 - Peter Cordes

Answer 3

除法也是如此。MULSS（a，RCPSS（b））比DIVSS（a，b）快得多。事实上，即使通过牛顿-拉夫森迭代增加其精度，它仍然更快。

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

- Spat

---

Broadwell及其后代具有更好的浮点除法性能，因此像clang这样的编译器选择不在最近的CPU上为标量使用倒数+牛顿方法，因为通常情况下它并不比其他方法更快。在大多数循环中，“div”不是唯一的操作，因此即使存在“divps”或“divss”，总的微操作吞吐量通常仍然是瓶颈。请参见浮点除法vs浮点乘法，我的答案中有一个关于为什么“rcpps”不再是吞吐量优势（或延迟优势）的部分，并附有除法吞吐量/延迟的数字。 - Peter Cordes

如果您的精度要求很低，可以跳过牛顿迭代，那么a * rcpss(b)可能会更快，但它仍然比a/b使用更多的uops！ - Peter Cordes

Answer 4

我不会提供一个可能是不正确的答案（我也不会检查或争论缓存和其他东西，假设它们是相同的），而是尝试指向可以回答你问题的来源。
差异可能在于如何计算sqrt和rsqrt。你可以在这里阅读更多信息http://www.intel.com/products/processor/manuals/。我建议从阅读你正在使用的处理器函数开始，有一些信息，特别是关于rsqrt（CPU使用内部查找表进行大量近似，使得获得结果变得简单得多）。看起来，rsqrt比sqrt快得多，以至于1个额外的乘法操作（成本并不高）可能不会改变这里的情况。

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

另外，澄清一下：我不确定我在链接的书籍中找到了它的位置，但我很确定我读过rsqrt使用了一些查找表，并且只有在结果不需要精确时才应该使用它，尽管-我也可能错了，因为那是一段时间前:)。

- Marcin Deptuła

Answer 5

牛顿-拉夫逊法使用增量等于-f/f'，其中f'为导数，收敛于f(x)的零点。

对于x=sqrt(y)，您可以尝试使用f(x) = x^2 - y来解决f(x) = 0以获得x，然后增量为：dx = -f/f' = 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，这也很容易。

- skal

Answer 6

它更快，因为这些指令忽略了舍入模式，并且不处理浮点异常或非规格化数。出于这些原因，它更容易流水线化、推测和执行其他 fp 指令。

- Witek

---

显然是错误的。FMA 取决于当前的舍入模式，但在 Haswell 及更高版本上每个时钟周期有两个吞吐量。Haswell 有两个完全流水线化的 FMA 单元，可以同时处理多达 10 个 FMAs。正确答案是 rsqrt 的精度要低得多，这意味着在查找表以获取起始猜测后需要做的工作要少得多（或者根本不需要？）。 - Peter Cordes