策划行业网站/怎样建网站
快速立即除法的乘法实现(通用算法)
原创：HAM
这里我们要讲的是如何使用整数乘法来完成除数为常数的整数除法运算，我们假定都是无符号操作，并且运算在Intel 32bit x86 CPU上。
在整数除法操作后，我们取的都是商和余数，于是我们猜想用乘法来完成这一操作，因为通常乘法操作比除法快得多。
在VC++编译器中已经专门为此做了很好的优化，但我依然执著于去研究一下。
先说明一下这里要使用的一些运算符号：
'>>' 逻辑右移
'<
└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语言程序：
#include "stdafx.h"
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