Elaborazione di alti 64 bit di un prodotto int 64x64 in C

Question

Vorrei che la mia funzione C calcolasse in modo efficiente gli alti 64 bit del prodotto di due firmati a 64 bit. So come farlo nell'assemblaggio x86-64, con imulq e estraendo il risultato da% rdx. Ma non riesco a scriverlo affatto in C, e tanto meno convincere il compilatore a farlo in modo efficiente.

Qualcuno ha qualche suggerimento per scrivere questo in C? Questo è sensibile alle prestazioni, quindi i "metodi manuali" (come le librerie di Russian Beasant o bignum) sono fuori.

Questa funzione assembly inline dorky ho scritto funziona ed è grosso modo la codegen che sto cercando:

static long mull_hi(long inp1, long inp2) { 
    long output = -1; 
    __asm__("movq %[inp1], %%rax;" 
      "imulq %[inp2];" 
      "movq %%rdx, %[output];" 
      : [output] "=r" (output) 
      : [inp1] "r" (inp1), [inp2] "r" (inp2) 
      :"%rax", "%rdx"); 
    return output; 
} 

Answer 1

La risposta generale è che x * y possono essere suddivisi in (a + b) * (c + d), dove a e c sono nell'ordine alto parti.

In primo luogo, si espandono per ac + ad + bc + bd

Ora, si moltiplicano le condizioni dei numeri a 32 bit memorizzati come long long (o meglio ancora, uint64_t), e basta ricordare che quando si moltiplicato un numero d'ordine più alto, è necessario scala di 32 bit. Quindi fai gli add, ricordando di rilevare carry. Tieni traccia del segno. Naturalmente, devi fare gli addendi a pezzi.

Answer 2

Aspetta, hai già una soluzione di assemblaggio perfettamente funzionante e ottimizzata, , per questo motivo, e vuoi eseguire il back-out e provare a scrivere in un ambiente che non supporta la matematica a 128 bit? Non sto seguendo.

Come ovviamente sapete, questa operazione è una singola istruzione su x86-64. Ovviamente niente di ciò che farai funzionerà meglio. Se vuoi davvero C portatile, devi fare qualcosa come il codice DigitalRoss sopra e spero che il tuo ottimizzatore capisca cosa stai facendo .

Se avete bisogno di architettura portabilità, ma sono disposti a limitarsi alle piattaforme gcc, ci sono __int128_t (e __uint128_t) tipi nei intrinseci compilatore che farà ciò che si desidera.

Answer 3

se si sta utilizzando un relativamente recente GCC su x86_64:

int64_t mulHi(int64_t x, int64_t y) { 
    return (int64_t)((__int128_t)x*y >> 64); 
} 

A -O1 e superiori, questa compila a ciò che si vuole:

_mulHi: 
0000000000000000 movq %rsi,%rax 
0000000000000003 imulq %rdi 
0000000000000006 movq %rdx,%rax 
0000000000000009 ret 

credo che clang e VC++ anche avere il supporto per il tipo __int128_t, quindi questo dovrebbe funzionare anche su quelle piattaforme, con i soliti avvertimenti su come provarlo tu stesso.

Answer 4

Per quanto riguarda la soluzione di montaggio, non codificare le istruzioni mov! Lascia che sia il compilatore a farlo per te. Ecco una versione modificata del codice:

static long mull_hi(long inp1, long inp2) { 
    long output; 
    __asm__("imulq %2" 
      : "=d" (output) 
      : "a" (inp1), "r" (inp2)); 
    return output; 
} 

di riferimento di aiuto: Machine Constraints

Answer 5

Dal momento che hai fatto un buon lavoro di risolvere il tuo problema con il codice macchina, ho pensato che si meritava un aiuto con la versione portatile.Vorrei lasciare un ifdef in cui devi solo usare l'assembly se in gnu su x86.

In ogni caso, ecco un'implementazione ... Sono abbastanza sicuro che sia corretto, ma non ci sono garanzie, l'ho appena battuto ieri sera ... probabilmente dovresti liberarti delle statistiche positive_result [] e result_negative, quelli sono solo artefatti del mio test unitario ...

#include <stdlib.h> 
#include <stdio.h> 

// stdarg.h doesn't help much here because we need to call llabs() 

typedef unsigned long long uint64_t; 
typedef signed long long int64_t; 

#define B32 0xffffffffUL 

static uint64_t positive_result[2]; // used for testing 
static int result_negative;   // used for testing 

static void mixed(uint64_t *result, uint64_t innerTerm) 
{ 
    // the high part of innerTerm is actually the easy part 

    result[1] += innerTerm >> 32; 

    // the low order a*d might carry out of the low order result 

    uint64_t was = result[0]; 

    result[0] += (innerTerm & B32) << 32; 

    if (result[0] < was) // carry! 
     ++result[1]; 
} 


static uint64_t negate(uint64_t *result) 
{ 
    uint64_t t = result[0] = ~result[0]; 
    result[1] = ~result[1]; 
    if (++result[0] < t) 
    ++result[1]; 
    return result[1]; 
} 

uint64_t higherMul(int64_t sx, int64_t sy) 
{ 
    uint64_t x, y, result[2] = { 0 }, a, b, c, d; 

    x = (uint64_t)llabs(sx); 
    y = (uint64_t)llabs(sy); 

    a = x >> 32; 
    b = x & B32; 
    c = y >> 32; 
    d = y & B32; 

    // the highest and lowest order terms are easy 

    result[1] = a * c; 
    result[0] = b * d; 

    // now have the mixed terms ad + bc to worry about 

    mixed(result, a * d); 
    mixed(result, b * c); 

    // now deal with the sign 

    positive_result[0] = result[0]; 
    positive_result[1] = result[1]; 
    result_negative = sx < 0^sy < 0; 
    return result_negative ? negate(result) : result[1]; 
}