高精度除法中，试商法的 3/2 算法
这几天花了非常长的时间阅读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的文件，包含我对上面那三个炒鸡难的定理的证明
已经下载，好好向您学，我看懂程序代码要很长的时间。 其实 试商 里面用 被除数高三位 除以 除数 高两位 进行试商 , 这个误差分析并没有你想得那么复杂.
可以说应该是很简单的.
比如说
被除数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
代入化简一下就行了啦 本帖最后由 knate 于 2018-4-6 17:42 编辑
knate 发表于 2018-4-6 17:32
其实 试商 里面用 被除数高三位 除以 除数 高两位 进行试商 , 这个误差分析并没有你想得那么复杂.
可以 ...
这里 我更偏向于 使用 极小值 来作为试商,
最小值, 能保证 试商一定是小于的, 即 被除数 - 除数 * 试商时 进行调整时, 一定是整数, 无需进行大小比较.
而且, 个人更感兴趣的是,这个试商, 在被除数 高两位 什么情况下,
就已经不需要用高三位除以高两位试商, 而是直接使用高两位除以高一位试商.
这个是我更感兴趣的. knate 发表于 2018-4-6 17:37
这里 我更偏向于 使用 极小值 来作为试商,
最小值, 能保证 试商一定是小于的, 即 被除数 - 除数 * 试 ...
确实不复杂，至少差1很好证，但这里的难度在于找到一种方法让64位寄存器处理64bit的三位除以二位而绝不会出现溢出。
你说的估商与真值至多差1是建立在已经有三位除以二位的商的基础上的，而我给的方法正是在计算这个除法本身。
我觉得这个方法已经很贴近实现了，毕竟它给出了只有几十行的汇编算法，只包含三次乘法和一次跳转。而判断什么情况下只用算2除1估商，这个判断的代价可能都比做一次这个算法要高。
问题或者关键字，高精度除法中的试商
高精度除法中的试商是重点也是难点，好的试商方法可减少调整商的次数，提高除法运算的速度。
我之前用到类似的算法，只是采用“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的数，在加的过程中，不会溢出。
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
不错 Ickiverar 发表于 2018-4-26 22:52
确实不复杂，至少差1很好证，但这里的难度在于找到一种方法让64位寄存器处理64bit的三位除以二位而绝不会 ...
好久没上来了!
其实觉得只要证明出3/2的试商误差不会超过2的话, 我觉得这个办法已经是可以结束了!
从整个除法来说. 这个3/2试商占的比例很低很低, 肯定不会超过1%时间. 甚至可能是千分之几,或者更低比例.
换句话来说. 即使试商优化至无需时间, 那实际对整体来说, 也是无关痛痒的.
至于是64位数据,还是128位甚至是192位, 这些在我看来, 都是非常细枝末节的事情.
能优化就优化, 没优化, 只要结果是准确,那就没有问题了!
既然提到优化, 我把所知道的3/2 算法一些优化简单提一下!
假定一下条件
现在求 abc/df(3位除以2位,结果是1位数字 ,)
(一般来说, 估商 被除数舍弃尾数往上进1,除数直接舍弃位数, 则估商偏大 ,
被除数直接舍弃尾数,除数舍弃尾数往上进1, 则估商偏小
后面我习惯是往偏小的估算, 如果往偏大估算,则反过来就好了,区别不大
)
方法1:
f = ab/(d + 1);
abc = abc - df * f ;
至 abc < df为止, 这个优点是逻辑简单, 但是ab/(d+1) ,当d很小时,误差很大,因此可能需要迭代多次才能得到准确结果.
方法2:
k = 1/(df + 1) (k实际是小数, 存储时,可以先放大若干倍, 存储为整数, 这样就不会出现浮点数)
f = ab * k;(f是整数)
abc = abc -df * f
至 abc < df为止, 这个里面有个小小的分支,
就是k到底使用 1位有效数据还是2位有效数据,
1位速度快, 但是误差稍微更大, 有时可以少迭代一些次数,
2位 速度稍慢, 毕竟计算f时需要计算2 * 2的乘法,但是估商误差会极小, 测试过, 基本上1次就能差不多达到要求.自己取舍.
方法3:
这个似乎没什么人提到的, 不知道为什么.(或者可能是这个限制条件比较多, 应用范围不大)
就是 在基R很小的时候, 比如 是2进制, 3进制, 10进制, 等等, 一个机器位长能能保持 基的若干次(起码超过2次)幂的情况下, 可以做一些优化处理.
打个比方,比如现在能存 10进制 的 4次幂,
那么 当d 很小, 比如 小于10的2次幂时,
那么可以把改写一下 abc , df这两个数,
ab = ab * 10^2+ c / 10^2;
d = d * 10^2 + f / 10^2 ;
f =ab / (d + 1) ;
abc = abc - df * f ;
这样做的根本原因在于 当除数变大时,自然而然, 试商时的误差就变小, 到时需要迭代的次数就自然变小.
有个原则就是 除数被除数放大时, 需要保证 被除数至多是2位数, 除数至多1位数,不一定是对半截断. 可以根据情况自行处理
PS: 这帖子居然吸引了梁宝宝过来! , 厉害厉害! 还有一个估计是更好的点子.
基本确定, 3/2的估商, 是准确的.
如果 后面的算法, 得到的商, 和3/2的商一致, 那可以认为 估商是准确的了.
根据 kunth 上提到的.
2/1 估商. ab/d (B进制)
只要除数 b >= B/2.
误差就不会超过2.
只要我把d的精度再提高一位以上, 后继使用2/1 估商,误差也就不会超过2/2 = 1 ,即误差也就不会超过1了.
即估商变成了精确的了.
abc/df(3/2估商)
一般来说. B进制一般不会使用超过2^31.(32位系统)
这里有两个方式处理.
一: 找出一个系数 k = (2^64 - 1) / df ,使得 2^63 <= df*k < 2^64
计算
(abc *k) / (df *k)
由于df*k < 2^64, 依然是 3/2,位 因此 使用2/1估商,
二: 找出一个系数 p = (2^64 - 1) /(ab+1) 使得 2^191 <= abc*p < 2^192
很容易知道 df * p < 2^64,
即 估商时,只要计算一次 int64/int64即可.
理论上, 方法二比方法一高不少.
但是由于误差都不会超过1, 实际结果没什么区别.
由于计算 k,p时实际无需精确的,
因此, 实际处理时, 可以简单实用移位接口能够实现.
实际计算,
大致计算量, 都是一次int64/int64,加3-5乘法次,和若干次加减和移位运算.
应该能够在30条指令内完成.
页:
[1]
2