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

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

在一些没有硬件实现乘除法指令的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)，那么上一节你显然没看。 [关于人肉优化的必要性] 由于人肉优化不具有一般性，不是这篇文章讨论的范畴。但是考虑到但凡需要使用移位来模拟乘除的场合，往往是性能关键的部分，所以还是很有必要的。