从著名的Q_rsqrt说起

凡是做过一段时间开发的人大多都知道算法优化的必要性。兴许是近期某社区的火爆,导致了很多人都会把卡马克在开发quark3时做的那段可谓是天神级别的“Q_rsqrt”作为一个经典的算法优化教程来说。然后就是某位仁兄就用了这段代码替换掉了原本math.h调用中的sqrtf函数,堂而皇之的提交进了代码库。

话说这个“Q_rsqrt”是对牛顿展开法求平方根“猜数”过程的一种特例,它的优势在于只要做一次递归就可以简单的得到一个数平方根倒数的近似值。
为了验证这个问题,我写了一段代码,主要是计算从1到0xfffffff所有数字的平方根的倒数,并计算这个过程的耗时。为了这个测试,又写了两段代码,分别是通过标准的math.h这种“笨办法”以及采用SSE指令直接出结果这中强依赖于系统平台的方法。
兴许作为程序员,你会很抵制SSE这种“新技术”,或者指手画脚它的兼容性。但这个指令集已经覆盖到了所有10年之内的x86主机上了,不管你是amd还是Intel,这个指令集已经是标配,几乎碰不上兼容性问题。arm平台就另当别论了。
#include <stdio.h>
#include <time.h>
//#define TEST_Q
//#define TEST_MATH 1
#define TEST_SSE 1


#ifdef TEST_Q
float new_rsqrt( float number ) //这是卡马克大神级别的优化
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;      // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration

    return y;
}
#endif

#ifdef TEST_MATH
#include <math.h>

float new_rsqrt(float f){ //这是调用标准的math.h函数库

    return 1/sqrtf(f);
}
#endif

#ifdef TEST_SSE
#include <xmmintrin.h>
float new_rsqrt(float f){ //这是调用用CPU SSE指令集中rsqrt函数直接得出结果

    __m128 m_a = _mm_set_ps1(f);
    __m128 m_b = _mm_rsqrt_ps(m_a);

    return m_b[0];
}
#endif


int main(int argc, const char * argv[]) {

    int i ;
    float r;
    printf("%f\n", new_rsqrt(9)); //因为结果应该是1/3,无限循环,用以对比精度
    clock_t start_time = clock();
    for (i=0;i< 0xfffffff; i++){
        r = new_rsqrt(i);
    }

    printf("Time cost: %ld ns\n", clock() - start_time);

}

 得出的结果如下:
LitrindeMacBook-Pro:~ litrin$ time ./test_q
0.332953
Time cost: 7746636 ns

real 0m7.773s
user 0m7.733s
sys 0m0.019s
LitrindeMacBook-Pro:~ litrin$ time ./test_sse
0.333252
Time cost: 1391658 ns

real 0m1.399s
user 0m1.391s
sys 0m0.005s
LitrindeMacBook-Pro:~ litrin$ time ./test_normal
0.333333
Time cost: 3882738 ns

real 0m3.892s
user 0m3.879s
sys 0m0.009s
  •  Q_rsqrt:耗时7.7秒,精度3位
  • sse: 耗时1.4秒,精度4位
  • math库:耗时3.9秒,精度6位

数据差异很明显了,那个经典的算法优化实际上是大大的拖累了整体的性能,别着急,我们可以从汇编的角度上看看这是为什么。
首先是Q_rsqrt,足足30步计算,尽管都是很简单的四则运算范畴,但步数超多。

Playground`new_rsqrt:
    0x100000e20 <+0>:   pushq  %rbp
    0x100000e21 <+1>:   movq   %rsp, %rbp
    0x100000e24 <+4>:   movss  0x150(%rip), %xmm1
    0x100000e2c <+12>:  movl   $0x5f3759df, %eax
    0x100000e31 <+17>:  movl   %eax, %ecx
    0x100000e33 <+19>:  movss  0x145(%rip), %xmm2
    0x100000e3b <+27>:  movss  %xmm0, -0x4(%rbp)
    0x100000e40 <+32>:  movss  %xmm1, -0x1c(%rbp)
    0x100000e45 <+37>:  mulss  -0x4(%rbp), %xmm2
    0x100000e4a <+42>:  movss  %xmm2, -0x14(%rbp)
    0x100000e4f <+47>:  movss  -0x4(%rbp), %xmm0
    0x100000e54 <+52>:  movss  %xmm0, -0x18(%rbp)
    0x100000e59 <+57>:  movq   -0x18(%rbp), %rdx
    0x100000e5d <+61>:  movq   %rdx, -0x10(%rbp)
    0x100000e61 <+65>:  movq   -0x10(%rbp), %rdx
    0x100000e65 <+69>:  sarq   $0x1, %rdx
    0x100000e69 <+73>:  subq   %rdx, %rcx
    0x100000e6c <+76>:  movq   %rcx, -0x10(%rbp)
->  0x100000e70 <+80>:  movss  -0x10(%rbp), %xmm0
    0x100000e75 <+85>:  movss  %xmm0, -0x18(%rbp)
    0x100000e7a <+90>:  movss  -0x18(%rbp), %xmm0
    0x100000e7f <+95>:  movss  -0x14(%rbp), %xmm2
    0x100000e84 <+100>: mulss  -0x18(%rbp), %xmm2
    0x100000e89 <+105>: mulss  -0x18(%rbp), %xmm2
    0x100000e8e <+110>: subss  %xmm2, %xmm1
    0x100000e92 <+114>: mulss  %xmm1, %xmm0
    0x100000e96 <+118>: movss  %xmm0, -0x18(%rbp)
    0x100000e9b <+123>: movss  -0x18(%rbp), %xmm0
    0x100000ea0 <+128>: popq   %rbp
    0x100000ea1 <+129>: retq

接下来是math库,10步,一个sqrtss命令得出结果,再divss除一下。

Playground`new_rsqrt:
    0x100000e80 <+0>:  pushq  %rbp
    0x100000e81 <+1>:  movq   %rsp, %rbp
    0x100000e84 <+4>:  movss  0xf0(%rip), %xmm1
    0x100000e8c <+12>: movss  %xmm0, -0x4(%rbp)
    0x100000e91 <+17>: movss  -0x4(%rbp), %xmm0
    0x100000e96 <+22>: sqrtss %xmm0, %xmm0
