当 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将适合 a uint。

但是，a * b本身的值会溢出 a 的大小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

针对众多重复的建议和评论：

我手头没有一些“更大的类型”，我可以用它来解决这个问题。这就是为什么这个问题的开场白是：

假设这uint是我定点平台上最大的整数类型

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

Answer 1

我已经建立了一个复杂的解决方案O(1)（无循环）：

typedef unsigned long long uint;

typedef struct
{
    uint n;
    uint d;
}
fraction;

uint func(uint a, uint b, uint c);
fraction reducedRatio(uint n, uint d, uint max);
fraction normalizedRatio(uint a, uint b, uint scale);
fraction accurateRatio(uint a, uint b, uint scale);
fraction toFraction(uint n, uint d);
uint roundDiv(uint n, uint d);

uint func(uint a, uint b, uint c)
{
    uint hi = a > b ? a : b;
    uint lo = a < b ? a : b;
    fraction f = reducedRatio(hi, c, (uint)(-1) / lo);
    return f.n * lo / f.d;
}

fraction reducedRatio(uint n, uint d, uint max)
{
    fraction f = toFraction(n, d);
    if (n > max || d > max)
        f = normalizedRatio(n, d, max);
    if (f.n != f.d)
        return f;
    return toFraction(1, 1);
}

fraction normalizedRatio(uint a, uint b, uint scale)
{
    if (a <= b)
        return accurateRatio(a, b, scale);
    fraction f = accurateRatio(b, a, scale);
    return toFraction(f.d, f.n);
}

fraction accurateRatio(uint a, uint b, uint scale)
{
    uint maxVal = (uint)(-1) / scale;
    if (a > maxVal)
    {
        uint c = a / (maxVal + 1) + 1;
        a /= c; // we can now safely compute `a * scale`
        b /= c;
    }
    if (a != b)
    {
        uint n = a * scale;
        uint d = a + b; // can overflow
        if (d >= a) // no overflow in `a + b`
        {
            uint x = roundDiv(n, d); // we can now safely compute `scale - x`
            uint y = scale - x;
            return toFraction(x, y);
        }
        if (n < b - (b - a) / 2)
        {
            return toFraction(0, scale); // `a * scale < (a + b) / 2 < MAXUINT256 < a + b`
        }
        return toFraction(1, scale - 1); // `(a + b) / 2 < a * scale < MAXUINT256 < a + b`
    }
    return toFraction(scale / 2, scale / 2); // allow reduction to `(1, 1)` in the calling function
}

fraction toFraction(uint n, uint d)
{
    fraction f = {n, d};
    return f;
}

uint roundDiv(uint n, uint d)
{
    return n / d + n % d / (d - d / 2);
}

这是我的测试：

#include <stdio.h>

int main()
{
    uint a = (uint)(-1) / 3;            // 0x5555555555555555
    uint b = (uint)(-1) / 2;            // 0x7fffffffffffffff
    uint c = (uint)(-1) / 1;            // 0xffffffffffffffff
    printf("0x%llx", func(a, b, c));    // 0x2aaaaaaaaaaaaaaa
    return 0;
}