How can I compute a * b / c when both a and b are smaller than c, but a * b overflows?

Question

Assuming that uint is the largest integral type on my fixed-point platform, I have:

uint func(uint a, uint b, uint c);

Which needs to return a good approximation of a * b / c.

The value of c is greater than both the value of a and the value of b.

So we know for sure that the value of a * b / c would fit in a uint.

However, the value of a * b itself overflows the size of a uint.

So one way to compute the value of a * b / c would be:

return a / c * b;

Or even:

if (a > b)
    return a / c * b;
return b / c * a;

However, the value of c is greater than both the value of a and the value of b.

So the suggestion above would simply return zero.

I need to reduce a * b and c proportionally, but again - the problem is that a * b overflows.

Ideally, I would be able to:

• Replace a * b with uint(-1)
• Replace c with uint(-1) / a / b * c.

But no matter how I order the expression uint(-1) / a / b * c, I encounter a problem:

• uint(-1) / a / b * c is truncated to zero because of uint(-1) / a / b
• uint(-1) / a * c / b overflows because of uint(-1) / a * c
• uint(-1) * c / a / b overflows because of uint(-1) * c

How can I tackle this scenario in order to find a good approximation of a * b / c?

---

Edit 1

I do not have things such as _umul128 on my platform, when the largest integral type is uint64. My largest type is uint, and I have no support for anything larger than that (neither on the HW level, nor in some pre-existing standard library).

My largest type is uint.

Edit 2

In response to numerous duplicate suggestions and comments:

I do not have some "larger type" at hand, which I can use for solving this problem. That is why the opening statement of the question is:

Assuming that uint is the largest integral type on my fixed-point platform

I am assuming that no other type exists, neither on the SW layer (via some built-in standard library) nor on the HW layer.

Answer 1

needs to return a good approximation of a * b / c
My largest type is uint
both a and b are smaller than c

Variation on this 32-bit problem:

Algorithm: Scale a, b to not overflow

SQRT_MAX_P1 as a compile time constant of sqrt(uint_MAX + 1)
sh = 0;
if (c >= SQRT_MAX_P1) {
  while (|a| >= SQRT_MAX_P1) a/=2, sh++
  while (|b| >= SQRT_MAX_P1) b/=2, sh++
  while (|c| >= SQRT_MAX_P1) c/=2, sh--
}
result = a*b/c

shift result by sh.

With an n-bit uint, I expect the result to be correct to at least about n/2 significant digits.

Could improve things by taking advantage of the smaller of a,b being less than SQRT_MAX_P1. More on that later if interested.

---

Example

#include <inttypes.h>

#define IMAX_BITS(m) ((m)/((m)%255+1) / 255%255*8 + 7-86/((m)%255+12))
// https://stackoverflow.com/a/4589384/2410359

#define UINTMAX_WIDTH (IMAX_BITS(UINTMAX_MAX))
#define SQRT_UINTMAX_P1 (((uintmax_t)1ull) << (UINTMAX_WIDTH/2))

uintmax_t muldiv_about(uintmax_t a, uintmax_t b, uintmax_t c) {
  int shift = 0;
  if (c > SQRT_UINTMAX_P1) {
    while (a >= SQRT_UINTMAX_P1) {
      a /= 2; shift++;
    }
    while (b >= SQRT_UINTMAX_P1) {
      b /= 2; shift++;
    }
    while (c >= SQRT_UINTMAX_P1) {
      c /= 2; shift--;
    }
  }
  uintmax_t r = a * b / c;
  if (shift > 0) r <<= shift;
  if (shift < 0) r >>= shift;
  return r;
}



#include <stdio.h>

int main() {
  uintmax_t a = 12345678;
  uintmax_t b = 4235266395;
  uintmax_t c = 4235266396;
  uintmax_t r = muldiv_about(a,b,c);
  printf("%ju\n", r);
}

Output with 32-bit math (Precise answer is 12345677)

12345600  

Output with 64-bit math

12345677  

Answer 2

Here is another approach that uses recursion and minimal approximation to achieve high precision.

First the code and below an explanation.

Code:

uint32_t bp(uint32_t a) {
  uint32_t b = 0;
  while (a!=0)
  {
    ++b;
    a >>= 1;
  };
  return b;
}

int mul_no_ovf(uint32_t a, uint32_t b)
{
  return ((bp(a) + bp(b)) <= 32);
}

uint32_t f(uint32_t a, uint32_t b, uint32_t c)
{
  if (mul_no_ovf(a, b))
  {
    return (a*b) / c;
  }

  uint32_t m = c / b;
  ++m;
  uint32_t x = m*b - c;
  // So m * b == c + x where x < b and m >= 2

  uint32_t n = a/m;
  uint32_t r = a % m;
  // So a*b == n * (c + x) + r*b == n*c + n*x + r*b where r*b < c

  // Approximation: get rid of the r*b part
  uint32_t res = n;
  if (r*b > c/2) ++res;

  return res + f(n, x, c);
}

Explanation:

The multiplication a * b can be written as a sum of b

a * b = b + b + .... + b

Since b < c we can take a number m of these b so that (m-1)*b < c <= m*b, like

(b + b + ... + b) + (b + b + ... + b) + .... + b + b + b
\---------------/   \---------------/ +        \-------/
       m*b        +        m*b        + .... +     r*b
     \-------------------------------------/
            n times m*b

so we have

a*b = n*m*b + r*b

where r*b < c and m*b > c. Consequently, m*b is equal to c + x, so we have

a*b = n*(c + x) + r*b = n*c + n*x + r*b

Divide by c :

a*b/c = (n*c + n*x + r*b)/c = n + n*x/c + r*b/c

The values m, n, x, r can all be calculated from a, b and c without any loss of 
precision using integer division (/) and remainder (%).

The approximation is to look at r*b (which is less than c) and "add zero" when r*b<=c/2
and "add one" when r*b>c/2.

So now there are two possibilities:

1) a*b = n + n*x/c

2) a*b = (n + 1) + n*x/c

So the problem (i.e. calculating a*b/c) has been changed to the form

MULDIV(a1,b1,c) = NUMBER + MULDIV(a2,b2,c)

where a2,b2 is less than a1,b2. Consequently, recursion can be used until 
a2*b2 no longer overflows (and the calculation can be done directly).