->  0x100000e9a <+26>: divss  %xmm0, %xmm1
    0x100000e9e <+30>: movaps %xmm1, %xmm0
    0x100000ea1 <+33>: popq   %rbp
    0x100000ea2 <+34>: retq

轮到SSE了,13步,初始化sse寄存器用了太多的步数,但一个rsqrtps命令直接得出结果,然后从寄存器中取出数据。从另一个方面将,128位的SSE操作实际上是可以支持在同一个时钟周期内对4个浮点数做rsqrt的,但由于基础测试代码的倾向性,我们只能发挥出1/4的SSE性能。

Playground`new_rsqrt:
    0x100000e80 <+0>:  pushq  %rbp
    0x100000e81 <+1>:  movq   %rsp, %rbp
    0x100000e84 <+4>:  movss  %xmm0, -0x38(%rbp)
->  0x100000e89 <+9>:  movss  %xmm0, -0x34(%rbp)
    0x100000e8e <+14>: movss  %xmm0, -0x14(%rbp)
    0x100000e93 <+19>: shufps $0x0, %xmm0, %xmm0
    0x100000e97 <+23>: movaps %xmm0, -0x30(%rbp)
    0x100000e9b <+27>: movaps %xmm0, -0x50(%rbp)
    0x100000e9f <+31>: movaps %xmm0, -0x10(%rbp)
    0x100000ea3 <+35>: rsqrtps %xmm0, %xmm0
    0x100000ea6 <+38>: movaps %xmm0, -0x60(%rbp)
    0x100000eaa <+42>: popq   %rbp
    0x100000eab <+43>: retq

说了这么多并不是要直接否定了Q_rsqrt这种优化的方式,开发这个方法的时代兴许很多PC都没有一个独立的浮点运算单元(x87),这种做法无疑会大大提升平台的整体性能。作为对算法的极致追求,这种做法值得我们去不断雕琢。但事实上随着这些年被摩尔定律推着走,我们已经可以拿出更多的精力去思考一些超出于基础算法的东西了。或者又回到那个所谓程序员的悖论:到底应该花一年半的时间将代码的效率提升到150%,还是索性等一年半让计算机的处理速度提升到200%。

其实在这个实验中,我还计划用AVX指令做一组数据,但AVX只是将同时进行计算的数量提升到了8个,但单次计算时间相比SSE没有太大的提升,如果不修改代码的话,我并不期待它在单个单个顺序操作rsqrt上的性能能大幅领先SSE。感兴趣的话可以通过这段代码实现:

#ifdef TEST_AVX
#include <immintrin.h&rt;

float new_rsqrt( float number )
{
    __m256 a = _mm256_set1_ps(number);
    __m256 b = _mm256_rsqrt_ps(a);
    return b[1];
}
#endif

需要注意的是,avx尽管已经覆盖了5年之内的x86CPU,但确实是存在平台兼容性不佳的问题,一般的编译器都不会默认打开AVX支持,需要增加-mavx编译参数。

如果还有什么要说的,那就是:这几年CPU的进化方向已经到了多核和重指令集上去了,那以后的优化思路就是能用并发就用并发,能用重指令集就用重指令集!不过话说重指令可能在某些情况下使CPU核心的降频从而导致性能的下降出现与本例相悖的结果。这就是另外一个故事了,本文不再累述。

推荐阅读:
事出前些日子有人咨询我:“在某
时延 latency(亦称为延
似乎每次开头都要讲述一下计算机

发表评论

电子邮件地址不会被公开。 必填项已用*标注

请补全下列算式: *

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据