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

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

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

《寻找更快的平方根倒数算法》上有2条评论

发表回复

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