110

我一直在分析我们在英特尔酷睿双核上的一些核心数学，在研究平方根的各种方法时，我注意到了一些奇怪的事情：使用 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 古怪的 Newton-Raphson 近似技术，它比硬件运行得更好，为 4.3ns/float，尽管误差为 2 10分之 1 （这对我的目的来说太多了）。

当我尝试使用 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 的内置硬件平方根操作码比从其他两个数学运算中合成它要慢？

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

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

（编辑：stephentyrone 正确地指出，对长字符串的操作应该使用矢量化 SIMD 打包操作，例如rsqrtps——但这里的数组数据结构仅用于测试目的：我真正想要测量的是代码中使用的标量性能不能向量化。）

4

6 回答 6

222

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

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

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

于 2009-10-06T23:52:07.290 回答
8

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

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

于 2011-07-12T14:32:14.897 回答
6

而不是提供一个答案，这实际上可能是不正确的（我也不会检查或争论缓存和其他东西，假设它们是相同的）我会尝试向您指出可以回答您问题的来源。
区别可能在于 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 正在使用一些查找表，它应该只在结果时使用不需要很准确，虽然 - 我也可能错了，就像前一段时间一样:)。

于 2009-10-06T23:55:34.477 回答
6

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

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

以下是共识出错的地方：

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

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

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

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

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

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

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，即使您不计算平方根，它们（再次）通常也很有用。

于 2019-12-05T01:08:40.727 回答
4

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

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

然后增量是：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这也很容易。

于 2012-08-02T22:20:35.850 回答
-3

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

于 2016-07-05T14:17:12.570 回答