GCC 部分 builtin 函数轮子

Posted on by xris

某个空虚寂寞的晚上，随便造点轮子，只是常用的二进制位运算的函数而已。实现不一定是常数最优的，可读性优先。

裸的打表、分块打表、不到 log 级别的表不写，不屑。

Counting Trailing Zeros

Counting Trailing Zeros

__builtin_ctz ，返回从最低位开始的连续的 0 的个数；如果传入 0 则行为未定义。

_BitScanForward ， Visual Studio 中的内建函数，等价于 GCC 的 __builtin_ctz 。

朴素

int ctz64(uint64_t x) {
  for (int i = 0; i != 64; ++i)
    if (x >> i & 1) return i;
  return 0;
}

二分

int ctz64(uint64_t x) {
  int r = 63;
  x &= ~x + 1;
  if (x & 0x00000000FFFFFFFF) r -= 32;
  if (x & 0x0000FFFF0000FFFF) r -= 16;
  if (x & 0x00FF00FF00FF00FF) r -= 8;
  if (x & 0x0F0F0F0F0F0F0F0F) r -= 4;
  if (x & 0x3333333333333333) r -= 2;
  if (x & 0x5555555555555555) r -= 1;
  return r;
}

浮点数

支持 IEEE-754 binary64 的架构下可以使用，代码中 double 需要保证为 binary64 。不过最低要求其实只需要使用 binary32 就可以达到目标。

不解释，自己复习 IEEE-754 定义的浮点数二进制表示。

int ctz64(uint64_t x) {
  union {
    double f;
    uint64_t i;
  } v = {.f = x & ~x + 1};
  return (v.i >> 52) - 1023;
}

哈希查表

关于代码中出现的魔数请参考 de Bruijn 序列。

int ctz64(uint64_t x) {
  static const int bit_perm[64] = {
     0,  1,  2,  7,  3, 13,  8, 19,  4, 25, 14, 28,  9, 34, 20, 40,
     5, 17, 26, 38, 15, 46, 29, 48, 10, 31, 35, 54, 21, 50, 41, 57,
    63,  6, 12, 18, 24, 27, 33, 39, 16, 37, 45, 47, 30, 53, 49, 56,
    62, 11, 23, 32, 36, 44, 52, 55, 61, 22, 43, 51, 60, 42, 59, 58,
  };
  return bit_perm[0x218A392CD3D5DBF * (x & ~x + 1) >> 58];
}

Find First Set

__builtin_ffs ，返回最低的非 0 位的下标，下标从 1 还是计数；如果传入 0 则返回 0 。

除 0 以外，满足恒等式 __builtin_ffs(x) == __builtin_ctz(x) + 1 ，懒得罗嗦了。

int ffs64(uint64_t x) {
  if (x == 0) return 0;
  return ctz64(x) + 1;
}

Count Leading Zeros

__builtin_clz ，返回从最高位开始连续的 0 的个数。当传入 0 则行为未定义。

_BitScanReverse ， Visual Studio 中的内建函数，等价于 __builtin_clz 。

朴素

int clz64(uint64_t x) {
  for (int i = 0; i != 64; ++i)
    if (x >> 63 - i & 1) return i;
}

二分

int clz64(uint64_t x) {
  int r = 0;
  if (!(x & 0xFFFFFFFF00000000)) r += 32, x <<= 32;
  if (!(x & 0xFFFF000000000000)) r += 16, x <<= 16;
  if (!(x & 0xFF00000000000000)) r += 8,  x <<= 8;
  if (!(x & 0xF000000000000000)) r += 4,  x <<= 4;
  if (!(x & 0xC000000000000000)) r += 2,  x <<= 2;
  if (!(x & 0x8000000000000000)) r += 1,  x <<= 1;
  return r;
}

浮点数

需要 IEEE-754 binary128 才能实现。实际上最小上界是尾数能完全存下 64-bit 的浮点数。 x86 的 extended 恰好差一位，但是如果我们截断尾部一位而手动判断 0 或 1 的情况的话， double 类型还是可以满足需求的。

以下为了方便，还是使用 binary128 ，假设 __float128 即为该类型。

