Cra*_*rks 106 floating-point performance x86 assembly sse
我一直在英特尔Core Duo上分析我们的一些核心数学,并且在研究平方根的各种方法时,我注意到一些奇怪的事情:使用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 /浮动.我也试过卡马克的古怪牛顿迭代逼近技术,其运行时间比硬件,为4.3ns /浮甚至更好,虽然以1比2的错误10(这是太多了,我的目的).
doozy是当我尝试SSE op的倒数平方根时,然后使用乘法得到平方根(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- 但这里的数组数据结构仅用于测试目的:我真正想要测量的是用于代码的标量性能无法矢量化的.)
Ste*_*non 211
sqrtss给出了正确的舍入结果. rsqrtss给出了倒数的近似值,精确到大约11位.
sqrtss当需要准确性时,产生更准确的结果. rsqrtss对于近似值足够的情况存在,但需要速度.如果您阅读英特尔的文档,您还会发现一个指令序列(倒数平方根逼近,后跟单个Newton-Raphson步骤),它几乎可以提供完全精度(如果我没记错的话,精度约为23位),并且仍然有点快于sqrtss.
编辑:如果速度很关键,并且你真的在循环中为许多值调用它,那么你应该使用这些指令的向量化版本,rsqrtps或者sqrtps两者都为每条指令处理四个浮点数.
分裂也是如此.MULSS(a,RCPSS(b))比DIVSS(a,b)快.实际上,即使使用Newton-Raphson迭代提高其精度,它仍然更快.
英特尔和AMD都在其优化手册中推荐了这种技术.在不需要符合IEEE-754标准的应用程序中,使用div/sqrt的唯一原因是代码可读性.
而不是提供答案,实际上可能是不正确的(我也不会检查或争论缓存和其他东西,让我们说它们是相同的)我会试着指出你可以回答你的问题的来源.
差异可能在于如何计算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正在使用一些查找表,并且只应该使用它,当结果不需要确切,虽然 - 我可能也错了,因为它是在不久前:).
几年前就已经有许多其他答案了。以下是共识的正确之处：
以下是共识出错的地方：
正如其他人所指出的，用于计算倒数平方根的 NR 算法具有此更新步骤：
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，即使您不计算平方根，它们（再次）通常也很有用。