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

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

#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次牛顿迭代就可以达到我们所需要的精度。