Skip to content

Compter chaque position de bit séparément sur plusieurs masques de bits de 64 bits, avec AVX mais pas AVX2.

Notre groupe de rédacteurs a longuement investigué la solution à vos recherches, nous vous donnons les réponses pour cela, notre souhait est qu'il vous soit d'un grand soutien.

Solution :

Sur mon système, un MacBook de 4 ans (2.7 GHz intel core i5) avec clang-900.0.39.2 -O3, votre code s'exécute en 500ms.

Il suffit de changer le test interne en if ((pLong[j] & m) != 0) permet d'économiser 30%, en s'exécutant en 350ms.

En simplifiant encore la partie interne à target[i] += (pLong[j] >> i) & 1; sans test la ramène à 280ms.

D'autres améliorations semblent nécessiter des techniques plus avancées comme le déballage des bits en blocs de 8 ulongs et l'ajout de ceux-ci en parallèle, en traitant 255 ulongs à la fois.

Voici une version améliorée utilisant cette méthode. elle s'exécute en 45ms sur mon système.

#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 technique pour gonfler un octet à une longue de 64 bits sont étudiés et expliqués dans la réponse : https://stackoverflow.com/a/55059914/4593267 . J'ai fait le target une variable locale, ainsi que le tableau inflate et j'imprime les résultats pour m'assurer que le compilateur n'optimisera pas les calculs. Dans une version de production, vous calculeriez le tableau inflate séparément.

L'utilisation directe de SIMD pourrait fournir des améliorations supplémentaires au détriment de la portabilité et de la lisibilité. Ce type d'optimisation est souvent mieux laissé au compilateur car il peut générer un code spécifique pour l'architecture cible. À moins que les performances soient critiques et que le benchmarking prouve que c'est un goulot d'étranglement, je privilégierais toujours une solution générique.

Une solution différente par njuffa fournit des performances similaires sans le besoin d'un tableau précalculé. Selon les spécificités de votre compilateur et de votre matériel, elle pourrait être plus rapide.

Relié :

  • un duplicata antérieur a quelques idées alternatives : Comment compter rapidement les bits dans des bacs séparés dans une série d'ints sur Sandy Bridge....
  • Réponse d'Harold sur l'algorithme de comptage de la population des colonnes AVX2 sur chaque colonne de bits séparément.
  • Transposition matricielle et comptage de population a quelques réponses utiles avec AVX2, y compris des benchmarks. Il utilise des morceaux de 32 bits au lieu de 64 bits.

Aussi : https://github.com/mklarqvist/positional-popcount a SSE blend, divers AVX2, divers AVX512 dont Harley-Seal qui est génial pour les grands tableaux, et divers autres algorithmes pour le popcount positionnel. Possiblement seulement pour uint16_t, mais la plupart pourraient être adaptés à d'autres largeurs de mots. Je pense que l'algorithme que je propose ci-dessous est ce qu'ils appellent... adder_forest.


Votre meilleure chance est SIMD, en utilisant AVX1 sur votre CPU Sandybridge. Les compilateurs ne sont pas assez intelligents pour auto-vectoriser votre boucle sur les bits pour vous, même si vous l'écrivez sans branche pour leur donner une meilleure chance.

Et malheureusement pas assez intelligents pour auto-vectoriser la version rapide qui s'élargit et s'ajoute progressivement.


Voir is there an inverse instruction to the movemask instruction in intel avx2 ? pour un résumé des méthodes bitmap -> vector unpack pour différentes tailles. La suggestion d'Ext3h dans une autre réponse est bonne : Le déballage des bits vers quelque chose de plus étroit que le tableau de comptage final vous donne plus d'éléments par instruction. Les octets sont efficaces avec SIMD, et ensuite vous pouvez faire jusqu'à 255 verticaux. paddb sans débordement, avant de déballer pour accumuler dans le tableau de comptage de 32 bits.

Il ne faut que 4x 16 octets __m128i pour contenir les 64 uint8_t éléments, donc ces accumulateurs peuvent rester dans les registres, ne s'ajoutant à la mémoire que lors de l'élargissement aux compteurs 32 bits dans une boucle externe.

