在科学运算、图形学、游戏等很多领域中,开方是很常见却又非常耗时的运算,因此必须使用快速(有时还要求准确)的开方算法。
说起开方算法我们一般想到的是牛顿迭代法,这里我介绍一种更好的方法——逐比特确认法。
逐比特确认法从数字的本质出发,关注结果的每一比特位。它从最高位开始,向低位逐一确认某位是0还是1。在数字很大时这种方法的速度比牛顿法快不少。
要理解这种方法,得先了解二进制乘法。例如,对于数字10(二进制为0b1010
),平方为100(二进制为1100100
),它的二进制平方运算过程为:
1010 X 1010 ___________ 1010x1000 + 1010x 10 =========== 1000x1000 (*1) + 10x1000 (*2) + 1000x 10 (*3) + 10x 10 (*4) =========== 1000000 + 10000 + 10000 + 100 =========== = 1100100
开方则需要我们反过来,已经有结果N = 1100100
,判断根sqrt
的二进制:
首先1100100
有8位,可以判断sqrt
起码有4位且不超过4位。如果sqrt有5位,那么仅最高位10000*10000 = 1 0000 0000就已经大于N;如果sqrt只有3位,即使sqrt为111结果110001也不超过6位。
现在判断sqrt
的第4位:如果第4位为1 , sqrt
平方运算中有上面(*1)这项
1000 X 1000 ========== 1000x1000 (*1) ========== = 1000000 (n4) < 1100100 (N)
结果 n4 < N
。容易判断,第4位一定为1。不然乘不出N
这么大的数。
现在判断sqrt
的第3位:如果第3位为1,则sqrt
为1100,它的平方为
1100 X 1100 ========== 1000x1000 (*1) 100x1000 1000x 100 100x 100 ========== = 10010000 (n43) > 1100100 (N)
结果n43 > N
,所以这一位不是1,只能是0。
到目前为止其实都是二分法的思路,先是2^3,然后是2^3 – (2^3 + 2^2),这样逐次将范围减半。
但是这里有个问题,后面每次都求了sqrt
的平方,其实重复求了之前求过的一部分,例如在第二步中,我们算了 1000x1000(*1)
,这其实就是第一步中的算的。如果我们每次算完平方,确认了这一位为1后,就从N
中减去这一部分的平方,那么下次比较大小的时候就可以少算这一位。
我们从第二步重新开始:
我们第一步确认了1000
, 从N中减去它的平方 N = 1100100 - n4
,结果为N
= 1000100
。
如果第3位为1, 那么sqrt
= 1100
, 已确认的为1000
, 正在确认的为100
, 平方为:
1100 X 1100 =========== (1000x1000)(已确认部分从N减去了,不计算) + 100x1000 (正在确认的*已确认的) + 1000x 100 (已确认的*正在确认的) + 100x 100 (正在确认的*正在确认的) =========== 2*(1000<<2) (1000*100 等于将1000左移2位) + 100<<2 =========== = 1010000 (n3) > 1000100 (N*)
和之前结果一样, 大了,所以第3位为0. 因为是0, 所以没必要从N*
里减去.
现在判断第2位: 如果为1 则sqrt = 1010.
1010 X 1010 =========== (1000x1000)(已确认部分从N减去了,不计算) + 10x1000 (正在确认的*已确认的 = 将已确认部分前移1位) + 1000x 10 (已确认的*正在确认的 = 将已确认部分前移1位) + 10x 10 (正在确认的*正在确认的 = 将正在确认的前移1位) =========== 2*(1000<<1) + 10<<1 =========== = 1000100 (n2) = 1000100 (N*)
n2 = N*
, 也就是说若这一位为1, sqrt
就是N
的根. 后面应该都是0,无需继续判断.
但我还想继续探究, 继续把N
减去新确认的部分: N*
= 1000100 – n2 = 0。
如果第1位为1,则sqrt= 1011, 平方运算为:
1011 X 1011 =========== (1000x1000) (已确认部分从N减去了,不计算) ( 10x1000) (已确认部分从N减去了,不计算) (1000x 10) (已确认部分从N减去了,不计算) ( 10x 10 ) (已确认部分从N减去了,不计算) + 1x1010 (正在确认的*已确认的 = 将已确认部分前移0位) + 1010x 1 (已确认的*正在确认的 = 将已确认部分前移0位) + 1x 1 (正在确认的*正在确认的 = 将正在确认的前移0位) =========== 2*(1010<<0) + 1<<0 =========== = 10101 (n2) > 0 (N*)
所以这一位肯定只能为0. 最终结果为sqrt
= 1010.
这就是逐比特确认法。
说了这么多,其实代码很简单:
1 int sqrt_bv(int n) 2 { 3 int sqrt = 0; 4 int shift = 15; 5 int sqrt2; //已确认部分的平方 6 while (shift >= 0) 7 { 8 sqrt2 = ((sqrt << 1) + (1 << shift)) << shift; 9 if (sqrt2 <= n) 10 { 11 sqrt += (1 << shift); 12 n -= sqrt2; 13 } 14 shift--; 15 } 16 return sqrt; 17 }
此文章首发于我的个人网站:三种高效的整数开平方算法