对32位整数除法优化原理的整理

shines

原帖是@liangbch发的：

http://bbs.emath.ac.cn/thread-521-3-1.html

这里我重新整理了一下，更清晰了一点，关于误差的分析更完整了。


在32位计算机里，计算 $a$ 除以 $c$，我们定义商为 $q$ 和余数为 $r$，即：$a / c = q  ...  (r)$。

推演步骤：

1. 找到一个数 $k$，另 $2^k < c < 2^(k+1)$；

2. 两边求倒数，即： $1/2^k > 1/c > 1/2^(k+1)$ ，两边都乘以$2^(k+32)$，

    得到：$2^32 > 2^(k+32)/c > 2^31$，

    即：   $2^31 < 2^(k+32)/c < 2^32$。

3. 另 $f = 2^(k+32)/c$，其中 $f$ 为浮点数（实数），对 $f$ 向上或向下取整，得到：$s = |f| = |2^(k+32)/c|$ ，这里 $|f|$ 定义为向上或向下取整；

4. 计算： $q'' = (a * s)$，这里 $q''$ 是一个被放大了约 $2^(k+32)$ 倍的 $q$ 的粗略值；

5. 计算： $q'  = q''$ >> $(k + 32)$，即把 $q''$ 向右位移 $(k+32)$ 位，得到的 $q'$ 是一个接近真实 $q$ 值的粗略值；

6. 到此为止，$q'$ 可能刚好等于真实值 $q$，或者略大于 $q$ (当 $s = |2^(k+32)/c|$ 为向上取整时)，或者略小于 $q$ (当 $s = |2^(k+32)/c|$ 为向下取整时)；

7. 计算余数：$r = a - q' * c$，如果 $r$ 值不在 $[0, c)$ 的范围内 (即 $r < 0$ 或 $r >= c$)，则修正 $q'$ 和 $r$ 的值；


下面我们来做误差分析：

$q' = (a * s) $ >> $(k + 32)$，而 $s = |f| = |2^(k+32)/c|$，假设为向上取整，设 $s$ 与 $f$ 的误差为 $e$，即 $e = s - f$；

其中 $q =  (a * f)$ >> $(k+32)$，

那么 $q$ 的误差 $E = q - q' = (a * (s - f))$ >> $(k+32) = (a * e)$ >> $(k+32)$，

由于 $a / c$ 的余数最小值为1，即只要满足总误差 $E < 1 / c$，则能保证 $q' = q$，也即：

[ $(a * e)$ >> $(k+32)$ ] $< 1 / c$，

由于我们的被除数是32位的整数 $a$，则有 $0 <= a < 2^32$，上式可简化为：

[ $e$ >> $k$ ] $< 1 / c$，即：$e / 2^k < 1 / c$，可得到：$Rightarrow$ $e < 2^k / c$；

我们知道，如果 $|f|$ 是向上取整，则有 $0 <= e < 1$，根据推演步骤1有 $2^k < c < 2^(k+1)$，则 $2^k / c < 1$，

即只要满足 $0 <= e < 2^k / c$ 即可，定义为（1）式。

由于只要 $c$ 一旦确定，然后确定是向上取整或向下取整后，$e$ 也是确定的，如果向上取整不能满足（1）式，
则可以尝试向下取整，如果向下取整也不能满足，则可以搜索下一个 $c$ 值，直到能够满足（1）式为止。

无心人

对最新的x86 CPU已经不需要这么麻烦了，调整的过程，已经比直接除法慢了

shines

不用调整的啊, 这里的除都是针对常量的除法, 对于非常量, 只能老老实实用除法的

shines

哦, 你是说乘法+位移的调整吗? 实际上并不慢, 原因是除法指令时钟周期是不固定的,
而且总体上还是比这个乘法+位移的过程慢, 所以还是比较有价值的

除法时钟周期不固定的原因是, 除法有两种实现方式, 一种是试商法, 一种是乘以除数的倒数,
即 $a * ( 1/c )$, 这个求$(1/c)$的过程花的指令也是不固定的(因为涉及精度的问题, 如果用这种方式,
必然是用浮点除法模拟整数除法的, 所以有精度的问题)

shines

更详细的推理和论述过程我写在这里了, 欢迎批评指正

http://www.cnblogs.com/shines77/p/4189074.html

shines

