Ho bisogno di un veloce 96-bit su algoritmo di divisione specifico a 64 bit per una libreria matematica a virgola fissa

Question

Attualmente sto scrivendo una libreria matematica rapida a 32.32 punti fissi. Sono riuscito a fare correttamente l'aggiunta, la sottrazione e la moltiplicazione, ma sono abbastanza bloccato alla divisione.

Un piccolo promemoria per coloro che non riescono a ricordare: un numero a 32.32 in virgola fissa è un numero con 32 bit di parte intera e 32 bit di parte frazionaria.

Il miglior algoritmo che ho trovato ha bisogno della divisione integer a 96 bit, che è qualcosa per cui i compilatori di solito non hanno built-in.

In ogni caso, qui va:

G = 2^32 

notation: x is the 64-bit fixed-point number, x1 is its low nibble and x2 is its high 

G*(a/b) = ((a1 + a2*G)/(b1 + b2*G))*G  // Decompose this 

G*(a/b) = (a1*G)/(b1*G + b2) + (a2*G*G)/(b1*G + b2) 

Come si può vedere, la (a2*G*G) è garantito per essere più grande del normale a 64 bit integer. Se uint128_t di stati effettivamente sostenuti dal mio compilatore, vorrei semplicemente fare quanto segue:

((uint128_t)x << 32)/y) 

Beh non lo sono e ho bisogno di una soluzione. Grazie per l'aiuto.

Answer 1

migliore auto-regolazione risposta:
perdonare il C# -ismo della risposta, ma il seguente dovrebbe funzionare in tutti i casi. È probabile che sia possibile una soluzione che trovi i giusti cambiamenti da usare più velocemente, ma dovrei pensare molto più in profondità di quanto possa fare adesso. Questo dovrebbe essere ragionevolmente efficiente, anche se:

int upshift = 32; 
ulong mask = 0xFFFFFFFF00000000; 
ulong mod = x % y; 
while ((mod & mask) != 0) 
{ 
    // Current upshift of the remainder would overflow... so adjust 
    y >>= 1; 
    mask <<= 1; 
    upshift--; 

    mod = x % y; 
} 
ulong div = ((x/y) << upshift) + (mod << upshift)/y; 

Risposta semplice ma non sicuro:
Questo calcolo può causare un overflow nel cambio marcia del x % y restante se questo resto ha dei bit impostati in alta 32 bit , causando una risposta errata.

((x/y) << 32) + ((x % y) << 32)/y 

La prima parte utilizza divisione intera e si dà le alte bit della risposta (spostarli back up).

La seconda parte calcola i bit bassi dal resto della divisione dei bit alti (il bit che non è stato possibile dividere ulteriormente), spostati verso l'alto e quindi divisi.

Answer 2

È possibile scomporre una divisione più grande in più blocchi che eseguono la divisione con meno bit. Come già menzionato un altro poster, l'algoritmo può essere trovato in TAOCP da Knuth.

Tuttavia, non è necessario acquistare il libro!

C'è un codice nel sito Web di hacker che implementa l'algoritmo in C. È scritto per fare divisioni a 64 bit senza segno utilizzando solo aritmetica a 32 bit, quindi non è possibile tagliare il codice direttamente. Per ottenere da 64 a 128 bit devi allargare tutti i tipi, maschere e costanti di due, ad es. un breve diventa un int, un 0xffff diventa 0xffffffffll ect.

Dopo questo facile cambiamento, dovresti essere in grado di fare divisioni a 128 bit.

Il codice è qui: http://www.hackersdelight.org/HDcode/divlu.c (potrebbe terminare male in un browser Web a causa di terminazioni di riga. In tal caso, è sufficiente salvare il codice e aprirlo con il blocco note o così via).

Poiché i valori massimi richiedono solo 96 bit, una delle divisioni a 64 bit sarà sempre a zero, in modo che sia possibile semplificare un po 'il codice.

Oh - e prima che lo dimentichi: il codice funziona solo con valori senza segno. Per convertire da firmato per dividere senza segno che si può fare qualcosa di simile (stile pseudo-codice):

fixpoint Divide (fixpoint a, fixpoint b) 
{ 
    // check if the integers are of different sign: 
    fixpoint sign_difference = a^b; 

    // do unsigned division: 
    fixpoint x = unsigned_divide (abs(a), abs(b)); 

    // if the signs have been different: negate the result. 
    if (sign_difference < 0) 
    { 
    x = -x; 
    } 

    return x; 
} 

Il sito è di per sé la pena di verificare così: http://www.hackersdelight.org/

Speranza che aiuta.

Btw - bel compito su cui stai lavorando. Ti dispiace dirci per cosa hai bisogno della libreria di punti fissi?

---

Btw - anche l'algoritmo di sottrazione e sottrazione ordinario per divisione funzionerebbe.

Se si target x86, è possibile implementarlo utilizzando gli intrinsechi MMX o SSE. L'algoritmo si basa solo su operazioni primitive, quindi potrebbe eseguire anche abbastanza velocemente.

Answer 3

Mi piace la risposta di Nils, che è probabilmente la migliore. E 'solo la divisione lungo, come abbiamo tutti imparato in grado di scuola, tranne le cifre sono di base 2^32, invece di base 10.

Tuttavia, si potrebbe anche considerare l'utilizzo metodo di approssimazione di Newton per la divisione:

x := x (N + N - N * D * x) 

dove N è il numeratore e D è il demoninatore.

Questo utilizza solo moltiplica e aggiunge, che hai già, e converge molto rapidamente a circa 1 ULP di precisione. D'altra parte, non sarà in grado di ottenere la risposta esatta 0,5-ULP in tutti i casi.

In ogni caso, la parte più difficile è il rilevamento e la gestione degli overflow.

Answer 4

Quick -n- sporco.

Fare il divario A/B con il doppio punto a virgola mobile di precisione. Questo ti dà C ~ = A/B. È approssimativamente a causa della precisione in virgola mobile e 53 bit di mantissa.

Arrotondare C a un numero rappresentabile nel sistema a punto fisso.

Calcolare ora (sempre con il punto fisso) D = A-C * B. Questo dovrebbe avere un'ampiezza significativamente inferiore a quella di A.

Ripetere, ora calcolare D/B con virgola mobile. Ancora una volta, arrotondare la risposta a un numero intero. Aggiungi ogni risultato di divisione insieme mentre procedi. Puoi fermarti quando il tuo resto è così piccolo che la tua divisione in virgola mobile restituisce 0 dopo l'arrotondamento.

Non hai ancora finito. Ora sei molto vicino alla risposta, ma le divisioni non erano esatte. Per finalizzare, dovrai fare una ricerca binaria. Usando la (molto buona) stima iniziale, vedi se aumentarla migliora l'errore .. fondamentalmente vuoi raggruppare la risposta corretta e continuare a dividere la gamma a metà con nuovi test.

Sì, è possibile eseguire l'iterazione di Newton qui, ma la ricerca binaria sarà probabilmente più semplice poiché è necessario solo moltiplicare e aggiungere semplici utilizzando il toolkit di precisione 32.32 esistente.

Questo è non il metodo più efficiente, ma è di gran lunga il più semplice da codificare.