数学研发论坛

 找回密码
 欢迎注册
查看: 797|回复: 4

[分享] 高精度除法中,试商法的 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, 下载次数: 10, 下载积分: 金币 -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估商,这个判断的代价可能都比做一次这个算法要高。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2018-5-27 21:54 , Processed in 0.056959 second(s), 19 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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