以乘法实现除法--缺失小数精度

先说明一下这里要使用的一些运算符号：
'>>'  逻辑右移
'<<'  逻辑左移
└x┘   对x向负无穷方向取整
┌x┐   对x向正无穷方向取整
{x}   x的小数部分
x^y   x的y次幂  或者  是x异或y(程序表达式中使用)
mod(x,y)  x/y的余数
log2(x)  取x以2为底的对数

uint32  定义unsigned int无符号32bit整数
uint64  定义unsigned long long无符号64bit整数

在这里我们假设a为被除数，b为除数，且b > 1
我们知道
a / b
= (1 / b) * a
因为b>1，所以1/b是小于1的，而CPU通用寄存器只能处理整数，所以我们必须想办法把(1/b)扩大到2^e倍，也就是左移e个bit，然后把得到的结果在右移e个bit即可得到商：
└a/b┘
= └((2^e / b) * a) / 2^e┘
这里e是一个大于0的整数，具体由b来确定。
因为当b不是2的整数次幂时2^e / b不为整数，所以我们设想：
(2^e+c)/ b，这样得到一个整数，
所以当b不为2的整数次幂是c就等于：
c = b - mod(2^e,b)
b为2的整数次幂时:
c = 0
那么必须满足：
└a/b┘
=└((2^e / b) * a) / 2^e┘
=└((2^e+c)/ b) * a / 2^e┘
=└a / b + c*a/(b*2^e)┘
很明显只有当：
c*a/(b*2^e) < 1/b
=>
c*a < 2^e 时有：
└a / b + c*a/(b*2^e)┘
=└a/b┘
例：
当：
b = 3 时
如果e = 32
c = 3 - mod(2^32,3)
= 2
此时a为最大值2^32-1时，不满足：
c*a < 2^e
如果e = 33
c = 3 - mod(2^33,3)
= 1
此时满足：
c*a < 2^e
并且因为3 > 2，所以当e = 33，不会使
((2^e+c)/b) * a
超出 64bit而溢出，所以e = 33成立；
C语言可描述为：
uint32 a,b,r;
a = 被除数;
b = 0xAAAAAAAB; // b = (2^33 + 1) / 3
r = ((uint64)a * b) >> 33; //此时r中就是a/3的商

当：
b = 5 时
可以使用e = 34；

当：
b = 7 时
e = 32,33,34时都无法满足
c*a < 2^e
如果e = 32，那么
c = 3
为使c*a < 2^e
a最大就只能为0x55555555
当超出时就有可能使结果比实际商大1，
但如果让a加上这个商试试，

C语言描述为：

uint32 a,b,t,r;
a = 被除数;
b = 0x24924925; // b = (2^32 + 3) / 7
t = ((uint64)a * b) >> 32;
r = (a + t) >> 3; //如果(a + t)不溢出，这将是结果
因为a >= 7*T,T为实际商，得：
7*T <= a <= 7*T+6
=>
7*T+t <= a+t <= 7*T+6+t
=>
8*T <= a+t <= 8*T+6+1 < 8*(T+1)
=>
└(a+t)/8┘ = T
所以(a + t) >> 3含义成立，
但(a + t)有可能产生溢出，改进一下得：
(a + t) >> 3
= ((a + t) >> 1) >> 2
= └(a + t) / 2┘ >> 2
= └(a-t + 2t) / 2┘ >> 2
= (└(a-t)/2┘+t) >> 2
= (((a-t) >> 1) + t) >> 2
这样就可以不担心溢出了，最后
r = (((a-t) >> 1) + t) >> 2;
看了上面几个数的除法实现，应该有了一个初步的了解，现在我们来构造一个通用的算法来求其e，知道了合适的e，一切都好办了。
一些数字可以对e满足((2^e+c)/ b) * a不超过64bit时有：
c*a < 2^e
但还有一些数则无法办到，所以我们退其次，我们让((2^e+c)/ b) * a可以为65bit，这样我们将永远满足：
c*a < 2^e
因为：
a*(2^e+c)/b <= 2^64-1
=>
(2^e+c)/b <= (2^64-1)/(2^32-1)
=>
(2^e+c)/b <= 2^32+1
当c=0时：
2^e / b为整数，
=>
2^e / b <= 2^32+1
=>
e <= └log2(b) + log2(2^32+1)┘
=>
e <= log2(b) + 32
当c!=0时：
=>
2^e / b < 2^32+1
=>
2^e < b*(2^32+1)
=>
e <= ┌log2(b) + log2(2^32+1)┐ - 1
=>
e <= └log2(b) + log2(2^32+1)┘
=>
e <= └log2(b)┘ + 32
结合两次分析最终得：
e <= └log2(b)┘ + 32
时
((2^e+c)/ b) * a
不超过64bit，
现在我们使用：
e = └log2(b)┘ + 33来完成
((2^e+c)/ b) * a
这样此式的结果最高就为65bit了，((2^e+c)/ b)就是33bit。

但((2^e+c)/ b)在一个寄存器中最多只能表示32bit的数，还有一个最高bit暂时忽略，只要其低32bit，然后我们在结果的高32上加一个a，这样我们就得到了正确结果的低64bit，C语言描述为：

uint32 a,b,t,r;
a = 被除数;
b = 0x********;//b = (2^(e+33)+c)/b, e = └log2(b)┘
t = ((uint64)a * (uint64)b) >> 32;
r = (a + t) >> (e+1);//商
刚刚已说(a + t)可能溢出，换一种表达式来表示：
(a + t) >> 1
= └(a + t) / 2┘
= └(a-t + 2t) / 2┘
= └(a-t)/2┘+t
= (a - t) >> 1 + t
最终r = ((a - t) >> 1 + t) >> e
这样就不会担心中间结果溢出了。
下面是为上述而设计的一个C语言程序：

typedef unsigned int uint32_t;typedef int int32_t;typedef unsigned long long uint64_t;#define is2exp(x) (((x - 1) & x) == 0)//检测x是否为2的整数次幂uint32_t log2x(uint32_t x){//└log2(x)┘ , x>0//求小于等于x以2为底的对数的最大整数 int32_t e; for (e = 31; e >= 0; e--) {  if ((int32_t)x < 0) return e;  x <<= 1; } return -1;}int _tmain(int argc, _TCHAR* argv[]){ uint32_t b,c,i; b = 21;//除数，b > 1 for (i = 0; true; i++) {  c = b - ((uint64_t)1 << (32+i)) % b;  if (c == b) {   c = 0;//此时b为2的整数次幂   break;  } else if (is2exp(c)) {   if (log2x(c) <= i) break;  } else {   if (log2x(c) < i) break;  } } uint32_t e,db; e = log2x(b); db = (((uint64_t)1 << (32+i)) + c) / b;//当i = e+1时取结果的低32bit printf("当除数为 %d 时:/n",b); if (i <= e) {  printf("商 = ((uint64_t)被除数 * 0x%x) >> %d",db,32+i); } else {//i = e + 1时就说明在 e <= └log2(b)┘内无法完全满足c*a < 2^(e+32)  printf("t = ((uint64_t)被除数 * 0x%x) >> 32/n",db);  printf("商 = ((((uint64_t)被除数 - t) >> 1) + t) >> %d",e); } getchar(); return 0;}

测试：
如果b = 21将输出：
当除数为 21 时:
t = ((uint64_t)被除数 * 0x86186187) >> 32

商 = ((((uint64_t)被除数 - t) >> 1) + t) >> 4

注：记不得是哪里的文章了，经过整理。