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

liangbch

正好我前一段时间在学习和研究 使用以乘代除法，即计算一个数除以一个不变的数的除法。我见到如下2篇文章，本打算通读的，但后来将注意力集中到自己实现除法上面来了，就没有深入。 现将我的方法简要介绍一下： 目标：在32位计算机，计算a除以c的商q和余数r。 步骤： 1. 找一个数k，使其满足 $2^k> (k+32) (将q右移k+32比特） 　　　到此为止，得到的q可能等于，略大于(当s为$2^(k+32)/c$ 向上取整时），略小于(当s为$2^(k+32)/c$ 向下取整时) 5. 计算余数$r=a-(qxxc)$，若r小于0或者大于c,调整余数r和商q。 当被除数和除数是单精度整数时，在一定条件下，得到的q可为准确值，不需要做步骤5. 其步骤和约束条件为： 若 $f= 2^(k+32)/c$，f是浮点数，将f向上取整得到s, e为误差，e=s-f 容易知道，在计算a*s的过程中，误差为$a**s - a**f= a** (s-f) = a**e$. 我们知道，浮点运算a/c得到的值是离散的，例如浮点运算a/10的值，小数部分总是0.0，0.1, 0.2, 0.3,0.4,0.5,0.6,0.7,0.8,0.9, 相邻的2个值的差为1/c=1/10=0.1. 故若q用 a*s>>(k+32)这样的做法来计算，如其误差$a**e<1/c$，则得到的q总是正确的。 再举例做进一步说明，如浮点运算$a/10=q+0.0$, 当其误差$a**e>=1$时，才会影响到q的值，若浮点运算$a/10=q+0.9$， 当误差$a*e>=0.1$时，才会影响q的值）。

liangbch

举2个实际的例子 1. 计算a/10的商q，VC的算法 q= (a*0xcccccccd)>>(32+3)，下面我们分析其原理 因为 $(1/10)**2^3=0.8<1$ 且 $(1/10) ** 2^4=1.6 >1$, 故取k=3 $f=0.1**(2^35)= 3435973836.8$, s=ceil(f)=3435973837，16进制表示为0xCCCCCCCD. 误差为 $(a**0.2)/(2^35)$, 解不等式$(a**0.2)/2^35<0.1$得，$a<2^34$。 故当a为32bit数时，这个方法得到的商总是正确的。 2. 计算a/3的商q，csdn中曾讨论过这个问题，我也参与其中，但我当时没有给出证明。 因为 $(1/3)**2^1=0.66667<1$ 且$(1/3)**2^2=1.33333>1$, 故取k=1 $f=0.33333333333333**(2^33)= 2863311530.6666667$, s=ceil(f)=2863311531. 误差为$(a**0.3333333)/(2^33)$, 解不等式$(a**0.333333333)/(2^33)<1/3$得，$a<2^33$。 故当a为32bit数时，这个方法得到的商总是正确的。

G-Spider

是的，但是以上两个例子可以得心应手的实施，是利用了mul指令特性，把被除数放到eax，除数的相关倒数放到edx，做乘法，可以看到商进到了edx，余数进到了eax，再做简单调整可以得到正确结果。问题是大数需要64位的被除数与32位的除数，得到商和32位余数。所以也不太容易解决，上面的附件讨论的应该就是这个问题。

liangbch

我的实验结果是，若忽略计算除数的倒数的时间，则以下三个方法均比直接使用除法指令更快,故使用此法计算多精度整数除以单精度整数比直接使用除法指令更快。 A. 计算UINT32 除以 UINT32的商 B. 计算UINT32 除以 UINT32的商和余数 C. 计算UINT64 除以 UINT32的商和余数 当然，需要花费的时间 t(A)

G-Spider

在G大数库中，用以乘代除的方法求mutil limb /one limb得到商和余数，经测试，有性能提升:
mpz_div_d性能测试.zip (6.4 KB, 下载次数: 10045)

2013-6-30 13:52 上传

点击文件名下载附件

.·.·.

liangbch 发表于 2013-6-6 00:34
举2个实际的例子
  1. 计算a/10的商q，VC的算法 q= (a*0xcccccccd)>>(32+3)，下面我们分析其原理
    因 ...

第一次知道论坛里竟然有这么多好东西

我一直好奇那些神奇的数字是怎么来的
我猜到了求逆然后除法变乘法
然后想着，如果整除，我们可以用求得的mod 2^32的逆乘以被除数，得到商
然后不整除就没办法了

原来办法这么简单