how to calculate (a times b) divided by c only using 32-bit integer types even if a times b would not fit such a type

Question

Consider the following as a reference implementation:

/* 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 am interested in an implementation (in C or pseudocode) that does not require a 64-bit integer type.

I started sketching an implementation that outlines like this:

/* 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 to pick values for d1 and d2 that manage to avoid the overflow ((a / d1) * (b / d2) <= UINT32_MAX) and minimize the error of the whole calculation.

Any thoughts?

Answer 1

I have adapted the algorithm posted by Paul for unsigned ints (by omitting the parts that are dealing with signs). The algorithm is basically Ancient Egyptian multiplication of a with the fraction floor(b/c) + (b%c)/c (with the slash denoting real 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 yield the exact answer as long as it fits in 32 bits. You can optionally also return the remainder r.

Answer 2

The simplest way would be converting the intermediar result to 64 bits, but, depending on value of c, you could use another approach:

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

The only problem is that the last term could still overflow for large values of c. still thinking about it..

Answer 3

Searching on www.google.com/codesearch turns up a number of implementations, including this wonderfuly obvious one. I particularly 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 4

You can first divide a by c and also get the reminder of the division, and multiply the reminder with b before dividing it by c. That way you only lose data in the last division, and you get the same result as making the 64 bit division.

You can rewrite the formula like this (where \ is integer division):

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

By making sure that a >= b, you can use larger values before they overflow:

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 to loop addition and subtraction instead of multiplying and dividing, but that is of course a lot more work:

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 calculate the result very simply using Egyptian fractions.

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

y = a / 25 + a / 2475

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

Having b and c fixed in advance might seem like a bit of a restriction, but this method is a lot simpler than the general case answered by others.

Answer 6

If b = 3000000000 => qn = 3000000000, qn*2 will be overflowed. So I edit the code of 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;

}

Answer 7

I suppose there are reasons you can't do

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

are there? And maybe add

y = b/c;
y = y*a;

if ( x != y )
    return ERROR_VALUE;

Note that, since you're using integer division, a*b/c and a/c*b might lead to different values if c is bigger than a or b. Also, if both a and b are smaller than c it won't work.