How to calculate (time b) divided by c only using 32-bit integer types, even if time b is not suitable for this type

Question

As a reference implementation, consider the following:

/* 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; } 

I'm interested in an implementation (in C or pseudo-code) that does not require a 64-bit integer type.

I began to outline the implementation, which is described as follows:

 /* 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); } 

But the difficulty is the choice of values ​​for d1 and d2, which manage to avoid overflow ((a / d1) * (b / d2) <= UINT32_MAX) and minimize the error of the whole calculation.

Any thoughts?

Answer 1

I adapted the algorithm published by Paul for unsigned ints (by excluding parts that deal with characters). The algorithm is basically Ancient Egyptian multiplication a with a fraction of floor(b/c) + (b%c)/c (with a slash indicating the actual division here).

 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; } 

This algorithm will give an exact answer if it matches 32 bits. You can also return the remainder r .

Answer 2

A search on www.google.com/codesearch includes a number of implementations, including this surprisingly obvious one. I especially like the extensive comments and well-chosen variable names

 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

Answer 3

The easiest way is to convert the result of the intermediate result to 64 bits, but depending on the value of c, you can use a different approach:

 ((a/c)*b + (a%c)*(b/c) + ((a%c)*(b%c))/c 

The only problem is that the last member can still overflow for large c values. still thinking about it.

Answer 4

You can first divide a by c, and get a division reminder and multiply the reminder by b before dividing it by c. Thus, you only lose data in the last section, and you get the same result as 64-bit division.

You can rewrite the formula as follows (where \ is an integer division):

 a * b / c = (a / c) * b = (a \ c + (a % c) / c) * b = (a \ c) * b + ((a % c) * b) / c 

After making sure that a> = b, you can use large values ​​before overflowing:

 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; } 

Another approach would be the addition and subtraction of cycles instead of multiplication and division, but this, of course, works a lot more:

 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; } 

Answer 5

If b and c are both constants, you can simply calculate the result using Egyptian fractions.

For instance. y = a * 4/99 can be written as

 y = a / 25 + a / 2475 

You can express any fraction as the sum of the Egyptian fractions, as explained in the answers to the Egyptian fractions in C.

Having b and c fixed in advance may seem a bit limited, but this method is much simpler than the general case that others answer.

Answer 6

If b = 3000000000 => qn = 3000000000, qn * 2 will be full. Therefore, I am editing the Sven Marnach code.

 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; 

}

Answer 7

I believe there are reasons why you cannot do

 x = a/c; x = x*b; 

there is? And maybe add

 y = b/c; y = y*a; if ( x != y ) return ERROR_VALUE; 

Note that since you are using integer division, a*b/c and a/c*b can lead to different values ​​if c greater than a or b . Also, if both a and b less than c , this will not work.