Le déballage n'a pas besoin d'être dans l'ordre.: on peut toujours mélanger target[] une fois à la toute fin, après avoir accumulé tous les résultats.

La boucle interne pourrait être déroulée pour commencer avec un chargement de vecteur de 64 ou 128 bits, et déballer de 4 ou 8 façons différentes en utilisant. pshufb (_mm_shuffle_epi8).


Une stratégie encore meilleure est d'élargir progressivement

En commençant par des accumulateurs de 2 bits, puis en masquant/décalant pour élargir ceux-ci à 4 bits. Ainsi, dans la boucle la plus interne, la plupart des opérations travaillent avec des données "denses", ne les "diluant" pas trop tout de suite. Une densité d'information / entropie plus élevée signifie que chaque instruction fait un travail plus utile.

L'utilisation des techniques SWAR pour 32x 2-bit add à l'intérieur des registres scalaires ou SIMD est facile / bon marché parce que nous devons éviter la possibilité de porter le haut d'un élément de toute façon. Avec un SIMD approprié, nous perdrions ces comptes, avec SWAR nous corromprions l'élément suivant.

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

Ensuite, vous répétez jusqu'à 4 vecteurs d'éléments de 4 bits, puis 8 vecteurs d'éléments de 8 bits, puis vous devriez élargir jusqu'à 32 et accumuler dans le tableau en mémoire parce que vous serez à court de registres de toute façon, et ce travail de boucle externe est assez peu fréquent pour que nous n'ayons pas besoin de nous embêter à passer à 16 bits. (Surtout si nous vectorisons manuellement).

Le plus gros inconvénient : ceci ne auto-vectoriser, contrairement à la version de @njuffa. Mais avec gcc -O3 -march=sandybridge pour AVX1 (puis en exécutant le code sur Skylake), cette exécution scalaire 64 bits est en fait encore légèrement plus rapide que l'asm 128-bit AVX auto-vectorisé du code de @njuffa.

Mais c'est le timing sur Skylake, qui a 4 ports ALU scalaires (et l'élimination des mov), tandis que Sandybridge n'a pas d'élimination des mov et n'a que 3 ports ALU, donc le code scalaire se heurtera probablement à des goulots d'étranglement de port d'exécution back-end. (Mais le code SIMD peut être presque aussi rapide, parce qu'il y a beaucoup de AND / ADD mélangé avec les shifts, et SnB a des unités d'exécution SIMD sur ses 3 ports qui ont des ALUs sur eux. Haswell vient d'ajouter le port 6, pour les scalaires uniquement, y compris les décalages et les branches).

Avec une bonne vectorisation manuelle, cela devrait être un facteur de presque 2 ou 4 plus rapide.

Mais si vous devez choisir entre ce scalaire ou celui de @njuffa avec l'autovectorisation AVX2, celui de @njuffa est plus rapide sur Skylake avec... -march=native

Si la construction sur une cible 32 bits est possible/exigée, cela souffre beaucoup (sans vectorisation à cause de l'utilisation de uint64_t dans les registres 32 bits), alors que le code vectorisé souffre à peine (parce que tout le travail se passe dans des registres vectoriels de même largeur).

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

Nous ne nous soucions pas de l'ordre, donc accum4[0] a des accumulateurs de 4 bits pour chaque 4ème bit, par exemple. Le dernier correctif nécessaire (mais pas encore implémenté) à la toute fin est une transposition 8x8 de l'élément uint32_t target[64] tableau, ce qui peut être fait efficacement en utilisant unpck et vshufps avec seulement AVX1. (Transpose un float 8x8 en utilisant AVX/AVX2). Et aussi une boucle de nettoyage pour les derniers masques jusqu'à 251.

Nous pouvons utiliser n'importe quelle largeur d'élément SIMD pour implémenter ces décalages ; nous devons de toute façon masquer pour les largeurs inférieures à 16 bits (SSE/AVX n'a pas de décalages de granularité d'octet, seulement 16 bits minimum).

Résultats des benchmarks sur Arch Linux i7-6700k. à partir du harnais de test de @njuffa, avec cet ajout. (Godbolt) N = (10000000 / (3*4*21) * 3*4*21) = 9999864 (c'est-à-dire 10000000 arrondi à un multiple du facteur de "déroulement" de 252 itérations, donc mon implémentation simpliste fait la même quantité de travail, sans compter le ré-ordonnancement. target[] ce qu'elle ne fait pas, donc elle imprime des résultats non conformes.
Mais les comptes imprimés correspondent à une autre position du tableau de référence).

J'ai exécuté le programme 4x à la suite (pour m'assurer que le CPU était chauffé au turbo maximum) et j'ai pris une des exécutions qui semblait bonne (aucune des 3 fois anormalement élevée).

