数学研发论坛

 找回密码
 欢迎注册
查看: 1306|回复: 7

[分享] 高精度除法中,试商法的 3/2 算法

[复制链接]
发表于 2017-6-12 05:38:56 | 显示全部楼层 |阅读模式

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

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

x
这几天花了非常长的时间阅读GMP的源码,今天我来分享一下它里面的一小部分:除法的试商法。
我们这里用 $B$ 代表高精度数的基,在 GMP 里 $B=2^32$ 或者 $2^64$, 我写的是64位程序,用后者。
一个$m$位的大整数$a$(每一位满足$0<=a_i<2^64$,下同),要去除以一个$n$位的大整数$b$,我们把$b$左移,使它的最高位满足$B/2<=b_{n-1}<B$,也就是说相当于在二进制下把除数左移到占满一个寄存器的所有bit(如果将被除数同时左移同样的位数,则商不变,最终余数只需右移复原即可,这个过程叫做正则化)。
然后定义$b$的逆:$v=\floor{\frac{2^192-1}{b_{n-1}*2^64+b_{n-2}}}$.
也就是$(B^3-1)$除以 $(b的最高两位)$ 的商。
这个192bit除以128bit的大除法也有点难度,但GMP给出了非常好(也同样很难证明)的算法解决它。我会把这个算法的汇编实现发上来,具体有时间再解释。
结果有$0<=v-2^64<2^64$, 也就是说 $v-2^64$ 是刚好可以放到一个寄存器里的,我们实际保存的是这个值。
在一次高精度除法中,这个逆只需要在做除法之前算1次。这使得它的代价其实很小。
————————————————————
假设我们已经有这个逆。
一旦有了这个值,我们就可以证明下述算法:
对高精度的$a$除以$b$,我们取$a$的最高三位和$b$的最高两位进行试商。
(在这里,冒号表示$B$进制连续数字,比如$1:2:3=1*B^2+2*B^1+3*B^0$)
对三位整数$a=a_2:a_1:a_0$和两位整数$b=b_1:b_0$,已知它们的真商小于$B$,即$a_2:a_1<b_1:b_0$,而且b经过正则化
如果$b_1:b_0$的逆是$v=\floor{\frac{2^192-1}{b_1:b_0}}$已知,则有:
1.$a_2*v+a_1$的结果是两位数$q_h:q_l$,且$q_h$总是小于或等于$a$除以$b$的真商,且至多小$2$。
2.所以我们首先试商$q_h+1$,计算对应的余数$z=a-(q_h+1)*b$,真正的$z$可能小于0或者怎么样,但我们不用在意,只考虑它模$2^128$的值,记为$z mod 2^128=z_h:z_l$。
模$2^128$在64位CPU上是简单的,考虑低两位,不管溢出就行了。
3.神奇的事情来了!
炒鸡难证(我证了三天三夜)三定理:
(以下定理中$z$代表不取模的值$z=a-(q_h+1)*b$,只在理论证明中出现。$z_h$和$z_l$代表取模之后的对应位,是我们实际汇编能得到的值)
(1).如果真商比$q_h+1$多1,那么一定有$z_h<q_l$,且$2^127<=z<2^128$;
(2).如果真商比$q_h+1$少1,那么一定有$z_h>q_l$,且$0<=z+b<2^128$;
(3).如果真商就是$q_h+1$,那么$z_h$和$q_l$之间没有确定大小关系。但一定有$0<=z<2^128$,而且,此时如果$z_h>=q_l$成立,那么$2^127<=z+b<2^128$也一定成立。
经过艰难的思考,我们总结出:
(i).如果$z_h<q_l$(可能是情况(1)(3)),那么我们取$Q=q_h+1$,$Z=z_h:z_l$,那么一定有$Z=a-Qb$(没有溢出求模了!);
(ii).否则($z_h>=q_l$,可能是情况(2)(3)),我们取$Q=q_h$,$Z=z+b mod 2^128$,那么一定仍然有$Z=a-Qb$(没有溢出求模了!)。
简直神奇!
要知道$mod 2^128$基本上是完全溢出不要命的算法!就这样都能给溢回来!此时$0<=Z=a-Qb<2^128$是真实的,丝毫没有溢出的!
那么很显然,因为正则化$2^127<=b<2^128$,$Q$至多比真商小1,比较一下余数$Z$和$b$哪个大,做最后一次调整就OK了。
(iii)如果$Z>=b$,则真商是$Q+1$,真余数是$Z-b$,否则真商是$Q$,真余数是$Z$。
整个过程的汇编代码非常简洁,不超过 60 行,包含两次 mul 和一次 imul ,以及一次条件跳转。
————————————————————
回到原话,对高精度的$a$除以$b$,我们取$a$的最高三位和$b$的最高两位进行试商。
这个商是精确的高三位除以高两位的结果,我们可以证明新的定理:
(α).高三位除以高两位的试商要么是$a$$b$的真商,要么比真商大1.
(β).对随机的大整数$0<a<=B^{n+1}$和$B^n/2<=b<B^n$(记住,我们经过了正则化),试商比真商大1的概率只有
$\frac{1}{2B}$
你没看错!对$B=2^64$,我们只有$2^-65$的概率需要去做调整!每次调整都要将整整一长串的高精度除数加在差不多一样长的高精度部分余数上,也就是说每次调整试商都是浪费!但我们现在只有不到$1/36893488147419103232$的概率需要做这样的调整!
以至于为了检验我的调整代码没有写错,做测试的时候我必须要巧妙地给出被除数和除数才能触发那个调整试商的分支Orz……

