c - 当 a 和 b 都小于 c，但 a * b 溢出时，如何计算 a * b/c？

Question

假设 uint是我定点平台上最大的积分类型，我有:

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

这需要返回 a * b / c 的一个很好的近似值.c的值大于 a 的值和 b 的值.
所以我们肯定知道 a * b / c 的值适合 uint .
但是，a * b 的值本身溢出了 uint 的大小.
所以一种计算 a * b / c 值的方法将是:

return a / c * b;

甚至:

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

但是，c 的值大于 a 的值和 b 的值.
所以上面的建议只会返回零。
我需要减少 a * b和 c按比例，但同样 - 问题是 a * b溢出。
理想情况下，我将能够:

替换 a * b与 uint(-1)

替换 c与 uint(-1) / a / b * c .

但无论我如何订购表达式 uint(-1) / a / b * c ，我遇到了一个问题:

uint(-1) / a / b * c由于 uint(-1) / a / b 被截断为零

uint(-1) / a * c / b溢出，因为 uint(-1) / a * c

uint(-1) * c / a / b溢出，因为 uint(-1) * c

我如何处理这种情况才能找到 a * b / c 的一个很好的近似值?

编辑 1
我没有_umul128之类的东西在我的平台上，当最大的整数类型是uint64 .我最大的类型是 uint ，并且我不支持比这更大的任何内容(无论是在硬件级别，还是在某些预先存在的标准库中)。
我最大的类型是 uint .
编辑 2
针对众多重复的建议和评论:
我手头没有一些“更大的类型”，我可以用它来解决这个问题。这就是为什么这个问题的开场白是:

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

我假设不存在其他类型，无论是在 SW 层(通过一些内置标准库)还是在 HW 层。

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

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.

使用 n 位 uint ，我希望结果至少对 n/2 有效数字是正确的。
可以通过利用 a,b 小于 SQRT_MAX_P1 的较小值来改进事情。如果有兴趣，稍后再详细介绍。

例子

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

用 32 位数学输出(准确答案是 12345677)

12345600  

64 位数学输出

12345677