ref : la meilleure boucle de bits (section suivante).
rapide : Le code de @njuffa. (auto-vectorisé avec des instructions entières AVX de 128 bits).
graduel : ma version (pas auto-vectorisée par gcc ou clang, du moins pas dans la boucle interne.) gcc et clang déroulent complètement les 12 itérations internes.

  • gcc8.2 -O3 -march=sandybridge -fpie -no-pie
    ref : 0,331373 secs, rapide : 0,011387 secs, graduel : 0,009966 secs.
  • gcc8.2 -O3 -march=sandybridge -fno-pie -no-pie
    ref : 0.397175 secs, rapide : 0.011255 secs, graduel : 0.010018 secs
  • clang7.0 -O3 -march=sandybridge -fpie -no-pie
    ref : 0.352381 secs, rapide : 0.011926 secs, graduel : 0.009269 secs (très faible nombre d'uops pour le port 7, clang a utilisé l'adressage indexé pour les magasins)
  • clang7.0 -O3 -march=sandybridge -fno-pie -no-pie
    ref : 0.293014 secs, rapide : 0.011777 secs, graduel : 0.009235 sec

-march=skylake (permettant AVX2 pour les vecteurs entiers de 256 bits) aide les deux, mais celui de @njuffa le plus parce qu'une plus grande partie de celui-ci vectorise (y compris sa boucle la plus interne) :

  • gcc8.2 -O3 -march=skylake -fpie -no-pie
    ref : 0,328725 sec, fast : 0,007621 sec, gradual : 0,010054 sec (gcc ne montre aucun gain pour "gradual", seulement "fast").
  • gcc8.2 -O3 -march=skylake -fno-pie -no-pie
    ref : 0.333922 secs, rapide : 0.007620 secs, graduel : 0.009866 secs

  • clang7.0 -O3 -march=skylake -fpie -no-pie
    ref : 0.260616 secs, rapide : 0.007521 secs, graduel : 0.008535 secs (IDK pourquoi graduel est plus rapide que -march=sandybridg)e; il n'utilise pas BMI1 andn. Je suppose que c'est parce qu'il utilise 256-bit AVX2 pour la boucle externe k=0..20 avec... vpaddq)

  • clang7.0 -O3 -march=skylake -fno-pie -no-pie
    ref : 0.259159 secs, rapide : 0.007496 sec graduel : 0.008671 secs

Sans AVX, seulement SSE4.2 : (-march=nehalem), bizarrement le graduel de clang est plus rapide qu'avec AVX / tune=sandybridge. "rapide" est à peine plus lent qu'avec AVX.

  • gcc8.2 -O3 -march=skylake -fno-pie -no-pie
    ref : 0,337178 secs, rapide : 0,011983 secs, graduel : 0,010587 secs.
  • clang7.0 -O3 -march=skylake -fno-pie -no-pie
    ref : 0.293555 secs rapide : 0.012549 sec, graduel : 0.008697 sec

-fprofile-generate / -fprofile-use aide un peu pour GCC, surtout pour la version "ref" où il ne déroule pas du tout par défaut.

J'ai mis en évidence les meilleurs, mais souvent ils sont dans la marge de bruit de mesure les uns des autres. Il n'est pas surprenant que les -fno-pie -no-pie était parfois plus rapide : indexer des tableaux statiques avec [disp32 + reg] est pas un mode d'adressage indexé, juste base + disp32, donc il ne se décolle jamais sur les processeurs de la famille Sandybridge.

Mais avec gcc parfois -fpie était plus rapide ; je n'ai pas vérifié mais je suppose que gcc s'est juste tiré une balle dans le pied d'une manière ou d'une autre quand l'adressage absolu 32 bits était possible. Ou juste des différences d'apparence innocente dans le code-gen sont arrivées à causer des problèmes d'alignement ou de uop-cache ; je n'ai pas vérifié en détail.


Pour le SIMD, on peut simplement faire 2 ou 4x... uint64_t en parallèle, en accumulant seulement horizontalement dans l'étape finale où nous élargissons les octets à des éléments de 32 bits. (Peut-être en mélangeant dans le couloir et ensuite en utilisant pmaddubsw avec un multiplicateur de _mm256_set1_epi8(1) pour ajouter les paires d'octets horizontaux en éléments de 16 bits).

TODO : vectorisation manuelle __m128i et __m256i (et __m512i) de cette version. Devrait être près de 2x, 4x, ou même 8x plus rapide que les temps "graduels" ci-dessus. Probablement que HW prefetch peut encore suivre, sauf peut-être une version AVX512 avec des données provenant de la DRAM, surtout s'il y a de la contention d'autres threads. Nous faisons une quantité significative de travail par qword que nous lisons.


Code obsolète : améliorations de la boucle de bits.

Votre version scalaire portable peut aussi être améliorée, en l'accélérant de ~1.92 secondes (avec un taux de mauvaise prédiction de branche de 34% dans l'ensemble. avec les boucles rapides commentées !) à ~0.35sec (clang7.0 -O3 -march=sandybridge) avec une entrée correctement aléatoire sur 3.9GHz Skylake. Ou 1,83 sec pour la version branchée avec != 0 au lieu de == mparce que les compilateurs ne parviennent pas à prouver que m a toujours exactement 1 bit activé et/ou l'optimisent en conséquence.

(contre 0,01 sec pour la version rapide de @njuffa ou la mienne ci-dessus, donc c'est plutôt inutile dans un sens absolu, mais cela vaut la peine d'être mentionné comme un exemple d'optimisation générale pour savoir quand utiliser du code sans branche).

Si vous vous attendez à un mélange aléatoire de zéros et de uns, vous voulez quelque chose sans branche qui ne se trompe pas de prédiction. Faire += 0 pour les éléments qui étaient zéro évite cela, et signifie également que la machine abstraite C touche définitivement cette mémoire indépendamment des données.

Les compilateurs ne sont pas autorisés à inventer des écritures, donc s'ils voulaient auto-vectoriser votre élément if() target[i]++ version, ils auraient à utiliser un magasin masqué comme x86. vmaskmovps pour éviter une lecture / réécriture non-atomique d'éléments non modifiés de target. Ainsi, un futur compilateur hypothétique qui peut auto-vectoriser le code scalaire simple aurait un temps plus facile avec cela.

Quoi qu'il en soit, une façon d'écrire ceci est target[i] += (pLong[j] & m != 0);, en utilisant la conversion bool->int pour obtenir un entier 0 / 1.

Mais nous obtenons un meilleur asm pour x86 (et probablement pour la plupart des autres architectures) si nous décalons simplement les données et isolons le bit de poids faible avec... &1. Les compilateurs sont un peu bêtes et ne semblent pas repérer cette optimisation. Ils optimisent bien l'élimination du compteur de boucle supplémentaire et transforment m <<= 1 en add same,same pour décaler efficacement à gauche, mais ils utilisent toujours xor-zero / test / setne pour créer un entier 0 / 1.

Une boucle interne comme celle-ci compile légèrement plus efficacement (mais toujours. beaucoup beaucoup pire que ce que nous pouvons faire avec SSE2 ou AVX, ou même scalaire en utilisant la table de consultation de @chrqlie qui restera chaude en L1d lorsqu'elle est utilisée de manière répétée comme ceci, permettant à SWAR en 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

Notez que unsigned long n'est pas garanti d'être un type 64-bit, et ne l'est pas dans x86-64 System V x32 (ILP32 en mode 64-bit), et Windows x64. Ou dans les ABI 32 bits comme i386 System V.

Compilé sur le Godbolt compiler explorer par gcc, clang, et ICC, c'est 1 uops de moins dans la boucle avec gcc. Mais tous sont simplement scalaires, avec clang et ICC déroulant par 2.

# 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);

C'est légèrement mieux que ce que nous obtenons de test / setnz. Sans déroulage, bt / setc aurait pu être égal, mais les compilateurs sont mauvais dans l'utilisation de l'expression bt pour implémenter bool (x & (1ULL << n))ou bts pour mettre en oeuvre x |= 1ULL << n.

Si de nombreux mots ont leur bit activé le plus élevé bien en dessous du bit 63, le fait de boucler sur while(tmp) pourrait être une victoire.. Les erreurs de prédiction de branche font que cela ne vaut pas la peine si cela n'économise que ~0 à 4 itérations la plupart du temps, mais si cela économise souvent 32 itérations, cela pourrait vraiment valoir la peine. Peut-être dérouler dans la source pour que la boucle teste seulement tmp toutes les 2 itérations (parce que les compilateurs ne feront pas cette transformation pour vous), mais alors la branche de la boucle peut être shr rdx, 2 / jnz.

Sur la famille Sandybridge, cela représente 11 uops de domaine fusionné pour le front-end par 2 bits d'entrée. (add [mem], reg avec un mode d'adressage non-indexé micro-fuse le load+ALU, et le store-address+store-data, tout le reste est single-uop. add/jcc macro-fuses. Voir le guide d'Agner Fog, et https://stackoverflow.com/tags/x86/info). Donc il devrait fonctionner à quelque chose comme 3 cycles par 2 bits = un uint64_t par 96 cycles. (Sandybridge ne "déroule" pas en interne dans son tampon de boucle, donc les comptages d'uop non-multiples de 4 s'arrondissent fondamentalement vers le haut, contrairement à Haswell et plus tard).

contre la version non déroulée de gcc qui est de 7 uop par 1 bit = 2 cycles par bit. Si vous avez compilé avec gcc -O3 -march=native -fprofile-generate / test-run / gcc -O3 -march=native -fprofile-use, l'optimisation guidée par le profil permettrait de dérouler les boucles.

C'est probablement plus lent qu'une version branchée sur des données parfaitement prévisibles comme on obtient à partir de... memset avec n'importe quel modèle d'octet répétitif.. Je suggérerais de remplir votre tableau avec des données générées aléatoirement à partir d'un PRNG rapide comme un SSE2 xorshift+, ou si vous chronométrez juste la boucle de comptage, alors utilisez ce que vous voulez, comme rand().

Une façon d'accélérer cette opération de manière significative, même sans AVX, est de diviser les données en blocs de 255 éléments maximum, et d'accumuler le nombre de bits par octet dans des blocs ordinaires. uint64_t ordinaires. Comme les données sources ont 64 bits, nous avons besoin d'un tableau de 8 accumulateurs par octet. Le premier accumulateur compte les bits dans les positions 0, 8, 16, .... 56, le second accumulateur compte les bits en positions 1, 9, 17, .... 57; et ainsi de suite. Après avoir fini de traiter un bloc de données, nous transférons les comptes de l'accumulateur par octet dans l'accumulateur target comptes. Une fonction pour mettre à jour le target pour un bloc de 255 nombres maximum peut être codée de manière simple selon la description ci-dessus, où BITS est le nombre de bits dans les données sources :

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

Le programme ISO-C99 complet, qui devrait pouvoir fonctionner au moins sur les plates-formes Windows et Linux, est présenté ci-dessous. Il initialise les données sources avec un PRNG, effectue une vérification de l'exactitude par rapport à l'implémentation de référence du demandeur, et évalue le code de référence et la version accélérée. Sur ma machine (Intel Xeon E3-1270 v2 @ 3.50 GHz), lorsque compilé avec MSVS 2010 à pleine optimisation (/Ox), la sortie du programme est :

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

ref fait référence à la solution originale de l'auteur de la demande. L'accélération ici est d'environ un facteur 74x. Des accélérations différentes seront observées avec d'autres compilateurs (et notamment les plus récents).

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

notes et commentaires

Si ce message vous a été utile, il serait très utile que vous le partagiez avec d'autres programmeurs de cette manière, vous contribuez à diffuser nos informations.



Utilisez notre moteur de recherche

Ricerca
Generic filters

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.