好了,现在我们就有了非常好的,一下就能试出来64位商的,基本上不需要调整的,试商法了!
End。
————————————————————
附件包含
(1).对1位数$2^63<=x<2^64$,计算$\floor{\frac{2^128-1}{x}}-2^64$的汇编过程 ilmp_inv_1_
(2).对2位数$2^127<=x<2^128$,计算$\floor{\frac{2^192-1}{x}}-2^64$的汇编过程 ilmp_inv_2_1_
(3).对三位整数$a=a_2:a_1:a_0$和已经正则化的两位整数$b=b_1:b_0$,在它们的真商小于$B$,即$a_2:a_1<b_1:b_0$的情况下,已知$b$的逆$v$,计算商和余数的汇编过程 ilmp_div_3_2_ni_
(4).一个Mathematica 8的文件,包含我对上面那三个炒鸡难的定理的证明

div_algorithm.rar

2.35 KB, 下载次数: 15, 下载积分: 金币 -1 枚, 经验 1 点, 下载 1 次

除法算法

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-7-21 13:24:23 | 显示全部楼层
已经下载,好好向您学,我看懂程序代码要很长的时间。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-4-6 17:32:52 | 显示全部楼层
其实 试商 里面  用 被除数高三位 除以 除数 高两位 进行试商 , 这个误差分析并没有你想得那么复杂.
可以说应该是很简单的.

比如说
被除数  abcxxxxxx...
除数     dfyyyyy....

这里有一些简单的条件 , 比如 ab < df  (废话, 如果ab >= df ,那直接就 两位除以两位 ,这种情况试商更简单, 略过 )
只考虑 ab < df的情况.

这里考虑极限情况, 只考虑 最大max 和最小 min的情况.
这里如果能证明 ,max - min <= 1 那就说明 试商的误差是很小的, 最多需要调整一次.
这里稍微缩放一下
abc /(df + 1) <= min <=max <= (abc + 1) / df

很幸运
(abc + 1) / df - abc/(df+1)  < 1
这个等式是恒成立的.
这个证明应该很简单 .
这里用多项式表示一下
abc ,df 用多项式表示  ax^2 + b x + c , d*x + f
代入化简一下就行了啦
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-4-6 17:37:59 | 显示全部楼层
本帖最后由 knate 于 2018-4-6 17:42 编辑
knate 发表于 2018-4-6 17:32
其实 试商 里面  用 被除数高三位 除以 除数 高两位 进行试商 , 这个误差分析并没有你想得那么复杂.
可以 ...


这里 我更偏向于 使用 极小值 来作为试商,
最小值, 能保证 试商一定是小于的, 即 被除数 - 除数 * 试商时 进行调整时, 一定是整数, 无需进行大小比较.

而且, 个人更感兴趣的是,  这个试商, 在被除数 高两位 什么情况下,
就已经不需要用高三位除以高两位试商, 而是直接使用高两位除以高一位试商.
这个是我更感兴趣的.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2018-4-26 22:52:30 | 显示全部楼层
knate 发表于 2018-4-6 17:37
这里 我更偏向于 使用 极小值 来作为试商,
最小值, 能保证 试商一定是小于的, 即 被除数 - 除数 * 试 ...

确实不复杂,至少差1很好证,但这里的难度在于找到一种方法让64位寄存器处理64bit的三位除以二位而绝不会出现溢出。
你说的估商与真值至多差1是建立在已经有三位除以二位的商的基础上的,而我给的方法正是在计算这个除法本身。
我觉得这个方法已经很贴近实现了,毕竟它给出了只有几十行的汇编算法,只包含三次乘法和一次跳转。而判断什么情况下只用算2除1估商,这个判断的代价可能都比做一次这个算法要高。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-26 00:25:46 | 显示全部楼层
问题或者关键字,高精度除法中的试商

