B计划之大数乘以10，除以10的快速算法的讨论和实现

gxqcn

原帖由 gxqcn 于 2008-6-12 21:15 发表
我没有针对除10作特殊优化。

且我给的算法仅适用于基数是10的整数次幂的情形。
这是刚刚注意到的，难怪你在6#曾强调“长整数”，
因为我看到10，脑子里就默认用10的整数次幂做基数了。

上面的回复不太正确，当时脑子有点乱，也没有源代码可供参考。
昨天晚上在台式机上安装 Visual Studio 2008，我是在老婆的笔记本上登录论坛的。

我说一下我的做法，该算法通用于以任意数为基数的大数运算。

对于除数d为int型的大数除法，下面仅针对无符号型：
1、若d==0，则 HugeCalc::SetLastError( HCERR_DIV_ZERO )，return；
2、如果恰好等于大数内部基数，右移“一位”即可；
3、按小学生的除法算法，从高位向低位挺进：(上次的余数*基数+当前值)/d=商…本次的余数
但，当d为2的整数次幂时，上述除法指令可被移位指令替代优化。

对于“大整数*常规小整数”，类似，只是改为从低位向高位挺进模式进行而已。

mathe

你要做这个优化的目的是什么呢?
如果仅仅要做进制转化,是不需要这个优化的

gxqcn

进制转化确实不应这么做，因为这样做的算法复杂度太高了。
可能楼主因为经常需要碰到10进制下的“移位”处理吧。

mathe

如果这样，可能直接进行十进制计算更加好些

gxqcn

在 常规小整数<大整数基数 前提下，不同基数下的乘除法效率基本是一致的。

差异仅在于其中蕴涵的“*基数”过程是否可用移位优化，
以及因基数大小的不同，相同大整数被分割的次数。

用 n^b 为基数的大整数，在乘以（或除以）n^k 时，
当 k>>b 时，可先通过以基数的“移位”，而后再乘以（或除以）n^(k mod b) 即可。
前者是真正意义上的“移位”，是个O(1)的操作。

所以，如果算法中存在大量乘以（或除以）“10的整数次幂”运算时，
考虑以 10^b 为基数是有必要的。

最显著的例子，如大整数的输出效率。

无心人

乘以10的移位肯定是快于直接乘吧？

能有兴趣测试么？

不过，C代码在这个问题上应该不合适

gxqcn

这很难说，1次mul 与 2次shift+1次add 谁快？可能不同系统上有不同答案。

还有，大数基数*10 ≥ 2^32 ？这个在算法设计优化时也很重要。

无心人

下面仅是一些想法

如果32位系统，这个问题的话，用SSE2存在一次移位64位的算法
考虑循环展开的话，一次移位操作可以达到256位

加肯定存在一次加128位的4路循环展开算法的，这个我和liangbch已达到共识

乘也存在一次乘128位的的4路循环展开算法的，这个我也和liangbch已达到共识
但乘的逻辑复杂的多

考虑指令时间
乘是9个周期
加和移位是1－2个

如果安排合理，
按双字算，一个运算单位消耗的时钟
乘能达到12个（可能不准确）
加是5个
移位是0.5个

所以推测，至少是不输于直接乘
但运算时间最快只能达到直接乘的1/2左右

G-Spider

17# gxqcn
嗯，除了像P4之类的要考虑下，乘法其实很快了， 30# G-Spider
这里关于乘法逆元用在这个帖子里不合适，事实上其局限性很大。"Quotient will only be correct if d divides a exactly."
n_limb/1_limb，不知道郭用的什么方法，gmp中mpn_divrem_1 -- mpn by limb division.用到的方法，我没有看特别明白，可能是先求导数(选取有些讲究)，再乘法代替除法。
而这种方法用到的变量有些多，于是在x86 ALU asm中没有看到相关的实现，xmm和x64倒是有，如mpn\x86\p6\mmx\Divrem_1.asm。
mpn_invert_limb
按照它的说法
At 4 limbs the div is a touch faster than the mul (and of course  simpler), so start the mul from 5 limbs. (mul_by_inverse)
在vc6.0中也有关于一个常量的除法的优化方法：
如(unsigned int) n/3，优化为

1.         mov        eax, -1431655765                        ; aaaaaaabH
   
2.         mul        n
   
3.         shr        edx, 1

复制代码

aaaaaaabH 对应于2/3 (规范化，有助于得到高精度的表示)，于是上在的代码解释为(n*(2/3))/2.
又比如n/10

1.         mov        eax, -858993459                                ; cccccccdH
   
2.         mul        n
   
3.         shr        edx, 3

复制代码

1/10=(4/5)/8 如上，cccccccdH 为4/5，然后除以8(用右移3)。

G-Spider

19# G-Spider
上面vc做的优化是dividing a single-word number by a single-word number ,using a single-word approximate reciprocal.
这在大数除以1个limb的情况下不合适，大数中需要实现dividing a two-word number by a single-word number,using a single-word approximate reciprocal.的优化方案。
gmp有篇相关的文章Improved division by invariant integers
division-paper[1].pdf (138.97 KB, 下载次数: 1)

2013-6-5 23:05 上传

点击文件名下载附件