как рассчитать (a b), деленное на c, только с использованием 32-битных целочисленных типов, даже если a b не подходит для такого типа

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

Меня интересует реализация (в C или псевдокоде), которая не требует 64-битного целочисленного типа.

Iначал делать набросок реализации, которая выглядит следующим образом:

/* 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) и минимизировать ошибку всего расчета. </p>

Есть мысли?

Answer 1

Я адаптировал алгоритм, опубликованный Полом для неподписанных целых (исключая части, которые имеют дело со знаками).Алгоритм в основном умножение древнеегипетского из 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.

Answer 2

Поиск по 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

Answer 3

Простейшим способом будет преобразование промежуточного результата в 64 бита, но, в зависимости от значения c, вы можете использовать другой подход:

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

Единственная проблема заключается в том, что последний член все еще может переполниться для больших значений c. все еще думаю об этом ..

Answer 4

Вы можете сначала разделить 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;
}

Answer 5

Если b = 3000000000 => qn = 3000000000, qn * 2 будет переполнено.Поэтому я редактирую код Свена Марнача.

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;
        if (qn >= UINT32_MAX) {
            cout << "CO CO" << endl;
        }
        r += rn;
        if (r >= c)
        {
            q++;
            r -= c;
        }
    }
    a >>= 1;
    qn <<= 1;
    int temp = rn;
    if (rn > INT32_MAX) {
        // rn times 2: overflow
        rn = UINT32_MAX;// rn 
        temp = (temp - INT32_MAX) * 2; // find the compensator mean: rn * 2  = UINT32_MAX + temp
        qn++;
        rn = rn - c + temp;
    }
    else {
        rn <<= 1;
        if (rn >= c)
        {
            qn++;
            rn -= c;
        }
    }


}

//return r;
return q;

}

Answer 6

Если b и c обе константы, вы можете очень просто вычислить результат, используя египетские дроби.

Например. y = a * 4/99 можно записать как

y = a / 25 + a / 2475

Вы можете выразить любую дробь как сумму египетских дробей, как объяснено в ответах на Египетские дроби в C .

Предварительное исправление b и c может показаться ограничением, но этот метод намного проще, чем общий случай, на который отвечают другие.

Answer 7

Полагаю, есть причины, по которым вы не можете это сделать

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

есть?И, возможно, добавьте

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

if ( x != y )
    return ERROR_VALUE;

Обратите внимание: поскольку вы используете целочисленное деление, a*b/c и a/c*b могут привести к различным значениям, если c больше a или b.Кроме того, если и a, и b меньше, чем c, это не будет работать.