一种高效整数开平方算法：逐比特确认法

Posted 2021-03-14 h5l0

在科学运算、图形学、游戏等很多领域中，开方是很常见却又非常耗时的运算，因此必须使用快速（有时还要求准确）的开方算法。

说起开方算法我们一般想到的是牛顿迭代法，这里我介绍一种更好的方法——逐比特确认法。

逐比特确认法从数字的本质出发，关注结果的每一比特位。它从最高位开始，向低位逐一确认某位是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 }

 

此文章首发于我的个人网站：三种高效的整数开平方算法