寻找更快的平方根倒数算法

机器学习的模型相关计算中，有很多诡异的运算。单个运算的开销很不起眼，但如果这些运算的量足够大，也会对性能产生一定的影响。这里谈的就是一个简单的运算：

a = b / sqrt(c);

对于 C/C++ 语言的程序员来说，sqrt 已经是非常基础的库函数，它的底层实现也仅仅是简单的一句 FSQRT (双精度是 SQRTSD) 指令，看起来没有什么优化的余地。但事实上 intel 提供了一个更快的指令，那就是 SQRTSS，利用这条指令，平方根倒数的计算速度可以达到 sqrt 版本的两倍（实测，与[1]相同）。你可以这样使用它：

#include <xmmintrin.h>
...
__m128 in = _mm_load_ss(&c);
__m128 out = _mm_sqrt_ss(in);
_mm_store_ss(&c, out);
a = b/c;
...

但这就是优化的尽头了么？不，单就求平方根倒数来说，还有一个神奇的近似算法，叫做 Fast Inverse Square Root（平方根倒数速算法）。一个神人在 Quake III Arena 游戏中使用了一个神奇的数字 0x5f3759df，创造了这个神奇的算法，这个算法可以将平方根倒数的计算速度提升到 sqrt 的 3 倍多（实测，效果比[1]差）。

float Q_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
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
 
        return y;
}

但 3 倍就是优化的尽头了么？很不幸，邪恶的 Intel 提供了这样一条指令 RSQRTSS，从硬件上支持了这个近似算法。利用这条指令，平方根倒数的计算速度能够达到 sqrt 版本的 6 倍以上！！！

#include <xmmintrin.h>
...
    __m128 in = _mm_load_ss(&c);
    __m128 out = _mm_rsqrt_ss(in);
    _mm_store_ss(&c, out);
    a = b*c;
...

虽然平方根倒数速算法只是一种近似算法，并且只有单精度版本，但是对 RSQRTSS 指令的简单测试发现大部分情况下误差在万分之一以下，指令说明书中给出的误差是 ±1.5*2^-12[2]，在非精确数值计算的工程系统中已经足够用了。

它带来的一个更有趣的后果是：如果使用 RSQRTSS 计算出来 c 的平方根倒数，然后再乘以 c，就得到了 c 的平方根近似值。用它可以反过来加速 sqrt 的运算！[1]

注1：编译相关程序时，需要打开优化开关，以实现函数的 inline
注2：RSQRTSS 和 SQRTSS 均有一个向量版本，如 RSQRTPS，可以同时计算 4 个 float 的平方根倒数；

[1] Timing square root
[2] RSQRTSS