快速开平方取倒数的算法--嵌入式ARM转载

#include<stdio.h>
#include<string.h>
#include <stdlib.h>     /* atof */


/*
计算=1/sqrt(n)
*/
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
#ifndef Q3_VM
#ifdef __linux__
    assert( !isnan(y) );
    // bk010122 - FPE?
#endif
#endif
    return y;
}


int main(int argc, char const *argv[])
{
	float f9 = 81.0f;
	f9 = Q_rsqrt(f9);

	printf("f9=%f\n", f9);

	return 0;
}

　　运行结果：

f9=0.111086

和计算机1/sqrt（81）很接近1/9=0.111111

 

相比 sqrt() 函数，这套算法要快将近4倍，要知道，编译器自带的函数，可是经过严格仔细的汇编优化的啊！

牛顿迭代法的原理是先猜测一个值，然后从这个值开始进行叠代。因此，猜测的值越准，叠代的次数越少。卡马克选了0x5f3759df这个值作为猜测的结果，再加上后面的移位算法，得到的y非常接近1/sqrt(n)。这样，我们只需要2次牛顿迭代法就可以达到我们所需要的精度。

 

函数返回1/sqrt(x),这个函数在图像处理中比sqrt(x)更有用。

 

注意到这个正数只用了一次叠代!(其实就是根本没用叠代，直接运算)。编译、实验，这个团数不仅工作的很好，而且比标准的sqrt()函数快4倍！

这个简洁的定数，最核心，也是最让人费解的，就是标注了what the fuck的一句 i   = 0x5f3759df - ( i >> 1 )；再加上y   = y * ( threehalfs - ( x2 * y * y ) )。

 

两句话就完成了开方运算！而且注意到，核心那句是移位运算，速度极快！特别在很多没有乘法指令的RISC结构CPU上，这样做是极其高效的。

 

算法的原理就是使用牛顿迭代法,用 x-f(x)/f'(x) 来不断的逼近 f(x)=a 的根。

 

求平方根：f(x)=x^2=a ,f'(x)= 2*x, f(x)/f'(x)=x/2,把 f(x) 代入 x-f(x)/f'(x)后有(x+a/x)/2，

 

现在我们选 a=5,选一个猜测值比如 2，  那么我们可以这么算  5/2 = 2.5; (2.5+2)/2 = 2.25; 5/2.25 = ……  这样反复迭代下去，结果必定收敛于 sqrt(5)。

 

但是卡马克作者真正厉害的地方是他选择了一个神秘的常数 0x5f375a86来计算那个梦“值，

就是我们加注释的那一行那行算出的值非常接近1/sqrt(n)这样我们只需要2次牛顿迭代就可以达到我们所需要的精度。