如何计算(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大于ab ,则a*b/ca/c*b可能会导致不同的值。 此外,如果ab都小于c ,它将无法工作。