找回密码
 欢迎注册
查看: 20407|回复: 14

[分享] 对32位整数除法优化原理的整理

[复制链接]
发表于 2014-12-9 19:56:18 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
原帖是@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)式为止。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-12-12 16:32:50 | 显示全部楼层
对最新的x86 CPU已经不需要这么麻烦了,调整的过程,已经比直接除法慢了
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-12-16 22:25:49 | 显示全部楼层
不用调整的啊, 这里的除都是针对常量的除法, 对于非常量, 只能老老实实用除法的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-12-21 16:34:09 | 显示全部楼层
哦, 你是说乘法+位移的调整吗? 实际上并不慢, 原因是除法指令时钟周期是不固定的,
而且总体上还是比这个乘法+位移的过程慢, 所以还是比较有价值的

除法时钟周期不固定的原因是, 除法有两种实现方式, 一种是试商法, 一种是乘以除数的倒数,
即 $a * ( 1/c )$, 这个求$(1/c)$的过程花的指令也是不固定的(因为涉及精度的问题, 如果用这种方式,
必然是用浮点除法模拟整数除法的, 所以有精度的问题)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-12-28 04:21:59 | 显示全部楼层
更详细的推理和论述过程我写在这里了, 欢迎批评指正

http://www.cnblogs.com/shines77/p/4189074.html
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-12-29 00:33:10 | 显示全部楼层
在网友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)式。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-12-31 17:55:41 | 显示全部楼层
我们得到的最终结论是:

在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$ 则为向右位移的位数。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-1-5 11:18:50 | 显示全部楼层
我的几点建议:

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)更好。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-1-14 10:48:42 | 显示全部楼层
本帖最后由 只是呼吸 于 2015-1-14 12:13 编辑

与各位探究: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,在编程时需要特别留意。
          加减法的次数是大约数。
                                       
                                       


毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-1-14 12:24:42 | 显示全部楼层
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

点评

唉,看了你写的这个数据,除法真是慢恼火了。  发表于 2015-1-15 09:46
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 04:09 , Processed in 0.048288 second(s), 17 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表