高精度除法中的试商是重点也是难点,好的试商方法可减少调整商的次数,提高除法运算的速度。
我之前用到类似的算法,只是采用“2/1”的试商法,这个方法的缺点是试商的精度不够,导致有不小的概率需要调整商。当试商偏小时,会导致余数大于被除数,这时需要将余数减去除数,同时将商减少1。
当试商偏大时,会导致余数小于0(不够减),这时需要将余数加上除数,同时将商增加1。
我分别试用过这2种策略,实验表明,第一种方法更快一些。原因是,在我的实现中,第一种策略得到的商是真商的概率更大。

楼主的帖子最吸引人的地方是
对高精度的$a$除以$b$,我们取$a$的最高三位和$b$的最高两位进行试商。
这个商是精确的高三位除以高两位的结果,我们可以证明新的定理:
(α). 高三位除以高两位的试商要么是$a$与$b$的真商,要么比真商大1.
(β). 对随机的大整数 $ 0< a <= B^{n+1} $, 和 $\frac{B^n}{2} <= b < B^n$ (记住,我们经过了正则化),
试商比真商大1的概率只有$\frac{1}{2B}$


这2天想起了这笔欠债(早就想学习下,可一直没有去研读),重新看了楼主的帖子(尽管没有全部看懂),给出我自己的方法和误差分析,以飨读者。

1.高精度数的表示
为了简单起见,我们对被除数和除数做一些限定,如果被除数和除数的值不在下述范围,可相应的扩大或者缩小$B^i$倍,将被除数和除数的值的范围调到以下规定的范围。

我们用B表示高精度数的基,在32位系统,$B=2^32$,在64位系统,$B=2^64$, 没有特别说明的话,下文中的$B=2^32$。
对于高精度数x(被除数), 我们表示为 $x= x_n . B + x_{n-1} + x_r$, 这里x包含n位数字,这里的位不是2进制的位,也不是10进制的位,而是B进制的位。$x_n$表示最高位,$x_{n-1}$表示次高位。其余部分表示为$x_r$,显然,$x_r<1$
  
对于高精度数y(除数), 我们表示为 $y= y_m  + y_{m-1}/B + y_r$, 这里y包含m位数字。
同理$y_m$表示最高位,$y_{m-1}$表示次高位,$y_r$表示其余部分。

这里要求
1) n>m。 否则,可在x尾部补0,使其长度达到m+1。
2) 对x,y做正规化,对除数乘以$ 2^k( 0<=k<32)$ ,使得 $B/2 <= y_m < B$,同时对x也同样乘以$2^k$
3) x的高m位,小于y. 否则,可在最高位前补0.
4) y不能是2的方幂。
说明,当y是2的方幂时,x/y可以使用速度更快的移位来实现,故这个约束并不会对高精度除法带来不便。
楼主的算法中,并不要求y满足条件4,之所以附加这一要求,是因为,在我的方法中,计算y的倒数时,使用$B^2/y$,而不是$(B^2-1)/y$,这个改变可使得误差分析更加简单。当y是2的方幂时,对y做正规化后,$y_m=B/2$,其后的所有位都是0。这样计算$B^2/y$时,会导致结果多1比特,进而不能以统一的方法试商。

2.计算y的倒数
  令 $v= B^2/y$, $v/B^2$为y的倒数,我们在试商时,将使用乘法来代替除法,基于小学学过的数学知识,“除以一个数等于乘以这个数的倒数”。
  
在计算倒数时,我们不可能也无需求得倒数的精确值,我们只需要v的整数部分,当$B=2^32$,v的值的整数部分是1个33比特的数,当$B=2^64$时,v的值的整数部分是一个65比特的数。在计算v时,我们可
1. 只考虑$y_m$, 这样最简单,但得到的v偏大,需要调整商的概率较高。
2. 为了使得v的值更准确,我们这里需考虑$y_m$和$y_{m-1}$. 即使得
  $v*( y_m. B + y_{m-1})<=B^3$ 且 $(v+1)*( y_m. B + y_{m-1})>B^3$. 这样在使用“3/2”试商时,得到的商更准确,使得其后几乎不用再调整商。这个算法并不难,这里从略。

3.计算 x / y的商

3.1 误差分析基础
  若一个数的真实取值为x,估算值为y,$y<=x$,表示为=n+f(n是整数,$0 <= f<1.0$),误差$e=|x-y|$,则有
  如果误差$0.0 <= e <= 1$, 则$\floor(x)$要么是n,要么是n+1,下面是推理过程
  $0.0 <= f < 1.0$
  $0.0 <= e <= 1.0$,
  $ =>  0 <= f+e <2.0 $
  $ =>  n <= n+f+e < n +2.0 $
  $ =>  \floor(n+f+e)=n$ 或者 $\floor(n+f+e)=n+1$

