在编程开发中,经常需要去计算一个数的平方根,小编一般是直接调用系统函数 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/
暂无评论