Skip to content

Conta ogni posizione di bit separatamente su molte bitmaschere a 64 bit, con AVX ma non con AVX2

Indaghiamo in diversi forum per avere per te la risposta al tuo dubbio, in caso di domande puoi lasciarci la tua preoccupazione e ti risponderemo con piacere, perché siamo qui per servirti.

Soluzione:

Sul mio sistema, un MacBook di 4 anni fa (intel core i5 da 2,7 GHz) con clang-900.0.39.2 -O3il vostro codice viene eseguito in 500 ms.

Basta cambiare il test interno in if ((pLong[j] & m) != 0) si risparmia il 30%, con un'esecuzione in 350 ms.

Semplificando ulteriormente la parte interna in target[i] += (pLong[j] >> i) & 1; senza un test, si riduce a 280 ms.

Ulteriori miglioramenti sembrano richiedere tecniche più avanzate, come la scomposizione dei bit in blocchi di 8 ulong e l'aggiunta di questi in parallelo, gestendo 255 ulong alla volta.

Ecco una versione migliorata che utilizza questo metodo e che viene eseguita in 45 ms sul mio sistema.

#include 
#include 
#include 
#include 
#include 
#include 

double getTS() {
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + tv.tv_usec / 1000000.0;
}

int main(int argc, char *argv[]) {
    unsigned int target[64] = { 0 };
    unsigned long *pLong = malloc(sizeof(*pLong) * 10000000);
    int i, j;

    if (!pLong) {
        printf("failed to allocaten");
        exit(1);
    }
    memset(pLong, 0xff, sizeof(*pLong) * 10000000);
    printf("p=%pn", (void*)pLong);
    double start = getTS();
    uint64_t inflate[256];
    for (i = 0; i < 256; i++) {
        uint64_t x = i;
        x = (x | (x << 28));
        x = (x | (x << 14));
        inflate[i] = (x | (x <<  7)) & 0x0101010101010101ULL;
    }
    for (j = 0; j < 10000000 / 255 * 255; j += 255) {
        uint64_t b[8] = { 0 };
        for (int k = 0; k < 255; k++) {
            uint64_t u = pLong[j + k];
            for (int kk = 0; kk < 8; kk++, u >>= 8)
                b[kk] += inflate[u & 255];
        }
        for (i = 0; i < 64; i++)
            target[i] += (b[i / 8] >> ((i % 8) * 8)) & 255;
    }
    for (; j < 10000000; j++) {
        uint64_t m = 1;
        for (i = 0; i < 64; i++) {
            target[i] += (pLong[j] >> i) & 1;
            m <<= 1;
        }
    }
    printf("target = {");
    for (i = 0; i < 64; i++)
        printf(" %d", target[i]);
    printf(" }n");
    printf("took %f secsn", getTS() - start);
    return 0;
}

La tecnica per gonfiare un byte in un long a 64 bit è studiata e spiegata nella risposta: https://stackoverflow.com/a/55059914/4593267 . Ho fatto il target come variabile locale, così come l'array inflate e stampo i risultati per garantire che il compilatore non ottimizzi i calcoli. In una versione di produzione si calcolerebbe l'array inflate separatamente.

L'uso diretto di SIMD potrebbe fornire ulteriori miglioramenti a scapito della portabilità e della leggibilità. Questo tipo di ottimizzazione è spesso meglio lasciarla al compilatore, che può generare codice specifico per l'architettura di destinazione. A meno che le prestazioni non siano critiche e il benchmarking dimostri che questo è un collo di bottiglia, preferirei sempre una soluzione generica.

Una soluzione diversa di njuffa fornisce prestazioni simili senza la necessità di un array precalcolato. A seconda del compilatore e delle specifiche hardware, potrebbe essere più veloce.

Correlato:

  • un precedente duplicato presenta alcune idee alternative: Come contare rapidamente i bit in contenitori separati in una serie di ints su Sandy Bridge?
  • Risposta di Harold sull'algoritmo di conteggio della popolazione di colonne AVX2 su ogni colonna di bit separatamente.
  • Matrix transpose and population count ha un paio di risposte utili con AVX2, inclusi i benchmark. Utilizza pezzi a 32 bit invece che a 64 bit.

