c - 即使 a 乘以 b 不适合这种类型，如何仅使用 32 位整数类型计算 (a 乘以 b) 除以 c

Question

考虑以下作为参考实现：

/* 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;
}

我对不需要 64 位整数类型的实现（在 C 或伪代码中）感兴趣。

我开始草拟一个概述如下的实现：

/* 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);
}

但困难在于为 d1 和 d2 选择能够避免溢出的值 ((a / d1) * (b / d2) <= UINT32_MAX) 并最大限度地减少整个计算的错误。

有什么想法吗？

score 4 · Accepted Answer

我已经针对无符号整数调整了Paul发布的算法（通过省略处理符号的部分）。该算法基本上是古埃及a与分数的乘法（floor(b/c) + (b%c)/c这里的斜线表示真正的除法）。

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

只要它适合 32 位，该算法就会产生准确的答案。您还可以选择返回剩余部分r。

score 3 · Accepted Answer

您可以先将 a 除以 c 并获得除法的提示，然后将提示与 b 相乘，然后再除以 c。这样，您只会丢失最后一次除法中的数据，并且得到与进行 64 位除法相同的结果。

您可以像这样重写公式（其中 \ 是整数除法）：

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

通过确保 a >= b，您可以在溢出之前使用更大的值：

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

另一种方法是循环加减法而不是乘法和除法，但这当然需要更多的工作：

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

score 3 · Accepted Answer

在www.google.com/codesearch上搜索会发现许多实现，包括这个非常明显的实现。我特别喜欢广泛的评论和精心挑选的变量名

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

score 3 · Accepted Answer

最简单的方法是将中间结果转换为 64 位，但是，根据 c 的值，您可以使用另一种方法：

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

唯一的问题是，对于较大的值，最后一项仍然可能溢出c。还在想。。

score 1 · Accepted Answer

我将 Sven 的代码实现为 UINT16，以对其进行深入测试：

uint16_t muldiv16(uint16_t a, uint16_t b, uint16_t c);

int main(int argc, char *argv[]){
    uint32_t a;
    uint32_t b;
    uint32_t c;
    uint16_t r1, r2;

// ~167 days, estimated on i7 6700k, single thread.
// Split the 'a' range, to run several instances of this code on multi-cores processor
// ~1s, with an UINT8 implementation
    for(a=0; a<=UINT16_MAX; a++){
        for(b=0; b<=UINT16_MAX; b++){
            for(c=1; c<=UINT16_MAX; c++){
                r1 = uint16_t( a*b/c );
                r2 = muldiv16(uint16_t(a), uint16_t(b), uint16_t(c));
                if( r1 != r2 ){
                    std::cout << "Err: " << a << " * " << b << " / " << c << ", result: " << r2 << ", exected: " << r1 << std::endl;
                    return -1;
                }
            }
        }
        std::cout << a << std::endl
    }
    std::cout << "Done." << std::endl;
    return 0;
}

不幸的是，它似乎仅限于“b”（0-2147483647）的 UINT31。

这是我的更正，这似乎有效（未在 UINT16 上完成测试，但运行了很多。在 UINT8 上完成）。

uint32_t muldiv32(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;
    uint32_t r_carry;
    uint32_t rn_carry;
    while(a)
    {
        if (a & 1)
        {
            q += qn;
            r_carry = (r > UINT32_MAX-rn);
            r += rn;
            if (r >= c || r_carry)
            {
                q++;
                r -= c;
            }
        }
        a  >>= 1;
        qn <<= 1;
        rn_carry = rn & 0x80000000UL;
        rn <<= 1;
        if (rn >= c || rn_carry)
        {
            qn++;
            rn -= c;
        }
    }
    return q;
}

编辑：一项改进，返回余数，管理回合，警告溢出，当然，管理 a、b 和 c 的 UINT32 的全部范围：

typedef enum{
    ROUND_DOWNWARD=0,
    ROUND_TONEAREST,
    ROUND_UPWARD
}ROUND;

//remainder is always positive for ROUND_DOWN ( a * b = c * q + remainder )
//remainder is always negative for ROUND_UPWARD ( a * b = c * q - remainder )
//remainder is signed for ROUND_CLOSEST ( a * b = c * q + sint32_t(remainder) )
uint32_t muldiv32(uint32_t a, uint32_t b, uint32_t c, uint32_t *remainder, ROUND round, uint8_t *ovf)
{
    uint32_t q = 0;              // the quotient
    uint32_t r = 0;              // the remainder
    uint32_t qn = b / c;
    uint32_t rn = b % c;
    uint32_t r_carry;
    uint32_t rn_carry;
    uint8_t o = 0;
    uint8_t rup;
    while(a)
    {
        if (a & 1)
        {
            o |= (q > UINT32_MAX-qn);
            q += qn;
            r_carry = (r > UINT32_MAX-rn);
            r += rn;
            if (r >= c || r_carry)
            {
                o |= (q == UINT32_MAX);
                q++;
                r -= c;
            }
        }
        a  >>= 1;
        qn <<= 1;
        rn_carry = rn & 0x80000000;
        rn <<= 1;
        if (rn >= c || rn_carry)
        {
            qn++;
            rn -= c;
        }
    }
    rup = (round == ROUND_UPWARD && r);
    rup |= (round == ROUND_TONEAREST && ((r<<1) >= c || r & 0x80000000));
    if(rup)
    {   //round
        o |= (q == UINT32_MAX);
        q++;
        r = (round == ROUND_UPWARD) ? c-r : r-c;
    }
    if(remainder)
        *remainder = r;
    if(ovf)
        *ovf = o;
    return q;
}

也许可能存在另一种方法，也许效率更高：8 位、16 位和 32 位 MCU 能够计算 64 位计算（long long int）。任何人都知道编译器如何模拟它？

编辑2：

下面是一些有趣的时序，在 8 位 MCU 上：

UINT8 x UINT8 / UINT8：3.5µs

UINT16 x UINT16 / UINT16：22.5µs，muldiv8：29.9 至 45.3µs

UINT32 x UINT32 / UINT32：84µs，multidiv16：120 至 189µs

FLOAT32 * FLOAT32 / FLOAT32：40.2 到 135.5µs，muldiv32：1.193 到 1.764ms

在 32 位 MCU 上：

类型 - 优化代码 - 未优化

UINT32：521ns - 604ns

UINT64：2958ns - 3313ns

FLOAT32：2563ns - 2688ns

muldiv32：6791ns - 25375ns

所以，编译器比这个 C 算法聪明。并且使用浮点变量（即使没有 FPU）总是比使用大于本机寄存器的整数更好（即使浮点 32 的精度比 uint32 最差，从 16777217 开始）。

Edit3：好的，所以：我的 N 位 MCU 正在使用N-bits MUL N-bits本机指令，产生 2N 位结果，存储到两个 N 位寄存器中。

在这里，您可以找到一个 C 实现（更喜欢 EasyasPi 的解决方案）

但他们没有2N-bits DIV N-bits本地指令。相反，他们使用来自 gcc 的__udivdi3函数，带有循环和 2N 位变量（此处为 UINT64）。因此，这不能解决原始问题。

score 0 · Accepted Answer

如果 b 和 c 都是常数，您可以使用埃及分数非常简单地计算结果。

例如。y = a * 4 / 99 可以写成

y = a / 25 + a / 2475

您可以将任何分数表示为埃及分数的总和，如对C 中的埃及分数的回答中所述。

预先固定 b 和 c 似乎有点限制，但这种方法比其他人回答的一般情况要简单得多。

score -4 · Accepted Answer

我想有你不能做的理由

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

在那儿？也许添加

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

if ( x != y )
    return ERROR_VALUE;

请注意，由于您使用的是整数除法，a*b/c因此如果大于ora/c*b可能会导致不同的值。此外，如果两者和都小于它将不起作用。cab abc

c - 即使 a 乘以 b 不适合这种类型，如何仅使用 32 位整数类型计算 (a 乘以 b) 除以 c

7 回答 7

Related

Reference