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

Question

我需要一种快速获取64位整数中所有位为1的位置的方法。例如，给定x = 123703，我想填充一个数组idx[] = {0, 1, 2, 4, 5, 8, 9, 13, 14, 15, 16}。我们可以假设我们事先知道位数。这将被称为10¹² - 10¹⁵倍，因此速度至关重要。到目前为止，我想到的最快的答案是以下的怪物代码，它使用64位整数的每个字节作为索引，并在其中存储每个字节中位为1的数量和位置的表：

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 个字节，n 是一个字节数组，表示一个字节中设置的位数，tabofs 给出了偏移量，指向 tabX 中持有的 X-th 字节中设置的位的位置。这比在我的 Xeon E5-2609 上使用具有 __builtin_ctz() 的展开循环方法快约 3 倍。（见下文。）我目前正在按字典序遍历给定位数的 x。

还有更好的方法吗？

编辑：添加了一个示例（我随后已经修复）。完整代码可在此处找到：http://pastebin.com/79X8XL2P。注意：GCC 使用 -O2 似乎会优化掉它，但英特尔编译器（我用来编写它）不会……

此外，让我提供一些额外的背景信息，以回答以下评论。目标是对可能的解释变量 N 中的每个可能子集执行统计测试；目前的特定目标是 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;
}

我已经编写了一个版本，在Intel Phi板上运行，应该能在大约15天内完成N=41的情况，其中约5-10%的时间花费在朴素的getIndices()函数中，所以优化后的版本可以节省一天或更多时间。我正在为NVidia Kepler实现一个版本，但不幸的是，我的问题（大量小矩阵操作）并不完全适合硬件（大量大矩阵操作）。尽管如此，这篇论文提出了一种解决方案，似乎可以在与我相同大小的矩阵上实现数百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版本通过`cilk_for`循环解决模型，循环从model=0到2的N次方(或者为了测试而选定的一个子集)。但是，现在需要对每个K=1到41个比特设置的模型号进行按字典顺序迭代，以便为GPU批处理工作并分摊内核启动开销，正如doynax所指出的那样。

编辑2: 现在假期结束了，我使用icc 15版本在我的Xeon E5-2602上进行了一些基准测试。我用来进行基准测试的代码在这里：http://pastebin.com/XvrGQUat。我对恰好设置了K个比特位的整数执行位提取，因此在表格中"Base"列中测量了按字典顺序迭代的一些开销。这些操作将在N=48时重复执行2的30次方次(如有必要)。

"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是我的基于表格的原始代码，其中包含以下拷贝宏：

#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（我进行测试时）。马克的循环比我的原始代码更快，但是对于K > 8，在COPY2宏中消除分支是最好的选择。

 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

- Andrew

---

出于好奇，这有什么样的应用？ - Cameron

难道尝试一个单一的模型不比计算这个小玩意儿更加昂贵吗？或者你手头有一些令人印象深刻的计算机集群吗？我的意思是，仅仅数到10^12对于一般的计算机来说就需要相当长的时间，更不用说10^15了... - Thomas

我并不太喜欢那些（可能是不可预测的）用于进入展开复制循环的代码分支。你试过直接发送八个完整字节，然后按照表格中的计数推进dst吗？另外，你可能可以将生成的汇编代码发布出来供我们查阅。我们已经深入到微观优化领域，很容易对代码生成做出不必要的假设。 - doynax

你知道你发布的代码只能在64位中的48位上运行，对吧？如果你的限制是50位，那么你就短了2位。 - Mark Ransom

有一个包含一系列设置了特定位的比特位列表，相较于只有“int64_t x”中设置这些位，我想知道它更加实用的原因是什么？同时，了解数据集将有助于优化（例如：大多数输入是否为0？）...并且在进行10^12次调用时，优化可能非常重要。此外，您是否查看了（man-k bit）ffs（3），ffs（3p），ffsl（3），ffsll（3）等信息？它们可以使用为此设计的CPU指令...而C语言无法与手动编写的汇编内联函数或宏相匹配。 - 9mjb

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

虽然没有特别优化，但是你可以得到一般的想法。

- doynax

---

你的观点很有道理（尽管我没有在比赛更新中提到它...） - Andrew

啊哇..你让我提心吊胆的。这可能是作弊，但至少要单独发布数字。无论如何，正如你现在可能已经猜到的那样，真正的关键是进一步探索更好的（搜索？）算法，以避免需要在第一时间暴力评估所有排列。 - doynax

明天我回到工作岗位后，我会发布数字。不幸的是，理论规定我们必须对整个幂集进行求和（任何小于此的都是近似值），而这个练习的整个目的就是量化人们提出的“更好”的算法的影响... - Andrew

Answer 2

这里有一个非常简单的方法，也许会更快 - 没有测试就没有办法知道。很大程度上取决于设置的位数与未设置的位数。你可以展开它以完全消除分支，但是随着今天处理器的发展，我不知道是否会加速。

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

附注：你在问题中提供的样例输出是错误的，请参见http://ideone.com/2o032E

- Mark Ransom

---

感谢您发现这个示例错误，我不知道当时在想什么（或者更确切地说，没有想）。我将在周一与英特尔编译器进行比较。 - Andrew

完全展开循环应该有所帮助。如果位模式非常稀疏，以下替代方案可能会很好地工作，如果快速的计算末尾零位数操作ctz()可用：int i = 0; while (x) { i += ctz(x); *dst = i; i++; dst++; x >>= i; }。类似的，基于种群计数popc()的并行解决方案似乎也是可能的，例如在GPU上实现，其中每个线程可以处理一个输入位并写入dst[]的适当元素。这基本上是一个流压缩问题。 - njuffa