Inoltre: https://github.com/mklarqvist/positional-popcount ha una miscela SSE, vari AVX2, vari AVX512, tra cui Harley-Seal che è ottimo per gli array di grandi dimensioni, e vari altri algoritmi per il popcount posizionale. Possibilmente solo per uint16_tma la maggior parte potrebbe essere adattata ad altre larghezze di parola. Credo che l'algoritmo che propongo di seguito sia quello che chiamano adder_forest.


La scelta migliore è SIMD, utilizzando AVX1 sulla CPU Sandybridge. I compilatori non sono abbastanza intelligenti da auto-vettorizzare i bit di loop-over per voi, anche se lo scrivete senza rami per dargli maggiori possibilità.

E sfortunatamente non sono abbastanza intelligenti da auto-vettorizzare la versione veloce che gradualmente si allarga e si aggiunge.


Vedere Esiste un'istruzione inversa all'istruzione movemask in intel avx2? per un riepilogo dei metodi di disimballaggio bitmap -> vettoriale per diverse dimensioni. Il suggerimento di Ext3h in un'altra risposta è buono: Disimballare i bit in qualcosa di più stretto dell'array di conteggio finale consente di ottenere più elementi per istruzione. I byte sono efficienti con SIMD, e poi si può fare fino a 255 verticale paddb senza overflow, prima di spacchettare per accumulare nell'array di contatori a 32 bit.

Ci vogliono solo 4 volte 16 byte __m128i per contenere tutti i 64 uint8_t quindi gli accumulatori possono rimanere nei registri, aggiungendosi alla memoria solo quando si allarga ai contatori a 32 bit in un ciclo esterno.

Lo spacchettamento non deve essere in ordine sparso.: si può sempre rimescolare target[] una volta alla fine, dopo aver accumulato tutti i risultati.

Il ciclo interno potrebbe essere srotolato per iniziare con un carico vettoriale a 64 o 128 bit, e scompattare in 4 o 8 modi diversi usando pshufb (_mm_shuffle_epi8).


Una strategia ancora migliore è quella di allargare gradualmente

Iniziando con accumulatori a 2 bit, poi mascherare/spostarli per portarli a 4 bit. In questo modo, nel ciclo più interno la maggior parte delle operazioni lavora con dati "densi", senza "diluirli" troppo subito. Una maggiore densità di informazione/entropia significa che ogni istruzione svolge un lavoro più utile.

L'uso di tecniche SWAR per l'addizione 32x a 2 bit all'interno di registri scalari o SIMD è facile/economico perché dobbiamo comunque evitare la possibilità di eseguire la parte superiore di un elemento. Con una SIMD corretta, perderemmo quei conteggi, con SWAR corromperemmo l'elemento successivo.

uint64_t x = *(input++);        // load a new bitmask
const uint64_t even_1bits = 0x5555555555555555;  // 0b...01010101;

uint64_t lo = x & even_1bits;
uint64_t hi = (x>>1) & even_1bits;            // or use ANDN before shifting to avoid a MOV copy

accum2_lo += lo;   // can do up to 3 iterations of this without overflow
accum2_hi += hi;   // because a 2-bit integer overflows at 4

Quindi si ripetono fino a 4 vettori di elementi a 4 bit, poi 8 vettori di elementi a 8 bit, poi si dovrebbe allargare fino a 32 e accumulare nell'array in memoria perché si esauriranno comunque i registri, e questo lavoro di loop esterno è abbastanza infrequente da non doverci preoccupare di passare a 16 bit. (Soprattutto se vettorializziamo manualmente).

Il più grande svantaggio: questo non autovettorizza, a differenza della versione di @njuffa. Ma con gcc -O3 -march=sandybridge per AVX1 (quindi eseguendo il codice su Skylake), l'esecuzione di questo scalare a 64 bit è in realtà ancora leggermente più veloce di 128 bit AVX auto-vettorizzato dal codice di @njuffa.

