如何计算(a次b)除以c仅使用32位整数类型,即使b次不适合这种类型
请考虑以下内容作为参考实现:
/* calculates (a * b) / c */ uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) { uint64_t x = a; x = x * b; x = x / c; return x; }
我感兴趣的是一个不需要64位整数类型的实现(在C或伪代码中)。
我开始草拟一个如下概述的实现:
/* calculates (a * b) / c */ uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) { uint32_t d1, d2, d1d2; d1 = (1 << 10); d2 = (1 << 10); d1d2 = (1 << 20); /* d1 * d2 */ return ((a / d1) * (b /d2)) / (c / d1d2); }
但困难在于选择d1和d2的值来避免溢出((a / d1)*(b / d2)<= UINT32_MAX)并最小化整个计算的误差。
有什么想法吗?
我已经调整了Paul发布的无符号整数算法(通过省略处理符号的部分)。 该算法基本上是古埃及乘以 a
与分数floor(b/c) + (b%c)/c
(斜线表示实际分割)。
uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) { uint32_t q = 0; // the quotient uint32_t r = 0; // the remainder uint32_t qn = b / c; uint32_t rn = b % c; while(a) { if (a & 1) { q += qn; r += rn; if (r >= c) { q++; r -= c; } } a >>= 1; qn <<= 1; rn <<= 1; if (rn >= c) { qn++; rn -= c; } } return q; }
只要它适合32位,该算法将产生确切的答案。 您也可以选择返回余数r
。
最简单的方法是将intermediar结果转换为64位,但是,根据c的值,您可以使用另一种方法:
((a/c)*b + (a%c)*(b/c) + ((a%c)*(b%c))/c
唯一的问题是,对于较大的c
值,最后一个术语仍然可能溢出。 还在考虑它..
在www.google.com/codesearch上进行搜索会产生许多实现,包括这个非常明显的实现。 我特别喜欢广泛的评论和精心选择的变量名称
INT32 muldiv(INT32 a, INT32 b, INT32 c) { INT32 q=0, r=0, qn, rn; int qneg=0, rneg=0; if (c==0) c=1; if (a<0) { qneg=!qneg; rneg=!rneg; a = -a; } if (b<0) { qneg=!qneg; rneg=!rneg; b = -b; } if (c<0) { qneg=!qneg; c = -c; } qn = b / c; rn = b % c; while(a) { if (a&1) { q += qn; r += rn; if(r>=c) { q++; r -= c; } } a >>= 1; qn <<= 1; rn <<= 1; if (rn>=c) {qn++; rn -= c; } } result2 = rneg ? -r : r; return qneg ? -q : q; }
http://www.google.com/codesearch/p?hl=en#HTrPUplLEaU/users/mr/MCPL/mcpl.tgz|gIE-sNMlwIs/MCPL/mintcode/sysc/mintsys.c&q=muldiv%20lang:c
您可以先将c除以并得到除法的提示,然后将提示与b相乘,然后再除以c。 这样你只会在最后一个分区中丢失数据,并且得到与进行64位除法相同的结果。
您可以像这样重写公式(其中\是整数除法):
a * b / c = (a / c) * b = (a \ c + (a % c) / c) * b = (a \ c) * b + ((a % c) * b) / c
通过确保a> = b,您可以在溢出之前使用更大的值:
uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) { uint32_t hi = a > b ? a : b; uint32_t lo = a > b ? b : a; return (hi / c) * lo + (hi % c) * lo / c; }
另一种方法是循环加法和减法而不是乘法和除法,但这当然是更多的工作:
uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) { uint32_t hi = a > b ? a : b; uint32_t lo = a > b ? b : a; uint32_t sum = 0; uint32_t cnt = 0; for (uint32_t i = 0; i < hi; i++) { sum += lo; while (sum >= c) { sum -= c; cnt++; } } return cnt; }
如果b和c都是常数,则可以使用Egyptian分数非常简单地计算结果。
例如。 y = a * 4/99可以写成
y = a / 25 + a / 2475
您可以将任何分数表示为埃及分数的总和,如C中埃及分数的答案中所述 。
事先确定b和c可能看起来有点限制,但这种方法比其他人回答的一般情况简单得多。
如果b = 3000000000 => qn = 3000000000,则qn * 2将溢出。 所以我编辑了Sven Marnach的代码。
uint32_t muldiv(uint32_t a, uint32_t b, uint32_t c) { uint32_t q = 0; // the quotient uint32_t r = 0; // the remainder uint32_t qn = b / c; uint32_t rn = b % c; while (a) { if (a & 1) { q += qn; if (qn >= UINT32_MAX) { cout << "CO CO" << endl; } r += rn; if (r >= c) { q++; r -= c; } } a >>= 1; qn <<= 1; int temp = rn; if (rn > INT32_MAX) { // rn times 2: overflow rn = UINT32_MAX;// rn temp = (temp - INT32_MAX) * 2; // find the compensator mean: rn * 2 = UINT32_MAX + temp qn++; rn = rn - c + temp; } else { rn <<= 1; if (rn >= c) { qn++; rn -= c; } } } //return r; return q;
}
我想你有理由不能这样做
x = a/c; x = x*b;
在那儿? 也许可以添加
y = b/c; y = y*a; if ( x != y ) return ERROR_VALUE;
请注意,由于您使用整数除法,如果c
大于a
或b
,则a*b/c
和a/c*b
可能会导致不同的值。 此外,如果a
和b
都小于c
,它将无法工作。