在网友BillGan的提醒下, 重新修改了一下误差分析的部分:

根据余数公式有：$r = a - q * c$，两边除以 $c$，然后再写成标准的商和余数的形式，有：

$a$ $/$ $c = q + r / c$，这里 $q$ 和 $r$ 均为整数，且 $0 <= r < c$ 。

同样的， 我们也可以把 $q$ 换成 $q'$，由于是向上取整，所以有 $q' >= q$，

我们设 $q'$ 和 $q$ 之间的误差为 $E$，并称之为 $q$ 误差：

$ E = (q' - q)$，且 $E >= 0$。

为了解释方便，我们以 /10 为例：

当 $a$ $/$ $10 = q + 0 / 10 = q + 0.0$ 时，只有当误差 $E >= 1.0$ 时，$q$ 的值才会增大 1；
更为极端的情况：
当 $a$ $/$ $10 = q + 9 / 10 = q + 0.9$ 时，只有当误差 $E >= 0.1$ 时，$q$ 的值才会增大 1；

（这里我们不用考虑 $E < 0$ 的情况，因为定义的时候保证了 $E >= 0$）。

综上可得，只要满足 $E<1/c$，就能保证取整后的 $q'$ 的值跟 $q$ 值是一样的。

由前面的公式可得：

$q' = (a * s) $ >> $(k + 32)$，而 $s = |f| = |2^(k+32)/c|$，假设为向上取整，设 $s$ 与 $f$ 的误差为 $e$，即 $e = s - f$；

其中 $q =  (a * f)$ >> $(k+32)$，

那么 $q$ 的误差 $E = q‘ - q = (a * (s - f))$ >> $(k+32) = (a * e)$ >> $(k+32)$，

只要满足 $q$ 误差 $E < 1 / c$，则能保证 $q' = q$，也即：

[ $(a * e)$ >> $(k+32)$ ] $< 1 / c$，

由于我们的被除数是32位的整数 $a$，则有 $0 <= a < 2^32$，上式可简化为（这一步很关键，但不详细解释了）：

[ $e$ >> $k$ ] $< 1 / c$，即：$e / 2^k < 1 / c$，可得到：$Rightarrow$ $e < 2^k / c$；

我们知道，如果 $|f|$ 是向上取整，则有 $0 <= e < 1$，根据推演步骤1有 $2^k < c < 2^(k+1)$，则 $2^k / c < 1$，

即只要满足 $0 <= e < 2^k / c$ 即可，定义为（1）式。

由于只要 $c$ 一旦确定，然后确定是向上取整或向下取整后，$e$ 也是确定的，如果向上取整不能满足（1）式，

则可以尝试向下取整，后面我们会论述为何向上取整或向下取整必然会有一种会满足（1）式。

shines

我们得到的最终结论是：

在32位计算机里，计算 $a$ 除以 $c$，我们定义：商为 $q$ 和余数为 $r$，即：

$a$ $/$ $c = q  ...  (r)$，其中 $q, r$ 均为整数，且 $0 <= r < c$。


设相乘系数为 $m$，向右位移位数为 $k$：


先找到一个正整数 $k$，令 $2^k < c < 2^(k+1)$，

则 $m = |2^(k+32) / c|$，其中 $|x|$ 定义为向上取整或向下取整，

$k$ 则为向右位移的位数。

liangbch

我的几点建议：

1. 以乘代除法不仅适用于除数是常数的情况，也适用于除数是相对不变的数的情况，如被除数是一个32位多重精度数，除数为32位单精度数（下文中的单精度和双精度均值32位数单精度和32位双精度数）。
2. 在实践中，被除数是单精度数的情况很少见，在绝大多数情况下，被除数是双精度数，而楼主的算法没有覆盖到这种情况。
3. 一般的，若被除数是a,除数是c，则a/c的算法为


1. 找一个数s,使得s满足以下条件
    a) s = 2^k / c 向上或者向下取整，k是整数
    b) s 是一个32比特的数

2. 找一个数k1, 使得 (a>>k1) < 2^32。 k1的可根据a的最大值来选择，在通常情况下，a的最大值是可以事先确定的。
  
3. 计算 ((a >> k1) * s ) >> (k-k1) 或者 ( ((a >> k1)+1) * s ) >> (k-k1)来得到a/c的商的近似值