抱歉，应该是 int s, i = 0; while (x) { s = ctz(x); i += s; *dst = i; dst++; i++; s++; x >>= s; }。 - njuffa

哇...你的展开无分支循环甚至在设置位数较少时也能击败ctz。太酷了！ - Andrew

Answer 3

作为一种最小的修改方案：

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 决定了表格的大小。8、13和16是合理的选择。每个条目都是一个字符串，以零结尾，并包含具有1偏移量的位位置。即 tab[5] 是 "\x03\x01"。内部循环会修正这个偏移量。

稍微更有效率一些：将 strcat 和内部循环替换为：

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

循环展开如果循环包含分支语句会有些麻烦，因为复制这些分支语句并不能帮助分支预测器。我很乐意将这个决定留给编译器。

我正在考虑的一件事是，tab[y] 是一个指向字符串的指针数组。这些字符串非常相似："\x1" 是 "\x3\x1" 的后缀。实际上，每个不以 "\x8" 开头的字符串都是以一个以 "\x8" 开头的字符串结尾。我想知道你需要多少个唯一的字符串，并且需要多少程度上实际上需要 tab[y]。例如，根据上述逻辑，tab[128+x] == tab[x]-1。

[编辑]

没关系，你绝对需要 128 个以 "\x8" 开头的 tab 条目，因为它们从不是另一个字符串的后缀。但是，tab[128+x] == tab[x]-1 规则意味着你可以节省一半的条目，但代价是两个额外的指令：char const* ptr = tab[x & 0x7F] - ((x>>7) & 1)。（设置 tab[] 指向 \x8 之后）

- MSalters

Answer 4

使用char类型无法提高速度，事实上经常需要进行更多的AND运算和符号/零扩展计算。只有在应该适合缓存的非常大的数组的情况下，才应该使用较小的int类型

你可以改进的另一件事是COPY宏。如果可能，不要逐字节复制，而是整个单词复制

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]经常同时访问，因此尝试将它们放在内存中靠近的位置，以确保它们始终同时缓存在高速缓存中。

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代替，它更加精确。

- phuclv

---

好的，我会在周一上班时使用英特尔编译器进行检查。 - Andrew

Answer 5

假设位集中的位数是稀疏的，

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

- boolAeon

---

你能解释一下为什么这是解决原问题的最快方法，以及你尝试过哪些其他方法吗？ - mjuarez

如果设定的位数是稀疏的，那么这是最快的方法，因为它使用gcc内置函数来查找下一个设置的位，它会快速跳过零位。 - Marcelo Pacheco

Answer 6

你的代码使用了1字节（256个条目）的索引表。如果你使用2字节（65536个条目）的索引表，可以将其加速2倍。

不幸的是，你可能无法进一步扩展它 - 3字节的表大小将达到16MB，不太可能适合CPU本地缓存，并且只会使事情变得更慢。

- mvp

---

@OliCharlesworth：Haswell微架构的L1缓存是每个核心64KB（http://en.wikipedia.org/wiki/Haswell_%28microarchitecture%29）-这可能足以实际良好地工作。唯一的方法是进行测试来确定。 - mvp

没有必要坚持使用字节大小的密钥。8192个条目的表也可以工作。(13位，5*13 > 64) - MSalters

@mvp：它被分成两个32k块，一个用于数据，一个用于指令。因此对于随机数据，平均情况下您将错过L1一半的时间。但我同意，这可能仍然比256个条目的表总体上更快。 - Oliver Charlesworth

Answer 7

如果我字面理解"I need a fast way to get the position of all one bits in a 64-bit integer"的话...

我知道这个问题已经几个星期了，但是出于好奇，我还记得在使用CBM64和Amiga汇编语言的时代，我们可以使用算术移位然后检查进位标志 - 如果它被设置，则移位后的位为1，否则为0

例如，对于算术左移（从第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 (i.e. bit position) in datastruct indexed by counter
    jmp decctr

...我希望你能理解这个想法。

自那以后，我就没有再使用汇编语言了，但我在想我们是否可以使用一些类似于上面的C++内联汇编来在这里做类似的事情。我们可以在汇编中完成整个转换（非常少量的代码），构建一个适当的数据结构。C++只需检查答案即可。

如果这是可能的，那么我想它应该会非常快。

- Jim

Answer 8

问题是你打算如何处理位置集合？
如果你需要多次迭代它，那么像你现在所做的一样，将它们收集起来并进行多次迭代可能是有趣的。
但是，如果只是为了迭代一次或几次，那么你可以考虑不创建一个中间的位置数组，在位迭代时遇到每个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位表演示的，但如果适合缓存，您可以轻松将其扩展到13位或更多。

对于分支预测，内部循环可以重写为for(i=0;i<nbit;i++)，并带有一个额外的表格nbit=numBitTable[mask]，然后使用switch展开（编译器可以执行此操作？），但是我让您首先衡量其性能...

- aka.nice

Answer 9

这被发现太慢了吗？
虽然小而粗糙，但所有内容都在高速缓存和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
}

将循环展开并生成一个包含64个true/false值的数组仍然比较快（但这不完全是想要的结果）。

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

- 9mjb

Answer 10

这里有一些紧凑的代码，为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};

- abelenky