come calcolare (a volte b) diviso per c usando solo tipi interi a 32 bit anche se un tempo b non si adatta a un tipo

Considera quanto segue come un’implementazione di riferimento:

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

Sono interessato a un’implementazione (in C o pseudocodice) che non richiede un tipo intero a 64 bit.

Ho iniziato a delineare un’implementazione che delinea così:

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

Ma la difficoltà è quella di selezionare i valori per d1 e d2 che riescono a evitare l’overflow ((a / d1) * (b / d2) <= UINT32_MAX) e minimizzare l'errore dell'intero calcolo.

qualche idea?

Ho adattato l’algoritmo pubblicato da Paul per gli unsigned (omettendo le parti che hanno a che fare con i segni). L’algoritmo è fondamentalmente una moltiplicazione egiziana antica di a con il floor(b/c) + (b%c)/c della frazione floor(b/c) + (b%c)/c (con la barra che denota la vera divisione qui).

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

Questo algoritmo fornirà la risposta esatta finché si adatta a 32 bit. Opzionalmente puoi anche restituire il resto r .

Il modo più semplice sarebbe convertire il risultato intermedio in 64 bit, ma, in base al valore di c, potresti usare un altro approccio:

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

L’unico problema è che l’ultimo termine potrebbe ancora traboccare per grandi valori di c . ancora a pensarci ..

La ricerca su http://www.google.com/codesearch presenta una serie di implementazioni, inclusa questa meravigliosa e ovvia. Mi piacciono particolarmente i commenti estesi e i nomi delle variabili ben scelti

 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

Puoi dapprima dividere per c e ottenere anche il promemoria della divisione, e moltiplicare il promemoria con b prima di dividerlo per c. In questo modo perdi solo i dati nell’ultima divisione e ottieni lo stesso risultato della divisione a 64 bit.

Puoi riscrivere la formula come questa (dove \ è la divisione in interi):

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

Assicurandoti che a> = b, puoi utilizzare valori più grandi prima che si sovrappongano:

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

Un altro approccio potrebbe essere l’addizione e la sottrazione del ciclo invece di moltiplicare e dividere, ma questo è ovviamente molto più lavoro:

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

Se b e c sono entrambe le costanti, puoi calcolare il risultato molto semplicemente usando le frazioni egiziane.

Per esempio. y = a * 4/99 può essere scritto come

 y = a / 25 + a / 2475 

Puoi esprimere qualsiasi frazione come sum delle frazioni egiziane, come spiegato nelle risposte alle frazioni egiziane in C.

Avere b e c fissato in anticipo potrebbe sembrare un po ‘una restrizione, ma questo metodo è molto più semplice del caso generale a cui gli altri hanno risposto.

Se b = 3000000000 => qn = 3000000000, qn * 2 sarà in overflow. Quindi modifico il codice di Sven Marnach.

 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; 

}

Suppongo che ci siano dei motivi che non puoi fare

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

ci sono? E forse aggiungo

 y = b/c; y = y*a; if ( x != y ) return ERROR_VALUE; 

Nota che, dal momento che stai usando la divisione integer, a*b/c e a/c*b potrebbero portare a valori diversi se c è maggiore di a o b . Inoltre, se sia a che b sono più piccoli di c non funzionerà.