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

Question

我需要一种快速的方法来获取 64 位整数中所有一位的位置。例如，给定 x = 123703 , 我想填充一个数组 idx[] = {0, 1, 2, 4, 5, 8, 9, 13, 14, 15, 16} .我们可以假设我们先验地知道位数。这将被称为 1012 - 1015 次，因此速度至关重要。到目前为止，我想出的最快的答案是以下怪物，它使用 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是一个最多可复制 8 个字节的 switch 语句，n是一个字节中设置的位数的数组，tabofs给出 tabX 的偏移量，它保存第 X 个字节中设置位的位置。这比使用 __builtin_ctz() 的基于展开循环的方法快大约 3 倍。在我的至强 E5-2609 上。 (见下文。)我目前正在迭代 x对于给定数量的位集，按字典顺序排列。
有没有更好的办法？
编辑 :添加了一个示例(我随后已修复)。完整代码可在此处获得:http://pastebin.com/79X8XL2P .注意:带有 -O2 的 GCC 似乎对其进行了优化，但英特尔的编译器(我曾经编写过它)并没有......
另外，让我提供一些额外的背景来解决下面的一些评论。目标是对 N 个可能的解释变量中的 K 个变量的每个可能子集进行统计检验；现在的具体目标是 N=41，但我可以看到一些项目需要 N 高达 45-50。该测试主要涉及分解相应的数据子矩阵。在伪代码中，是这样的:

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

我为英特尔 Phi 板编写了一个版本，它应该在大约 15 天内完成 N=41 的情况，其中约 5-10% 的时间花在了一个简单的 getIndices() 上。所以立即使用更快的版本可以节省一天或更多时间。我也在研究 NVidia Kepler 的实现，但不幸的是我遇到的问题(小矩阵运算的数量可笑)并不理想地适合硬件(可笑的大矩阵运算)。也就是说，this paper提出了一个解决方案，通过积极展开循环并在寄存器中执行整个分解，似乎在我的大小的矩阵上实现了数百 GFLOPS/s，但需要注意的是，矩阵的维度在编译时定义。 (这个循环展开应该有助于减少开销并改进 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 到 2N(或者，更确切地说，是测试的子集)的 `cilk_for' 循环中解决每个模型，但现在为了为 GPU 批量工作并分摊内核启动开销，我必须迭代模型为 K=1 到 41 位集(如 doynax 指出的)中的每一个按字典顺序排列的数字。
编辑 2:现在假期结束了，下面是我的 Xeon E5-2602 使用 icc 版本 15 的一些结果。我用来进行基准测试的代码在这里:http://pastebin.com/XvrGQUat .我对恰好设置了 K 位的整数执行位提取，因此在下表的“基本”列中测量的词典迭代存在一些开销。这些执行 230 次，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 循环:

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

Tab1 是我的原始基于表格的代码，带有以下复制宏:

#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 的代码相同，但复制宏只是将 8 个字节作为单个副本移动(从 doynax 和 Lưu Vĩnh Phúc 获取想法...但请注意，这并不能确保对齐):

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

这是结果。我想我最初声称 Tab1 比 CTZ 快 3 倍仅适用于大 K(我正在测试的地方)。 Mark 的循环比我的原始代码快，但摆脱了 COPY2 中的分支宏为 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

Answer 1

我相信这里的性能关键是关注更大的问题，而不是从随机整数中提取位位置的微观优化。

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

如果是这样，那么每次迭代不是重建位位置，而是尝试直接增加位数组中的位置。一半的时间这将涉及单个循环迭代和增量。

沿着这些路线的东西:

// 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');
}

没有特别优化，但你得到了一般的想法。