最直接的思路是二分法或者牛顿迭代法,不过新搜索到了更多厉害的算法,能够更快的计算“浮点数开方运算”。
二分法
注意事项是,退出while循环的时候,要用两次的mid值的差去跟eps比较,如果使用low或者up的话,由于精度的问题,可能mid计算完之后仍然是low或者up,这样while退出条件就没改变,可能会陷入死循环;所以这里一个重要的trick是使用两次的mid的值的差去判断是否达到精度。
1 #define eps 1e-9 2 float SqrtByBisection(float n){ 3 float low, up, mid, last; 4 low = 0, up=(n<1? 1:n); 5 mid = (low+up)/2.; 6 7 do{ 8 if(mid*mid>n) 9 up = mid; 10 else 11 low = mid; 12 last = mid; 13 mid = (up+low)/2.; 14 }while( fabsf(mid-last) > eps ); 15 return mid; 16 }
牛顿法
根据y=x2-a曲线,很容易推导出牛顿法的迭代公式为
xi = (xi-1 + a/xi-1) / 2
代码实现细节如下(仍然使用判断两次迭代结果的差值,来判断程序是否结束)
1 float SqrtByNewton(float x){ 2 float val = x; 3 float last; 4 do{ 5 last = val; 6 val = (val + x/val)/2.; 7 }while( fabsf(val-last) > eps ) 8 return val; 9 }
牛顿法虽然比二分法快很多,但是还有改进空间。二分法一般选取原始值作为搜索开始点,牛顿法的朴素写法也是以原始值作为搜索起点的,但是更好的初始值,可以大大减少迭代次数,加速计算。
假如要求x的开方,那么
int temp = (((*(int *)&x)&0xff7fffff)>>1)+(64<<23); float val=*(float*)&temp;
不失为一种较好的估计。具体原理与IEEE 的float的存储机制有关系。
参考博客:https://blog.csdn.net/yutianzuijin/article/details/78839981