int clz64(uint64_t x) {
  union {
    __float128 f;
    uint64_t i[2];
  } v = {.f = x};
  return 16446 - (v.i[1] >> 48);
}

利用了浮点数的指数部分，和 ctz 浮点数做法本质相同。

哈希查表

关于代码中出现的魔数请参考 de Bruijn 序列。

int clz64(uint64_t x) {
  static const int bit_perm[64] = {
    63, 62, 61, 56, 60, 50, 55, 44, 59, 38, 49, 35, 54, 29, 43, 23,
    58, 46, 37, 25, 48, 17, 34, 15, 53, 32, 28,  9, 42, 13, 22,  6,
     0, 57, 51, 45, 39, 36, 30, 24, 47, 26, 18, 16, 33, 10, 14,  7,
     1, 52, 40, 31, 27, 19, 11,  8,  2, 41, 20, 12,  3, 21,  4,  5,
  };
  x |= x >> 32;
  x |= x >> 16;
  x |= x >> 8;
  x |= x >> 4;
  x |= x >> 2;
  x |= x >> 1;
  x ^= x >> 1;
  return bit_perm[0x218A392CD3D5DBF * x >> 58];
}

Count Leading Redundant Sign Bits

__builtin_clrsb ，返回最高位开始的连续的和符号位相同的位数。

懒得罗嗦了。

int clrsb64(int64_t x) {
  if (!(x ^ x >> 63)) return 63;
  return clz64(x ^ x >> 63);
}

Population Count

__builtin_popcount ，计算二进制表示中 1 的数量。

比较常用，人们研究的比较多，算法比较多，不少骚优化。

朴素

int popcount64(uint64_t x) {
  int r = 0;
  do r += x & 1;
  while (x >>= 1);
  return r;
}

一个优化

int popcount64(uint64_t x) {
  int r = 0;
  for (; x; x &= x - 1) ++r;
  return r;
}

二分

给出最易懂的代码，以下代码具有较大优化潜能：

int popcount64(uint64_t x) {
  x = (x & 0x5555555555555555) + (x >> 1  & 0x5555555555555555);
  x = (x & 0x3333333333333333) + (x >> 2  & 0x3333333333333333);
  x = (x & 0x0F0F0F0F0F0F0F0F) + (x >> 4  & 0x0F0F0F0F0F0F0F0F);
  x = (x & 0x00FF00FF00FF00FF) + (x >> 8  & 0x00FF00FF00FF00FF);
  x = (x & 0x0000FFFF0000FFFF) + (x >> 16 & 0x0000FFFF0000FFFF);
  x = (x & 0x00000000FFFFFFFF) + (x >> 32 & 0x00000000FFFFFFFF);
  return x;
}

一个优化是在 x = (x & 0x00ff00ff00ff00ff) + (x >> 8 & 0x00ff00ff00ff00ff) 这一步过后直接 % 255 就是结果。因为 m^x\equiv1\pmod{m-1} ，所以 m 进制的数字 \sum d_km^k\equiv\sum_k d_k\pmod{m-1}。这里 m 是幂指数以幂次增长的，而 256 是最小的大于 64 的满足条件的 m 。

类似的，在算到第四轮之后完全不需要维护更高位的部分，每次只需要低位部分加上高位部分，最后截取低 8 位即可。该变种仍可继续优化，最后求和的部分可以一个乘法达成：

int popcount64(uint64_t x) {
  x -= (x >> 1) & 0x5555555555555555;
  x = (x & 0x3333333333333333) + (x >> 2 & 0x3333333333333333);
  x = x + (x >> 4) & 0x0F0F0F0F0F0F0F0F;
  return x * 0x0101010101010101 >> 56;
}

Parity

__builtin_parity ，等价于 __builtin_popcount(x) & 1 。

偷懒

int parity64(uint64_t x) {
  return popcount64(x) & 1;
}

二分 

int parity64(uint64_t x) {
  x ^= x >> 1;
  x ^= x >> 2;
  x ^= x >> 4;
  x ^= x >> 8;
  x ^= x >> 16;
  x ^= x >> 32;
  return x & 1;
}

一个优化是部分打表，在 x 的值比 16 小了之后直接用表即可。由于 parity 只会返回 01 ，因此可以压位。

参考资料

No comments yet

发表评论