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

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-битного целочисленного типа.

Я начал набрасывать реализацию, которая обрисовывает в общих чертах, как это:

/* 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, которым удается избежать переполнения ((а/d1) * (б/d2) < = UINT32_MAX) и минимизировать ошибку всего расчета.

Любые мысли?

Answer 1

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

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

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

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

Answer 3

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

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, это не сработает.

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 

Убедившись, что> = Ь, вы можете использовать более высокие значения, прежде чем они переполнения:

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

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

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

Например. у = А * 4/99 может быть записана в виде

y = a/25 + a/2475 

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

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

Answer 7

Если 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; 

}