43

我需要一种快速的方法来获取 64 位整数中所有一位的位置。例如，给定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是一个最多复制 8 个字节的 switch 语句，n是一个字节中设置的位数的数组，并tabofs给出了偏移量tabX，它保存了第 X 个字节中设置的位的位置。 __builtin_ctz()这比在我的 Xeon E5-2609 上展开的基于循环的方法快大约 3 倍。（见下文。）我目前正在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;
}

我为 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（或者，更确切地说，是用于测试的子集），但现在为了为 GPU 批量工作并分摊内核启动开销，我必须迭代每个 K=1 到 41 位集合的型号按字典顺序排列（如 doynax 所述）。

编辑 2： 现在假期结束了，这里是我的 Xeon E5-2602 使用 icc 版本 15 的一些结果。我用来基准测试的代码在这里：http ://pastebin.com/XvrGQUat 。我对恰好设置了 K 位的整数执行位提取，因此下表中“基本”列中测量的词典迭代存在一些开销。这些以 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 循环：

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
4

11 回答 11

7

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

从您的示例代码和之前的 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');
}

没有特别优化，但您大致了解。

于 2013-12-21T08:42:50.647 回答
6

这是一些非常简单的东西，可能会更快 - 没有测试就无法知道。很大程度上取决于设置的位数与未设置的位数。您可以展开它以完全删除分支，但对于今天的处理器，我不知道它是否会加速。

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

于 2013-12-21T03:07:42.130 回答
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". 内循环修复了这个偏移量。

效率更高一点：将strcatand 内部循环替换为

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

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

我正在考虑的一件事是tab[y]指向字符串的指针数组。它们非常相似："\x1"是 . 的后缀"\x3\x1"。实际上，每个不以开头"\x8"的字符串都是以开头的字符串的后缀。我想知道您需要多少个独特的字符串，以及tab[y]实际上需要多少。例如，按照上面的逻辑，tab[128+x] == tab[x]-1.

[编辑]

没关系，您肯定需要 128 个选项卡条目以开头，"\x8"因为它们绝不是另一个字符串的后缀。尽管如此，该tab[128+x] == tab[x]-1规则意味着您可以节省一半的条目，但代价是两个额外的指令：char const* ptr = tab[x & 0x7F] - ((x>>7) & 1). （设置tab[]到点后\x8）

于 2013-12-21T13:40:38.883 回答
2

使用 char 不会帮助您提高速度，但实际上在计算时通常需要更多的 ANDing 和符号/零扩展。只有在应该适合缓存的非常大的数组的情况下，才应该使用较小的 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，它更精确

于 2013-12-21T03:04:39.027 回答
1

您的代码使用 1 字节（256 个条目）索引表。如果使用 2 字节（65536 个条目）索引表，则可以将速度提高 2 倍。

不幸的是，您可能无法进一步扩展 - 对于 3 字节的表大小将是 16MB，不太可能适合 CPU 本地缓存，并且只会让事情变慢。

于 2013-12-21T02:18:08.123 回答
1

假设集合位数的稀疏性，

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;
}
于 2015-05-30T06:03:54.103 回答
0

问题是您将如何处理职位集合？
如果您必须对其进行多次迭代，那么是的，像现在这样收集它们一次并进行多次迭代可能会很有趣。
但是，如果它仅用于迭代一次或几次，那么您可能会考虑不创建中间位置数组，而只是在迭代位时在每个遇到的 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++)用一个附加表重写为 anbit=numBitTable[mask]然后用一个开关展开（编译器可以做到吗？），但我让你先测量它的执行情况......

于 2013-12-21T13:46:54.650 回答
0

有没有发现这太慢了？
小而粗，但都在缓存和 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 个真/假值的数组仍然快 3 倍（这不是我们想要的）

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);
}
于 2013-12-21T14:00:26.930 回答
0

如果我从字面上理解“我需要一种快速的方法来获取 64 位整数中所有一位的位置”......

我意识到这已经有几个星期了，但是出于好奇，我记得在我组装 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 (i.e. bit position) in datastruct indexed by counter
    jmp decctr

...我希望你明白。

从那时起我就没有使用过汇编，但我想知道我们是否可以使用一些类似于上面的 C++ 内联汇编在这里做类似的事情。我们可以在汇编中进行整个转换（很少的代码行），建立一个合适的数据结构。C++ 可以简单地检查答案。

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

于 2014-02-01T00:55:08.357 回答
0

这是一些紧凑的代码，为 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};
于 2013-12-22T02:53:34.413 回答
0

一个简单的解决方案，但可能不是最快的，具体取决于 log 和 pow 函数的时间：

#include<math.h>

void getSetBits(unsigned long num){

    int bit;
    while(num){
        bit = log2(num);
        num -= pow(2, bit);

        printf("%i\n", bit); // use bit number
    }
}

复杂度 O(D) | D 是设置的位数。

于 2021-07-13T15:27:42.887 回答