Ma questo è il timing su Skylake, che ha 4 porte ALU scalari (e la mov-elimination), mentre Sandybridge non ha la mov-elimination e ha solo 3 porte ALU, quindi il codice scalare probabilmente si scontrerà con i colli di bottiglia delle porte di esecuzione back-end. (Ma il codice SIMD potrebbe essere quasi altrettanto veloce, perché c'è un sacco di AND / ADD mescolati con gli shift, e SnB ha unità di esecuzione SIMD su tutte e 3 le sue porte che hanno ALU. Haswell ha appena aggiunto la porta 6, per il solo scalare, compresi shift e branch).

Con una buona vettorizzazione manuale, questo dovrebbe essere un fattore di quasi 2 o 4 volte più veloce.

Ma se si deve scegliere tra questo scalare e quello di @njuffa con l'autovettorizzazione AVX2, quello di @njuffa è più veloce su Skylake con -march=native

Se è possibile/richiesto costruire su un target a 32 bit, questo soffre molto (senza vettorizzazione a causa dell'uso di uint64_t nei registri a 32 bit), mentre il codice vettorizzato soffre appena (perché tutto il lavoro avviene in registri vettoriali della stessa larghezza).

// TODO: put the target[] re-ordering somewhere
// TODO: cleanup for N not a multiple of 3*4*21 = 252
// TODO: manual vectorize with __m128i, __m256i, and/or __m512i

void sum_gradual_widen (const uint64_t *restrict input, unsigned int *restrict target, size_t length)
{
    const uint64_t *endp = input + length - 3*4*21;     // 252 masks per outer iteration
    while(input <= endp) {
        uint64_t accum8[8] = {0};     // 8-bit accumulators
        for (int k=0 ; k<21 ; k++) {
            uint64_t accum4[4] = {0};  // 4-bit accumulators can hold counts up to 15.  We use 4*3=12
            for(int j=0 ; j<4 ; j++){
                uint64_t accum2_lo=0, accum2_hi=0;
                for(int i=0 ; i<3 ; i++) {  // the compiler should fully unroll this
                    uint64_t x = *input++;    // load a new bitmask
                    const uint64_t even_1bits = 0x5555555555555555;
                    uint64_t lo = x & even_1bits; // 0b...01010101;
                    uint64_t hi = (x>>1) & even_1bits;  // or use ANDN before shifting to avoid a MOV copy
                    accum2_lo += lo;
                    accum2_hi += hi;   // can do up to 3 iterations of this without overflow
                }

                const uint64_t even_2bits = 0x3333333333333333;
                accum4[0] +=  accum2_lo       & even_2bits;  // 0b...001100110011;   // same constant 4 times, because we shift *first*
                accum4[1] += (accum2_lo >> 2) & even_2bits;
                accum4[2] +=  accum2_hi       & even_2bits;
                accum4[3] += (accum2_hi >> 2) & even_2bits;
            }
            for (int i = 0 ; i<4 ; i++) {
                accum8[i*2 + 0] +=   accum4[i] & 0x0f0f0f0f0f0f0f0f;
                accum8[i*2 + 1] +=  (accum4[i] >> 4) & 0x0f0f0f0f0f0f0f0f;
            }
        }

        // char* can safely alias anything.
        unsigned char *narrow = (uint8_t*) accum8;
        for (int i=0 ; i<64 ; i++){
            target[i] += narrow[i];
        }
    }
    /* target[0] = bit 0
     * target[1] = bit 8
     * ...
     * target[8] = bit 1
     * target[9] = bit 9
     * ...
     */
    // TODO: 8x8 transpose
}

Non ci interessa l'ordine, quindi accum4[0] ha accumulatori a 4 bit per ogni quarto bit, per esempio. La correzione finale necessaria (ma non ancora implementata) alla fine è una trasposizione 8x8 dell'elemento uint32_t target[64] , che può essere fatto in modo efficiente usando unpck e vshufps con solo AVX1. (Trasporre un float 8x8 usando AVX/AVX2). E anche un ciclo di pulizia per le ultime maschere fino a 251.

Possiamo usare qualsiasi larghezza di elemento SIMD per implementare questi spostamenti; dobbiamo comunque mascherare per larghezze inferiori a 16 bit (SSE/AVX non ha spostamenti a granularità di byte, ma solo a 16 bit minimo).

Risultati dei benchmark su Arch Linux i7-6700k dal test harness di @njuffa, con l'aggiunta di questo. (Godbolt) N = (10000000 / (3*4*21) * 3*4*21) = 9999864 (cioè 10000000 arrotondato per difetto a un multiplo del fattore di "srotolamento" di 252 iterazioni, quindi la mia implementazione semplicistica sta facendo la stessa quantità di lavoro, senza contare il riordino target[] che non viene fatto, quindi stampa risultati non corrispondenti.
Ma i conteggi stampati corrispondono a un'altra posizione dell'array di riferimento).

Ho eseguito il programma per 4 volte di seguito (per assicurarmi che la CPU fosse riscaldata al massimo del turbo) e ho preso una delle esecuzioni che sembrava buona (nessuna delle 3 volte era anormalmente alta).

rif: il miglior bit-loop (sezione successiva)
veloce: Il codice di @njuffa. (autovettorizzato con istruzioni intere AVX a 128 bit).
graduale: la mia versione (non auto-vettorizzata da gcc o clang, almeno non nel ciclo interno). gcc e clang srotolano completamente le 12 iterazioni interne.

  • gcc8.2 -O3 -march=sandybridge -fpie -no-pie
    ref: 0,331373 secondi, fast: 0,011387 secondi, gradual: 0,009966 secondi
  • gcc8.2 -O3 -march=sandybridge -fno-pie -no-pie
    rif: 0,397175 sec, veloce: 0,011255 sec, graduale: 0,010018 sec
  • clang7.0 -O3 -march=sandybridge -fpie -no-pie
    ref: 0,352381 sec, veloce: 0,011926 sec, graduale: 0,009269 sec (conteggi molto bassi per gli uops della porta 7, clang ha usato l'indirizzamento indicizzato per gli store)
  • clang7.0 -O3 -march=sandybridge -fno-pie -no-pie
    rif: 0,293014 sec veloce: 0,011777 sec, graduale: 0,009235 sec.

-marcia=lago di cielo (che consente AVX2 per vettori di interi a 256 bit) aiuta entrambi, ma quello di @njuffa di più perché ne vettorizza di più (compreso il suo ciclo più interno):

  • gcc8.2 -O3 -march=skylake -fpie -no-pie
    ref: 0,328725 sec, fast: 0,007621 sec, gradual: 0,010054 sec (gcc non mostra alcun guadagno per "gradual", solo "fast")
  • gcc8.2 -O3 -march=skylake -fno-pie -no-pie
    rif: 0,333922 sec, veloce: 0,007620 sec, graduale: 0,009866 sec

  • clang7.0 -O3 -march=skylake -fpie -no-pie
    rif: 0,260616 sec, veloce: 0,007521 sec, graduale: 0,008535 sec (non capisco perché graduale sia più veloce di -march=sandybridg)e; non usa BMI1 andn. Credo perché utilizza AVX2 a 256 bit per il ciclo esterno k=0..20 con vpaddq)

  • clang7.0 -O3 -march=skylake -fno-pie -no-pie
    rif. 0,259159 sec., veloce: 0,007496 secondi graduale: 0,008671 sec.

Senza AVX, solo SSE4.2: (-march=nehalem), stranamente il graduale di clang è più veloce che con AVX / tune=sandybridge. "veloce" è solo appena più lento che con AVX.

  • gcc8.2 -O3 -march=skylake -fno-pie -no-pie
    ref: 0,337178 secondi, fast: 0,011983 secondi, gradual: 0,010587 secondi
  • clang7.0 -O3 -march=skylake -fno-pie -no-pie
    rif: 0,293555 sec veloce: 0,012549 sec, graduale: 0,008697 sec.

-fprofile-generate / -fprofile-use aiuta un po' per GCC, specialmente per la versione "ref" che non esegue alcuno scorrimento per impostazione predefinita.

Ho evidenziato i migliori, ma spesso si trovano all'interno di un margine di rumore di misura l'uno dall'altro. Non sorprende che il -fno-pie -no-pie sia a volte più veloce: l'indicizzazione di array statici con [disp32 + reg] è non una modalità di indirizzamento indicizzata, solo base + disp32, quindi non si sbloccherà mai sulle CPU Sandybridge.

Ma con gcc a volte -fpie era più veloce; non ho controllato, ma presumo che gcc si sia dato la zappa sui piedi in qualche modo quando l'indirizzamento assoluto a 32 bit era possibile. Oppure è capitato che differenze apparentemente innocenti nella generazione del codice causassero problemi di allineamento o di uop-cache; non ho controllato in dettaglio.


Per quanto riguarda le SIMD, possiamo semplicemente fare 2 o 4x uint64_t in parallelo, accumulando orizzontalmente solo nel passaggio finale in cui allarghiamo i byte a elementi a 32 bit. (Forse rimescolando in corsia e poi usando pmaddubsw con un moltiplicatore di _mm256_set1_epi8(1) per aggiungere coppie di byte orizzontali in elementi a 16 bit).

TODO: vettorizzato manualmente __m128i e __m256i (e __m512i) di questa versione. Dovrebbe essere quasi 2x, 4x o addirittura 8x più veloce dei tempi "graduali" di cui sopra. Probabilmente il prefetch HW può ancora tenere il passo, tranne forse una versione AVX512 con dati provenienti dalla DRAM, specialmente se c'è contesa da altri thread. Facciamo una quantità significativa di lavoro per ogni qword letta.


Codice obsoleto: miglioramenti al bit-loop

Anche la versione scalare portatile può essere migliorata, velocizzandola da ~1,92 secondi (con un tasso di mancata previsione del ramo del 34%. con i loop veloci commentati!) a ~0,35 secondi (clang7.0 -O3 -march=sandybridge) con un input adeguatamente casuale su Skylake a 3,9GHz. O 1,83 secondi per la versione con rami con != 0 invece di == mperché i compilatori non riescono a dimostrare che m ha sempre esattamente 1 bit impostato e/o ottimizzano di conseguenza.

(contro 0,01 sec per la versione di @njuffa o la mia versione veloce, quindi è abbastanza inutile in senso assoluto, ma vale la pena di menzionarlo come esempio di ottimizzazione generale di quando usare il codice senza rami).

Se ci si aspetta un mix casuale di zeri e uno, si vuole qualcosa di non ramificato che non faccia previsioni errate. Fare += 0 per gli elementi che erano zero evita questo problema e significa anche che la macchina astratta C tocca sicuramente quella memoria, indipendentemente dai dati.

I compilatori non sono autorizzati a inventare le scritture, quindi se volessero auto-vettorizzare il vostro if() target[i]++ dovrebbero usare un archivio mascherato come quello di x86 vmaskmovps per evitare una lettura/riscrittura non atomica degli elementi non modificati di target. Quindi un ipotetico compilatore futuro in grado di autovettorizzare il codice scalare semplice avrebbe un tempo più facile con questo.

Comunque, un modo per scrivere questo è target[i] += (pLong[j] & m != 0);, utilizzando la conversione bool->int per ottenere un intero 0 / 1.

Ma otteniamo un asm migliore per x86 (e probabilmente per la maggior parte delle altre architetture) se spostiamo i dati e isoliamo il bit basso con &1. I compilatori sono piuttosto stupidi e non sembrano accorgersi di questa ottimizzazione. Ottimizzano invece bene il contatore di loop in più e trasformano in m <<= 1 in add same,same per spostare in modo efficiente a sinistra, ma utilizzano ancora xor-zero / test / setne per creare un intero 0 / 1.

Un ciclo interno come questo viene compilato in modo leggermente più efficiente (ma ancora molto peggiore di quello che possiamo fare con SSE2 o AVX, o anche scalare usando la tabella di lookup di @chrqlie che rimarrà calda in L1d se usata ripetutamente come questa, permettendo a SWAR in uint64_t):

    for (int j = 0; j < 10000000; j++) {
#if 1  // extract low bit directly
        unsigned long long tmp = pLong[j];
        for (int i=0 ; i<64 ; i++) {   // while(tmp) could mispredict, but good for sparse data
            target[i] += tmp&1;
            tmp >>= 1;
        }
#else // bool -> int shifting a mask
        unsigned long m = 1;
        for (i = 0; i < 64; i++) {
            target[i]+= (pLong[j] & m) != 0;
            m = (m << 1);
        }
#endif

Si noti che unsigned long non è garantito che sia un tipo a 64 bit e non lo è in x86-64 System V x32 (ILP32 in modalità a 64 bit) e Windows x64. O in ABI a 32 bit come i386 System V.

Compilato sul compilatore Godbolt explorer da gcc, clang e ICC, con gcc c'è 1 uops in meno nel ciclo. Ma sono tutti scalari, con clang e ICC che li srotolano di 2 unità.

# clang7.0 -O3 -march=sandybridge
.LBB1_2:                            # =>This Loop Header: Depth=1
   # outer loop loads a uint64 from the src
    mov     rdx, qword ptr [r14 + 8*rbx]
    mov     rsi, -256
.LBB1_3:                            #   Parent Loop BB1_2 Depth=1
                                    # do {
    mov     edi, edx
    and     edi, 1                              # isolate the low bit
    add     dword ptr [rsi + target+256], edi   # and += into target

    mov     edi, edx
    shr     edi
    and     edi, 1                              # isolate the 2nd bit
    add     dword ptr [rsi + target+260], edi

    shr     rdx, 2                              # tmp >>= 2;

    add     rsi, 8
    jne     .LBB1_3                       # } while(offset += 8 != 0);

È un risultato leggermente migliore di quello ottenuto con test / setnz. Senza lo scorrimento, bt / setc avrebbe potuto essere uguale, ma i compilatori non sono bravi ad usare bt per implementare bool (x & (1ULL << n))oppure bts per implementare x |= 1ULL << n.

Se molte parole hanno il loro bit più alto impostato molto al di sotto del bit 63, il ciclo su while(tmp) potrebbe essere una vittoria. I mispredicts del ramo non valgono la pena se si risparmiano solo ~0-4 iterazioni la maggior parte delle volte, ma se si risparmiano spesso 32 iterazioni, potrebbe valerne davvero la pena. Forse si potrebbe srotolare nel sorgente in modo che il ciclo verifichi solo tmp ogni 2 iterazioni (perché i compilatori non faranno questa trasformazione per voi), ma poi il ramo del ciclo può essere shr rdx, 2 / jnz.

Su Sandybridge-family, si tratta di 11 uops a dominio fuso per il front-end per 2 bit di input. (add [mem], reg con una modalità di indirizzamento non indicizzata microfusibilizza il load+ALU, e lo store-address+store-data, tutto il resto è single-uop. add/jcc macrofusibilizza. Vedere la guida di Agner Fog e https://stackoverflow.com/tags/x86/info). Quindi dovrebbe funzionare a qualcosa come 3 cicli per 2 bit = un uint64_t per 96 cicli. (Sandybridge non esegue lo "srotolamento" interno del suo loop buffer, quindi i conteggi degli uop non multipli di 4 si arrotondano per eccesso, a differenza di Haswell e successivi).

contro la versione non srotolata di gcc che è di 7 uop per 1 bit = 2 cicli per bit. Se si compila con gcc -O3 -march=native -fprofile-generate / test-run / gcc -O3 -march=native -fprofile-usel'ottimizzazione guidata dal profilo consentirebbe lo srotolamento del ciclo.

Questo è probabilmente più lento di una versione ramificata su dati perfettamente prevedibili come quella ottenuta da memset con qualsiasi schema di byte ripetuto. Suggerirei di riempire il vostro array con dati generati casualmente da un PRNG veloce come un SSE2 xorshift+, o se state solo cronometrando il ciclo di conteggio allora usate qualsiasi cosa vogliate, come ad esempio rand().

Un modo per velocizzare significativamente questa operazione, anche senza AVX, è quello di suddividere i dati in blocchi di massimo 255 elementi e accumulare i conteggi dei bit in senso byte in un ordinario uint64_t in variabili ordinarie. Poiché i dati di partenza hanno 64 bit, abbiamo bisogno di un array di 8 accumulatori byte-wise. Il primo accumulatore conta i bit nelle posizioni 0, 8, 16, ... 56, il secondo accumulatore conta i bit nelle posizioni 1, 9, 17, ... 57; e così via. Una volta terminata l'elaborazione di un blocco di dati, si trasferiscono i conteggi dell'accumulatore byte-wise nell'accumulatore target conteggi. Una funzione per aggiornare l'accumulatore target per un blocco di massimo 255 numeri può essere codificata in modo semplice secondo la descrizione precedente, dove BITS è il numero di bit nei dati di origine:

/* update the counts of 1-bits in each bit position for up to 255 numbers */
void sum_block (const uint64_t *pLong, unsigned int *target, int lo, int hi)
{
    int jj, k, kk;
    uint64_t byte_wise_sum [BITS/8] = {0};
    for (jj = lo; jj < hi; jj++) {
        uint64_t t = pLong[jj];
        for (k = 0; k < BITS/8; k++) {
            byte_wise_sum[k] += t & 0x0101010101010101;
            t >>= 1;
        }
    }
    /* accumulate byte sums into target */
    for (k = 0; k < BITS/8; k++) {
        for (kk = 0; kk < BITS; kk += 8) {
            target[kk + k] += (byte_wise_sum[k] >> kk) & 0xff;
        }
    }
}

L'intero programma ISO-C99, che dovrebbe essere in grado di funzionare almeno su piattaforme Windows e Linux, è mostrato di seguito. Inizializza i dati di origine con un PRNG, esegue un controllo di correttezza rispetto all'implementazione di riferimento del richiedente e confronta sia il codice di riferimento che la versione accelerata. Sulla mia macchina (Intel Xeon E3-1270 v2 @ 3,50 GHz), quando viene compilato con MSVS 2010 alla massima ottimizzazione (/Ox), l'output del programma è:

p=0000000000550040
ref took 2.020282 secs, fast took 0.027099 secs

dove ref si riferisce alla soluzione originale del richiedente. L'aumento di velocità in questo caso è di circa un fattore 74x. Con altri compilatori (e soprattutto con quelli più recenti) si osserveranno accelerazioni diverse.

#include 
#include 
#include 
#include 

#if defined(_WIN32)
#if !defined(WIN32_LEAN_AND_MEAN)
#define WIN32_LEAN_AND_MEAN
#endif
#include 
double second (void)
{
    LARGE_INTEGER t;
    static double oofreq;
    static int checkedForHighResTimer;
    static BOOL hasHighResTimer;

    if (!checkedForHighResTimer) {
        hasHighResTimer = QueryPerformanceFrequency (&t);
        oofreq = 1.0 / (double)t.QuadPart;
        checkedForHighResTimer = 1;
    }
    if (hasHighResTimer) {
        QueryPerformanceCounter (&t);
        return (double)t.QuadPart * oofreq;
    } else {
        return (double)GetTickCount() * 1.0e-3;
    }
}
#elif defined(__linux__) || defined(__APPLE__)
#include 
#include 
double second (void)
{
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return (double)tv.tv_sec + (double)tv.tv_usec * 1.0e-6;
}
#else
#error unsupported platform
#endif

/*
  From: geo <[email protected]>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, 
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, 
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), 
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

#define N          (10000000)
#define BITS       (64)
#define BLOCK_SIZE (255)

/* cupdate the count of 1-bits in each bit position for up to 255 numbers */
void sum_block (const uint64_t *pLong, unsigned int *target, int lo, int hi)
{
    int jj, k, kk;
    uint64_t byte_wise_sum [BITS/8] = {0};
    for (jj = lo; jj < hi; jj++) {
        uint64_t t = pLong[jj];
        for (k = 0; k < BITS/8; k++) {
            byte_wise_sum[k] += t & 0x0101010101010101;
            t >>= 1;
        }
    }
    /* accumulate byte sums into target */
    for (k = 0; k < BITS/8; k++) {
        for (kk = 0; kk < BITS; kk += 8) {
            target[kk + k] += (byte_wise_sum[k] >> kk) & 0xff;
        }
    }
}

int main (void) 
{
    double start_ref, stop_ref, start, stop;
    uint64_t *pLong;
    unsigned int target_ref [BITS] = {0};
    unsigned int target [BITS] = {0};
    int i, j;

    pLong = malloc (sizeof(pLong[0]) * N);
    if (!pLong) {
        printf("failed to allocaten");
        return EXIT_FAILURE;
    }
    printf("p=%pn", pLong);

    /* init data */
    for (j = 0; j < N; j++) {
        pLong[j] = KISS64;
    }

    /* count bits slowly */
    start_ref = second();
    for (j = 0; j < N; j++) {
        uint64_t m = 1;
        for (i = 0; i < BITS; i++) {
            if ((pLong[j] & m) == m) {
                target_ref[i]++;
            }
            m = (m << 1);
        }
    }
    stop_ref = second();

    /* count bits fast */
    start = second();
    for (j = 0; j < N / BLOCK_SIZE; j++) {
        sum_block (pLong, target, j * BLOCK_SIZE, (j+1) * BLOCK_SIZE);
    }
    sum_block (pLong, target, j * BLOCK_SIZE, N);
    stop = second();

    /* check whether result is correct */
    for (i = 0; i < BITS; i++) {
        if (target[i] != target_ref[i]) {
            printf ("error @ %d: res=%u ref=%un", i, target[i], target_ref[i]);
        }
    }

    /* print benchmark results */
    printf("ref took %f secs, fast took %f secsn", stop_ref - start_ref, stop - start);
    return EXIT_SUCCESS;
}

Se hai qualche antipatia o un modo per correggere la nostra sezione, puoi suonare una nota e noi la interpreteremo volentieri.



Utilizzate il nostro motore di ricerca

Ricerca
Generic filters

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.