3.2 初步的估商方法
   若我们已经得到v(只需要整数部分),当$B=2^32$,v的整数部分是一个33比特的数,故可表示为B+w,这里 $0 <= w <= B-1$  
   则,$ x/y ~= q_0 = ( x_n.B + x_{n-1} + x_n.w)/B$
   用这个方法得到的值偏小,因为 $q_0 < x/y< B$, 故 $( x_n.B + x_{n-1} + x_n.w)$ 是1个64bit的数,在加的过程中,不会溢出。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-26 01:38:26 | 显示全部楼层
3.3 误差分析
  现在我们详细的推导出上面的公式,并做误差分析
  在前面的计算中,我们忽略了v的小数部分,现在我们令
  v= B + w + r , 这里r为v的小数部分, r<1

则 x /y
   $ =  x * (v/{B^2})$
   $ = (x_n.B + x_{n-1} + x_r)(B +w +r)/{B^2} $
   
   $ = ( x_n.B^2 + x_{n-1}.B + x_r.B  + x_n.B.w + x_{n-1}.w + x_r.w + x_n.B.r + x_{n-1}.r + x_r.r) /{B^2} $
   
   $= ( x_n.B^2 + x_{n-1}.B + x_n.B.w)/B^2 +  (x_r.B + x_{n-1}.w + x_r.w + x_n.B.r + x_{n-1}.r + x_r.r)/B^2 $
   
   $= (x_n.B + x_{n-1} + x_n.w)/B $
      $ + (x_r.B + x_{n-1}.w + x_r.w + x_n.B.r + x_{n-1}.r + x_r.r)/B^2 $
   上式分为2部分,
   第一部分为 $(x_n.B + x_{n-1} + x_n.w)/B $ 就是我们3.3中的估商法
   第二部分为 $(x_r.B + x_{n-1}.w + x_r.w + x_n.B.r + x_{n-1}.r + x_r.r)/B^2 $ 就是误差部分
   因为
      $ x_n <= B-1 , x_{n-1} <= B-1, w<= B-1 $
      $ r<1, x_r<1 $
   
   故误差的最大值
     $ me < (B + (B-1)*(B-1) + (B-1) + (B-1).B + (B-1) + 1)/B^2 $       
     $ = (B + (B^2-2.B+1) + (B-1) + (B^2-B) + (B-1) + 1)/B^2 $
     $ = (2B^2)/B^2 = 2 $
  
故误差 $ 0 <= 误差 < 2 $, 结合 3.1,在最坏情况下, 估整商可比真整商小2,真整商可能的取值为
  $ q_0$, $q_0+1$,$q_0+2$
   
3.4 估商算法的细化
    1. 取$q=q_0$,
    2. 计算商与除数的乘积 $p= q(y_m .B + y_{m-1})$,
    3. 计算 $z = (x_n*B^2 + x_{n-1}B + x_{n-2} )-p$,
    4. 若z大于除数,做 $ z = z- ( y_m .B + y_{m-1}),  q=q+1 $
    5. 若z依然大于除数, 做 $  z = z-(y_m .B + y_{m-1}), q=q+1 $
  上述过程,在最坏情况下,调整商2次,主要代价为 2次 3位数 减 3位数的减法。

3.5 估商算法的进一步改进
   重新考虑误差部分 $(x_r.B + x_{n-1}.w + x_r.w + x_n.B.r + x_{n-1}.r + x_r.r)/B^2 $
   因为
      $ w <= B-1,r<1, x_r<1 $

故误差的绝对值的最大值
    $ (B + x_{n-1}.(B-1) + (B-1) + x_n.B + x_{n-1} + 1)/B^2 $
    $ = (B + x_{n-1}.B- x_{n-1} + (B-1) + x_n.B + x_{n-1} + 1)/B^2 $
   $ = ( ( x_{n-1} + x_n)B  + 2.B   )/B^2 $
   为了使得这个值 <= 1,则需
     $ => ( ( x_{n-1} + x_n)B  + 2.B )/B^2 <= 1 $
     $ => (  x_{n-1} + x_n)B  <= B^2-2.B $
     $ => ( x_{n-1} + x_n) <= B-2 $
  故当 $x_n + x_{n-1} <= B-2 $ 时,误差部分不超过1
   
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 4 天前 | 显示全部楼层
不错
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2018-12-11 01:03 , Processed in 0.055593 second(s), 19 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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