误差分析，
  在整个过程中，有2处近似，
    1#, 求s的过程，因s是一个整数，2^k / c是一个有理数，s为 2^k / c的一个近似。
    2#, 因移位运算是整数运算，故 (a >> k1) 或者 (a >> k1)+1 是 a/(2^k1)的一个近似。
这两处运算中，每步运算的误差都可能是正的或者负的。可根据需要选择合适的取整方式.
如果1#中使用向上取整，在2#中亦使用向上取整，且总的误差的最大值小于1/c时，则不需要调整商。

注：
1. 在除数为单精度常数时，步骤1和2可根据被除数的范围和除数的值预先确定，具体到程序，就是hard code.
2. 在多重精度数除以单精度数变量时，在被除数位数很多的情况下，以除代乘法仍是值得做的。需要使用步骤1和步骤2来计算确定s,k和k1的值。
   故步骤1和步骤2只需1次，而步骤3需要重复多次。
3. ((a >> k1)+1) 的另一种形式形式是 (a + 2^k1-1)>> k1, 这样，当a的最末k1位是全0时，误差为0，比((a >> k1)+1)更好。

只是呼吸

与各位探究：32位（二进制）以内的两数相除，由于数的大小在机器字长之内，直接除就行。
超过32位的大整数除法可以用数组，方法和学生的竖式除法一样。例如：8000位十进制数除以4000位十进制数。那么除的结果商就约4000位十进制数。（更大的数应该采用牛顿法。）
步骤：1）输入8000位十进制数，存入字串A中。
      2）输入4000位十进制数。存入字串B中。
      3）将字串A中的数转化为数字，按顺序存入数组AA[2000]中，每个AA[J]存入四个有效数字，这种方法就是传说中的万进制。
      4）同3），将字串B存入数组BB[1000]中（万进制），但是这个数组的第一个元素（最高位BB[0]）是五位十进制数字，其余的元素是四位十进制数字。
      5）除法开始：
         A。取试商。取数组AA的前二位，AA[0]，AA[1]，数组BB的第一位，BB[0]。由计算机计算（AA[0]X10000+AA[1]）/BB[0]。这个计算是成立的，因为没有超过计算机的字长。得到的结果记为Q1。它有一个特点：它等于AA/BB的前四位商，或者等于AA/BB的前四位商减1。
             B。做乘法。用Q1乘数组BB的各位，将结果保存在数组CC中。共有1001项，要做1000次机器字长内的乘法，以及2000次机器字长的除法（必要的进位开销）。
            C。做减法。取AA数组的前1001项与CC相减。差保存在数组DD中，（这个DD数组最多只有1000项）。这一步要做1001次机器字长内的减法（以及适当的进位开销）。
            D。做判断。IF（DD大于等于0）则Q1就是精确商。
                             把DD数组与AA中还没有参与计算的数，组成新的AA，返回5）。
                       IF（DD小于0）则让Q1=Q1-1。这时Q1就是精确商。
                            做加法：把这个负的DD与BB相加，结果存入DD中，把DD数组与AA中还没有参与计算的数，组成新的AA，返回5）。这一步要做1000次机器字长内的加法（以及适当的进位开销）。

   说明：在转到5）的时候，要注意一个细节，那就是有可能会商0，在编程时需要特别留意。
          加减法的次数是大约数。
                                       
                                       

liangbch

32位（二进制）以内的两数相除，由于数的大小在机器字长之内，直接除就行。

我们讨论的不是可不可以直接除的问题，而是如何做，可使得除法的速度更快。
在32位模式，"64位除以32位也可以直接除，只要商的范围小于2^32。
在64位模式，"128位除以64位也可以直接除，只要商的范围小于2^64。

在绝大多数CPU，乘法速度是远快于除法的，这正是我们讨论“以乘代除法”的意义所在。例如：

在 Intel Sandy Bridge 系列CPU
32bits 乘以 32 bits 寄存器，延时为4个周期，每2个周期可发射1条指令
64bits 乘以 64 bits 寄存器，延时为3个周期，每1个周期可发射1条指令

与之相对的，
64bits   除以 32 bits  寄存器，延时为20-28个周期，每11-18个周期可发射1条指令
128bits 除以 64 bits  寄存器，延时为30-94个周期，每22-76个周期可发射1条指令

关于各个指令的周期等信息，请参阅Agner Fog 的 《Instruction tables 》，
该书 Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs