返回64位整数中所有设置位的位置的最快方法是什么?

我需要一个快速的方法来获得一个64位整数的所有位的位置。 例如,给定x = 123703 ,我想填充一个数组idx[] = {0, 1, 2, 4, 5, 8, 9, 13, 14, 15, 16} x = 123703 idx[] = {0, 1, 2, 4, 5, 8, 9, 13, 14, 15, 16} 。 我们可以假设我们知道比特数的先验。 这将被称为10 ^ 12 – 10 ^ 15倍,所以速度是至关重要的。 到目前为止,我提出的最快答案是下面的怪异问题,它使用64位整数的每个字节作为表中的索引,给出该字节中设置的位数和位置:

 int64_t x; // this is the input unsigned char idx[K]; // this is the array of K bits that are set unsigned char *dst=idx, *src; unsigned char zero, one, two, three, four, five; // these hold the 0th-5th bytes zero = x & 0x0000000000FFUL; one = (x & 0x00000000FF00UL) >> 8; two = (x & 0x000000FF0000UL) >> 16; three = (x & 0x0000FF000000UL) >> 24; four = (x & 0x00FF00000000UL) >> 32; five = (x & 0xFF0000000000UL) >> 40; src=tab0+tabofs[zero ]; COPY(dst, src, n[zero ]); src=tab1+tabofs[one ]; COPY(dst, src, n[one ]); src=tab2+tabofs[two ]; COPY(dst, src, n[two ]); src=tab3+tabofs[three]; COPY(dst, src, n[three]); src=tab4+tabofs[four ]; COPY(dst, src, n[four ]); src=tab5+tabofs[five ]; COPY(dst, src, n[five ]); 

其中COPY是一个switch语句,可复制多达8个字节, n是在一个字节中设置的位数的数组, tabofs给出了tabofs的偏移量, tabX保存了第X个字节中设置位的位置。 这比使用Xeon E5-2609上的__builtin_ctz()展开的基于循环的方法要快3倍。 (见下面。)我目前正在按照字典顺序迭代x给定的位数。

有没有更好的办法?

编辑 :增加了一个例子(我后来修复)。 完整的代码可在这里: http : //pastebin.com/79X8XL2P 。 注意:带有-O2的GCC似乎优化了它,但是英特尔的编译器(我用它来编写它)不是…

另外,让我给出一些额外的背景来解决下面的一些评论。 目标是从N个可能的解释variables的宇宙中对K个variables的每个可能子集进行统计检验; 目前的具体目标是N = 41,但是我可以看到有些项目需要N达到45-50。 testing基本上涉及分解相应的数据子matrix。 在伪代码,这样的事情:

 double doTest(double *data, int64_t model) { int nidx, idx[]; double submatrix[][]; nidx = getIndices(model, idx); // get the locations of ones in model // copy data into submatrix for(int i=0; i<nidx; i++) { for(int j=0; j<nidx; j++) { submatrix[i][j] = data[idx[i]][idx[j]]; } } factorize(submatrix, nidx); return the_answer; } 

我为Intel Phi开发板编写了一个版本,在大约15天内完成N = 41个案例,其中5-10%的时间花在了一个天真的getIndices() ,版本可以节省一天或更多。 我正在为NVidia Kepler开发一个实现,但是不幸的是我的问题(可笑的数字小matrix运算)并不理想地适用于硬件(非常大的matrix运算)。 也就是说, 本文提出了一个解决scheme,通过积极展开循环和执行寄存器中的整个因式分解,似乎在我的大小的matrix上实现了数百个GFLOPS / s,但需要在编译时定义matrix的维数。 (这个循环展开应该有助于减less开销,并改善向量化的Phi版本,所以getIndices()将变得更加重要!)所以现在我想我的内核应该看起来更像:

 double *data; // move data to GPU/Phi once into shared memory template<unsigned int K> double doTestUnrolled(int *idx) { double submatrix[K][K]; // copy data into submatrix #pragma unroll for(int i=0; i<K; i++) { #pragma unroll for(int j=0; j<K; j++) { submatrix[i][j] = data[idx[i]][idx[j]]; } } factorizeUnrolled<K>(submatrix); return the_answer; } 

Phi版本从model = 0到2 ^ N(或者testing的一个子集)在`cilk_for'循环中解决了每一个模型,但是现在为了批处理GPU并分摊内核启动开销,我必须按照字典顺序对K = 1到41位中的每一个进行迭代(如doynax所述)。

编辑2:现在假期结束了,这里是我的Xeon E5-2602使用icc版本15的一些结果。我用于基准testing的代码是在这里: http : //pastebin.com/XvrGQUat 。 我对具有完全K位的整数执行位提取,因此在下表中的“基本”列中测量的字典迭代会有一些开销。 这些进行2 ^ 30次,N = 48(必要时重复)。

“CTZ”是一个使用gcc内部__builtin_ctzll获取最低位设置的循环:

 for(int i=0; i<K; i++) { idx[i] = __builtin_ctzll(tmp); lb = tmp & -tmp; // get lowest bit tmp ^= lb; // remove lowest bit from tmp } 

Mark是Mark的无路循环:

 for(int i=0; i<K; i++) { *dst = i; dst += x & 1; x >>= 1; } 

Tab1是我原始的基于表格的代码,具有以下副本macros:

 #define COPY(d, s, n) \ switch(n) { \ case 8: *(d++) = *(s++); \ case 7: *(d++) = *(s++); \ case 6: *(d++) = *(s++); \ case 5: *(d++) = *(s++); \ case 4: *(d++) = *(s++); \ case 3: *(d++) = *(s++); \ case 2: *(d++) = *(s++); \ case 1: *(d++) = *(s++); \ case 0: break; \ } 

Tab2是与Tab1相同的代码,但复制macros只是作为一个副本移动8个字节(从doynax和LưuVĩnhPhúc的想法…但注意这不能确保alignment):

 #define COPY2(d, s, n) { *((uint64_t *)d) = *((uint64_t *)s); d+=n; } 

这是结果。 我猜我最初声称Tab1比CTZ快3倍只适用于大K(我正在testing的地方)。 马克的循环比我的原始代码快,但是摆脱COPY2macros中的分支需要K> 8的蛋糕。

  K Base CTZ Mark Tab1 Tab2 001 4.97s 6.42s 6.66s 18.23s 12.77s 002 4.95s 8.49s 7.28s 19.50s 12.33s 004 4.95s 9.83s 8.68s 19.74s 11.92s 006 4.95s 16.86s 9.53s 20.48s 11.66s 008 4.95s 19.21s 13.87s 20.77s 11.92s 010 4.95s 21.53s 13.09s 21.02s 11.28s 015 4.95s 32.64s 17.75s 23.30s 10.98s 020 4.99s 42.00s 21.75s 27.15s 10.96s 030 5.00s 100.64s 35.48s 35.84s 11.07s 040 5.01s 131.96s 44.55s 44.51s 11.58s 

我相信这里性能的关键是把重点放在更大的问题上,而不是从随机整数中微观提取位位置。

根据您的示例代码和之前的SO问题,您将按顺序设置K位的所有单词,并从中提取位索引。 这大大简化了事情。

如果是这样的话,而不是重build每个迭代的位位置尝试直接递增位arrays中的位置。 有一半时间这将涉及单循环迭代和增量。

沿着这些线路的东西:

 // Walk through all len-bit words with num-bits set in order void enumerate(size_t num, size_t len) { size_t i; unsigned int bitpos[64 + 1]; // Seed with the lowest word plus a sentinel for(i = 0; i < num; ++i) bitpos[i] = i; bitpos[i] = 0; // Here goes the main loop do { // Do something with the resulting data process(bitpos, num); // Increment the least-significant series of consecutive bits for(i = 0; bitpos[i + 1] == bitpos[i] + 1; ++i) bitpos[i] = i; // Stop on reaching the top } while(++bitpos[i] != len); } // Test function void process(const unsigned int *bits, size_t num) { do printf("%d ", bits[--num]); while(num); putchar('\n'); } 

不是特别优化,但你得到了一般的想法。

这里有一些非常简单的东西可能会更快 – 没有testing就无法知道。 很大程度上取决于设置的位数与未设置的数量。 你可以展开这个去掉分支,但是对于今天的处理器,我不知道它是否会加速。

 unsigned char idx[K+1]; // need one extra for overwrite protection unsigned char *dst=idx; for (unsigned char i = 0; i < 50; i++) { *dst = i; dst += x & 1; x >>= 1; } 

PS您的示例输出的问题是错误的,请参阅http://ideone.com/2o032E

作为最小的修改:

 int64_t x; char idx[K+1]; char *dst=idx; const int BITS = 8; for (int i = 0 ; i < 64+BITS; i += BITS) { int y = (x & ((1<<BITS)-1)); char* end = strcat(dst, tab[y]); // tab[y] is a _string_ for (; dst != end; ++dst) { *dst += (i - 1); // tab[] is null-terminated so bit positions are 1 to BITS. } x >>= BITS; } 

BITS的select决定了表的大小。 8,13和16是合乎逻辑的select。 每个条目都是一个string,零终止并包含具有1个偏移量的位位置。 即选项卡[5]是"\x03\x01" 。 内部循环修复了这个偏移量。

稍微高效一点:用strcat和inner loop取代

 char const* ptr = tab[y]; while (*ptr) { *dst++ = *ptr++ + (i-1); } 

如果循环包含分支,循环展开可能会有点痛苦,因为复制这些分支语句不利于分支预测器。 我会很乐意把这个决定留给编译器。

我正在考虑的一件事是tab[y]是一个指向string的指针数组。 这些高度相似: "\x1""\x3\x1"的后缀。 实际上,每个不以"\x8"开头的string都是一个string的后缀。 我想知道你需要多less个独特的string,以及tab[y]实际上需要多less程度。 例如通过上面的逻辑, tab[128+x] == tab[x]-1

[编辑]

没关系,你肯定需要128个以"\x8"开始的标签条目,因为它们永远不是另一个string的后缀。 尽pipe如此, tab[128+x] == tab[x]-1规则意味着您可以保存一半的条目,但是需要两个额外的指令: char const* ptr = tab[x & 0x7F] - ((x>>7) & 1) 。 (设置tab[]指向\x8

使用char不会帮助你提高速度,但实际上在计算时经常需要更多的ANDing和符号/零延伸。 只有在应该适合caching的非常大的数组的情况下,应该使用较小的inttypes

你可以改进的另一件事是COPYmacros。 如果可能,请复制整个单词,而不是逐字节复制

 inline COPY(unsigned char *dst, unsigned char *src, int n) { switch(n) { // remember to align dst and src when declaring case 8: *((int64_t*)dst) = *((int64_t*)src); break; case 7: *((int32_t*)dst) = *((int32_t*)src); *((int16_t*)(dst + 4)) = *((int32_t*)(src + 4)); dst[6] = src[6]; break; case 6: *((int32_t*)dst) = *((int32_t*)src); *((int16_t*)(dst + 4)) = *((int32_t*)(src + 4)); break; case 5: *((int32_t*)dst) = *((int32_t*)src); dst[4] = src[4]; break; case 4: *((int32_t*)dst) = *((int32_t*)src); break; case 3: *((int16_t*)dst) = *((int16_t*)src); dst[2] = src[2]; break; case 2: *((int16_t*)dst) = *((int16_t*)src); break; case 1: dst[0] = src[0]; break; case 0: break; } 

另外,由于tabofs [x]和n [x]通常是彼此靠近的,请尝试将它放在内存中以确保它们始终同时处于caching中

 typedef struct TAB_N { int16_t n, tabofs; } tab_n[256]; src=tab0+tab_n[b0].tabofs; COPY(dst, src, tab_n[b0].n); src=tab0+tab_n[b1].tabofs; COPY(dst, src, tab_n[b1].n); src=tab0+tab_n[b2].tabofs; COPY(dst, src, tab_n[b2].n); src=tab0+tab_n[b3].tabofs; COPY(dst, src, tab_n[b3].n); src=tab0+tab_n[b4].tabofs; COPY(dst, src, tab_n[b4].n); src=tab0+tab_n[b5].tabofs; COPY(dst, src, tab_n[b5].n); 

最后但并非最不重要的是, gettimeofday不是用于性能计数。 使用QueryPerformanceCounter,而是更精确

您的代码使用1个字节(256个条目)索引表。 如果使用2个字节(65536个条目)的索引表,则可以加快系数2。

不幸的是,你可能无法进一步扩展 – 对于3字节的表大小将是16MB,不太可能适合CPU本地caching,它只会让事情变慢。

问题是你将如何收集职位?
如果你不得不多次迭代它,那么是的,像现在这样收集它们可能会很有趣,并且迭代很多。
但是,如果只迭代一次或几次,那么您可能会考虑不创build一个中间数组的位置,而只是在每个遇到的位1上调用一个处理块闭包/函数,同时在位上进行迭代。

下面是我在Smalltalk中写的位迭代器的一个简单的例子:

 LargePositiveInteger>>bitsDo: aBlock | mask offset | 1 to: self digitLength do: [:iByte | offset := (iByte - 1) << 3. mask := (self digitAt: iByte). [mask = 0] whileFalse: [aBlock value: mask lowBit + offset. mask := mask bitAnd: mask - 1]] 

LargePositiveInteger是由字节数组成的任意长度的整数。 lowBit回答最低位的等级,并且被实现为具有256个条目的查找表。

在C ++ 2011中,您可以轻松地传递闭包,所以应该很容易翻译。

 uint64_t x; unsigned int mask; void (*process_bit_position)(unsigned int); unsigned char offset = 0; unsigned char lowBitTable[16] = {0,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0}; // 0-based, first entry is unused while( x ) { mask = x & 0xFUL; while (mask) { process_bit_position( lowBitTable[mask]+offset ); mask &= mask - 1; } offset += 4; x >>= 4; } 

这个例子是用一个4位的表格来演示的,但是如果它适合caching,你可以很容易地把它扩展到13位或更多。

对于分支预测,内部循环可以用一个额外的表nbit=numBitTable[mask]重写为for(i=0;i<nbit;i++) ,然后用开关展开(编译器可以这样做?让你先量一下它的performance

这被发现太慢了吗?
小而粗糙,但全部在caching和CPU寄存器中;

 void mybits(uint64_t x, unsigned char *idx) { unsigned char n = 0; do { if (x & 1) *(idx++) = n; n++; } while (x >>= 1); // If x is signed this will never end *idx = (unsigned char) 255; // List Terminator } 

它仍然是展开循环的3倍,并产生一个64真/假值的数组(这是不是很想要的)

 void mybits_3_2(uint64_t x, idx_type idx[]) { #define SET(i) (idx[i] = (x & (1UL<<i))) SET( 0); SET( 1); SET( 2); SET( 3); ... SET(63); } 

这里有一些紧密的代码,写成1字节(8位),但它应该很容易,显然扩展到64位。

 int main(void) { int x = 187; int ans[8] = {-1,-1,-1,-1,-1,-1,-1,-1}; int idx = 0; while (x) { switch (x & ~(x-1)) { case 0x01: ans[idx++] = 0; break; case 0x02: ans[idx++] = 1; break; case 0x04: ans[idx++] = 2; break; case 0x08: ans[idx++] = 3; break; case 0x10: ans[idx++] = 4; break; case 0x20: ans[idx++] = 5; break; case 0x40: ans[idx++] = 6; break; case 0x80: ans[idx++] = 7; break; } x &= x-1; } getchar(); return 0; } 

输出数组应该是:

 ans = {0,1,3,4,5,7,-1,-1}; 

如果我拿“我需要一个快速的方法来获得在一个64位整数的所有位”位置字面上…

我意识到这是几个星期的老,但是出于好奇,我记得在我的组装date间,用CBM64和Amiga使用算术移位,然后检查进位标志 – 如果设置了,那么移位位是1,如果清除那么它是零

例如,对于算术左移(从位64到位0检查)….

 pseudo code (ignore instruction mix etc errors and oversimplification...been a while): move #64+1, counter loop. ASL 64bitinteger BCS carryset decctr. dec counter bne loop exit carryset. //store #counter-1 (ie bit position) in datastruct indexed by counter jmp decctr 

…我希望你明白了。

从那时起,我还没有使用过程序集,但是我想知道是否可以使用类似于上面的一些C ++内嵌程序集来做类似的事情。 我们可以在汇编中完成整个转换(很less的代码行),build立一个合适的数据结构。 C ++可以简单地检查答案。

如果这是可能的,那么我想它会很快。

假设置位数的稀疏性,

 int count = 0; unsigned int tmp_bitmap = x; while (tmp_bitmap > 0) { int next_psn = __builtin_ffs(tmp_bitmap) - 1; tmp_bitmap &= (tmp_bitmap-1); id[count++] = next_psn; }