在编程开发中，经常需要去计算一个数的平方根，小编一般是直接调用系统函数 Math.Sqrt()。
当然除了系统函数，我们也可以自己编写函数进行求解，常用的方法有二分法和牛顿迭代法。
二分法
二分法的基本思想：对于区间 \([m,n]\) 上连续不断且 \(f(m)·f(n) < 0\) 的函数 \(y=f(x)\)，通过不断地把函数 \(f(x)\) 的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值。
显然，\(\sqrt{a}\) 是函数 \( f(x) =x^2 - a\) 的一个零点，但需要注意的是，如果直接将区间初始化为 \([0,a]\) 是无法保证 \( f(0)·f(a) = a^2 (1-a) < 0 \) 恒成立的。因为当 \( a < 1\) 时，\(f(x)\) 的正实根位于区间 \([0,a]\) 之外，此时需要将区间上限扩大为大于 1 的值。
|
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 |
double SqrtByBisection(double x) //用二分法 { if(x < 0) return Double.NaN; double eps = 1e-5;//允许误差 double low,up; low = 0;//区间下限 if (x < 1) up = 10 * x;//确保区间上限大于 1 else up = x; double mid,last; mid=(low+up)/2; do { if(mid*mid > x) up=mid; else low=mid; last=mid; mid=(up+low)/2; }while(Math.Abs(mid-last) > eps);//误差控制 return mid; } |
牛顿迭代法
牛顿迭代法的基本思想：以曲线的切线与 \(x\) 轴的交点逼近曲线的零点。
具体来说：对于函数 \( f(x)\)，选取一个初始值 \(x_0\)，过点 \((x,f(x_0))\) 作曲线 \( f(x)\) 的切线 \(L_1\)，其方程为 \( y= f(x_0) + f'(x_0)(x-x_0)\)，计算切线 \(L_1\) 与 \(x\) 轴交点的横坐标 \(x_1=x_0 - \frac{f(x_0)}{f'(x_0)}\) 作为 \( f(x) = 0\) 根的近似值，然后过点 \((x,f(x_1))\) 再作曲线 \( f(x)\) 的切线 \(L_2\)，其方程为 \( y= f(x_1) + f'(x_1)(x-x_1)\)，计算切线 \(L_2\) 与 \(x\) 轴交点的横坐标 \(x_2=x_1 - \frac{f(x_1)}{f'(x_1)}\) 作为 \( f(x) = 0\) 根的近似值。重复以上过程，直至满足精度要求。
对于函数 \( f(x) =x^2 - a\)，有 \( f'(x) =2x\)，因此近似根的迭代公式为： \[x_{n+1}=x_{n} - \frac{f(x_{n})}{f'(x_{n})} = =x_{n} - \frac{x^2_{n}-a}{2x_{n}} = \frac{1}{2}(x_n + \frac{a}{x_n})\]
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
double SqrtByNewton(double x) { double eps = 1e-5;//允许误差 double sqrt = x; //初始值 double last; //保存上一个计算的值 do { last = sqrt; sqrt =(sqrt + x/sqrt) / 2; }while(Math.Abs(sqrt-last) > eps);//误差控制 return sqrt; } |
不过从运行效率上来看，似乎并没有自己编写函数的需求，因为系统函数的运行效率明显优于以上两个自编函数。
直到某一天，小编在网络上看到一个神奇的方法，心中万马奔腾，如同代码中注释一样 ：“what the fuck?”
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
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; } |
这段代码来自于 90 年代的经典游戏之一 Quake-III Arena（雷神之锤 3），它的作用是将一个数开平方并取倒数。
有网友测试发现，这段代码竟然比 (float)(1.0/sqrt(x)) 快 4 倍！
毋庸置疑的“神奇”！
至于卡马克（雷神之锤 3 的作者）是怎么发现这个奇妙的数字的，我们不得而知。
不过，普渡大学的数学家 Chris Lomont 看到这个问题后曾有一些研究，他从理论上推导出一个最佳猜测值：0x5f37642f，和卡马克的数字非常接近，但是计算精度和速度不及卡马克的数字。面对这样的而结果，Chris Lomont 很受伤，不过好在通过暴力试算，总算找到了一个比卡马克的数字要好上那么一丁点的数字：0x5f375a86，虽然实际上这两个数字所产生的结果非常近似。
从 Chris Lomont 的分析来看，这个神奇的方法是巧妙利用了单精度浮点数的相关特性。
我们知道，单精度浮点数占用 4 个字节（32 位）存储空间，具体由符号位（1 位），阶码（8 位）和尾数（23 位）构成。以下图为例，符号位为 \(s=0\)，阶码 \(E=124\)，尾数 \(M=2097152\)，对应的数字为:\[x=(-1)^s(1+\frac{M}{2^{23}})2^{E-127} =(-1)^0(1+0.25)2^{-3}\]
鉴于以上原因，在 C# 中实现该神奇算法时，不能像其他算法一样将 float 改写为 double，同时还必须使用 unsafe 关键字。
以下是基于神奇算法的求平方根函数（C#）：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
float InvSqrt(float x) { float xhalf = 0.5f * x; long i; unsafe //C# 不安全代码 { i = *(long*) & x; // 二进制 i = 0x5f3759df - (i >> 1); // 初始值 x = *(float*) & i; // 转化为 float } x = x * (1.5f - xhalf * x * x); // 一次迭代，提高精度 x = x * (1.5f - xhalf * x * x); // 二次迭代，提高精度 x = x * (1.5f - xhalf * x * x); // 三次迭代，提高精度 return 1/x; } |
对于神奇数字，感兴趣的朋友可以到文后下载 Chris Lomont 的论文《FAST INVERSE SQUARE ROOT》，或者到这里寻找答案。
相关下载
除特别注明外，本站所有文章均为交通人原创，转载请注明出处来自http://www.hijtr.com/fast-sqrt/
暂无评论