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

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