移位实现常量乘除的简单优化

移位实现常量乘除的简单优化

在一些没有硬件实现乘除法指令的CPU中，我们常用移位操作来取代常量的乘除操作，运算结果向下取整。因此，位移的次数和运算的精度成为我们关注的重点。这里来谈谈几个显而易见，但是你也许会忽略的问题。

整数乘法和小数乘法

---

用位移实现常量整数乘法、小数乘法，都不会遇到精度问题，这是根据二进制数本身的原理来的。运算次数则根据转换算法的不同而不同。

[算法1：二进制数展开，只用加法]
最简单的方法就是把乘数展开，逐位相加。
步骤1：乘数转换成二进制
步骤2：凡是为“1”的位置，作为展开后的一项，移位的方向和次数与“1”的位置相同
步骤3：相加

例：9 = 0000 1001.      => (x << 3)       + (x << 0)例：0.625 = 0.1010      => (x >> 1)       + (x >> 3)

[算法2：加减法逐次逼近]
步骤1：从最高位（MSB）开始逐次递减，考察该项与乘数的距离最高位（MSB）开始逐次递减，考察该项与乘数的距离
步骤2：选取最接近当前值得位（Bit）
步骤3：如果该位大于当前值，符号位选取“+”号，否则选取“-”号
步骤4：求和

例：7 = 8 - 1（按算法1则是1 + 2 + 4）      =>  (x << 3)        - (x << 0)例：0.21875 = 0.25 - 0.03125（按算法1则是0.125 + 0.0625 + 0.03125）       =>  (x >> 2)        - (x >> 5)

[算法1与算法2的评价]
直观的看过去，算法2利用了四舍五入的原理，因此优于算法1。我们来验证一下：

用例1：0到65535，计算位宽uint16_t。

Better: 46736
Worse : 0
Equal : 18800

用例2：1/1到1/65535，计算位宽uint16_t。

Better: 4847
Worse : 0
Equal : 60688

整数除法和小数除法

---

通常认为除法是乘法的逆运算，通过取得除数的倒数与原数相乘，就相当于除法。但是这一规则并不完全适用于整数除法。单独利用算法2实现除法会遇到运算精度问题，这是由于求整数的倒数会丢失精度。

[精度问题]
例如我们要计算x / 800。利用算法2，在uin16_t下我们得到：

    y = (x >> 10) + (x >> 12)

但这不是x / 800的精确解。原因在于1 / 800的完整展开式为：

0.000000000101000111101011100001010001111010111000010100011110101110000...

这个误差是多大呢？以0-65535的uint16_t来说：

Error 0: 6752
Error 1: 27328
Error 2: 25984
Error 3: 5472
Error>3: 0

当然对于精度要求不高的场合，这也足够了。

[算法3：预先左移]
好在完成除法之后，我们只需要整数部分，所以过高的运算精度毫无意义。那么多高的精度能够满足要求呢？事实上只需要预先左移适当的位数，就能大大提高运算精度。

步骤1：除数转换为二进制
步骤2：以第1次出现1的位置作为预处理阶段左移的数量
步骤3：除数左移
步骤4：使用算法2
步骤5：运算结果右移

以x / 800为例，首次出现1是在小数点后第10位，所以选择左移9位（如果左移10位就溢出了）。然后回到算法2即可。运算结果：

    x =   (x >> 1)  + (x >> 3)  + (x >> 6)        - (x >> 11) - (x >> 13);    x >>= 9;

此时在0-65535以及uint16_t下的误差有所下降：

Error 0: 65524
Error 1: 12
Error>1: 0

你也许会问，预先左移是否导致被除数溢出？这个问题太幼稚回家自己想，想不出来打屁股。

[算法4：循环节优化]
我们用膝盖想也知道有理数除法 + 二进制转换不可能产生无理数。细心的你也许会发现，我上面举的1 / 800的例子并非毫无规律：

1 / 800 = 0.000000000
10100011110101110000
10100011110101110000
10100011110101110000...

[作用1：提高运算精度]
我们来试着处理uint32_t的情况。使用算法3得到：

    x =   (x >> 1)  + (x >> 3)  + (x >> 6)        - (x >> 11) - (x >> 13) - (x >> 16)        + (x >> 21) + (x >> 23) + (x >> 26)        - (x >> 31);    x >>= 9;

在0-4294967295之中的误差数目是：

Error 0: 4293047518
Error 1: 1919778（万分之4）
Error>1: 0

由于意识到循环节的存在，我们改为

    x =   (x >> 1)   + (x >> 3)  + (x >> 6)        - (x >> 11)  - (x >> 13) - (x >> 16);    x =   (x << 0)        + (x >> 20);    x >>= 9;

现在0-4294967295之中的误差数目是：

Error 0: 4294567686
Error 1: 399610（十万分之9）
Error>1: 0

看来 精度提高了4.8倍。

[作用2：缩短移位次数]
循环节越短，缩短就越明显。以1 / 3为例，在uint_16的情况下：

1 / 3 = 0.010101010...

使用算法3得到：

    x =   (x >> 1)  + (x >> 3)   + (x >> 5)  + (x >> 7)        + (x >> 9)  + (x >> 11)  + (x >> 13) + (x >> 15);    x >>= 1;

现在利用算法4可以二分展开：

    x = (x >> 1) + (x >> 3);    x = x + (x >> 4);    x = x + (x >> 8);    x >>= 1;

在这个例子中，算法4的精度同样高于算法3，但我想这应该是运气因素。如果你问我为什么不直接(x >> 2) + (x >> 4)，那么上一节你显然没看。

[关于人肉优化的必要性]

---

由于人肉优化不具有一般性，不是这篇文章讨论的范畴。但是考虑到但凡需要使用移位来模拟乘除的场合，往往是性能关键的部分，所以还是很有必要的。