GCC 部分内建位运算函数轮子

Table of Contents

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

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

Counting Trailing Zeros

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

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

朴素

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

二分

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
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;
}

这里用到了一个位运算技巧 x & ~x + 1 可以直接提取出 x 的最低设置位(Lowest Set Bit)。这是一个比较巧妙的方法,在 CTZ 中被广泛使用。它在树状数组中被也应用于定位前继节点(向前走,求前缀和)或者父节点(向后走,向上更新状态)。一个小细节是在补码表示的情况下,~x + 1-x 是等价的,本文用 ~x + 1 更清晰。

后面几行就是二分判断最低设置位在二进制表达的哪个位置。第 4 行判断是前半截还是后半截,第 5 行将 64 位二进制数分成两个 32 位,判断在 32 位的前半截还是后半截,依次类推。顺序反过来或者从 0 开始往上加都是可以的。

IEEE-754 性质

1
2
3
4
5
6
7
8
9
int ctz64(uint64_t x) {
  union {
    double f;
    uint64_t i;
  } v = {
    .f = x & ~x + 1,
  };
  return (v.i >> 52) - 1023;
}

利用了 IEEE-754 中定义的 binary 格式,在遇到 2 的幂次时,它指数域的值减去置偏值就是 2 右上那个指数,也就是我们需要求的值。上面的代码里没有管符号位,因为我们用的是 uint64_t,如果不是无符号整数可能要针对负数进行处理一下。

de Bruijn 序列性质(映射查表)

1
2
3
4
5
6
7
8
9
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[0x218A392CD3D5DBFUL * (x & ~x + 1) >> 58];
}

这种方法本质上是找到了一个巧妙的映射函数,利用了 de Bruijn 序列的一个性质将 64 个 2 的幂 \(\displaystyle\left\{2^k\middle|k\in[0,64)\cap\mathbb Z\right\}\) 双射回 \([0,64)\cap\mathbb Z\) 区间内,结合之前的 x & x + 1 技巧,我们就可以直接得到答案。

代码中的魔数 0x218A392CD3D5DBFUL 就是 B(2, 6) 的 de Bruijn 序列,如果将它转化为二进制并且头尾相接成环,使用 6 位的窗口滑动过去,可以发现每个窗口都不同。所以将它左移 0~63 位,每次都能拿到不同的 6 位数字(无论取哪部分的 6 位)。这样的话只需要获取到 x 的最低设置位,因为它本身是一个 2 的幂,它乘以任何数都等价于左移一个 CTZ。又因为是双射,所以肯定是覆盖满值域定义域的,只需要把 0~63 对应 2 的几次幂全部求出来,存到数组里再映射就行。

关于 de Bruijn 具体可以参考维基百科 de Bruijn sequence,我这里将生成 128 位 de Bruijn 序列和它的双射表的代码放上来,64 位改个参数就行:GitHub Gists

Find First Set

__builtin_ffs ,返回最低的非 0 位的下标,下标从 1 开始计数;如果传入 0 则返回 0 。

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

1
2
3
4
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

一个数学特性是 \(\mathrm{clz}\ x+\lfloor\log_2 x\rfloor+1=类型位宽\)。

朴素

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

二分

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

这个还是比较好理解的,不多解释了。

IEEE-754 性质

因为 IEEE-754 会自动调整成科学计数法,所以直接将数字转为浮点数就可以使用了。可能需要担心的就是非规格化浮点数,但是我们处理的都是整数,碰不到。

1
2
3
4
5
6
7
8
9
int clz64(uint64_t x) {
  union {
    double f;
    uint64_t i;
  } v = {
    .f = x,
  };
  return 1086 - (v.i >> 52);
}

首先也是拿到指数域,减去置偏值,得到的是 \(\left\lfloor\log_2 x\right\rfloor\),我们要求前缀 0 个数,还要用 63 减一下,也就是 63 - ((v.i >> 52) - 1023) 化简后就是上面的。

同样这里也没有考虑负数的符号位,如果有负数直接返回 0 好了。

de Bruijn 序列性质(映射查表)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
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];
}

和 ctz 本质是一样的,先单独拎出最高设置位,中间那一段 |= 就是干这个的,逻辑是把最高设置位右边填满 1,然后减一下,就只剩最高设置位了。

其它的逻辑和 de Bruijn 是一样的,代码生成改一下之前的脚本就可以了,这里就略过了。

Count Leading Redundant Sign Bits

__builtin_clrsb,返回从最高位(MSB)开始和符号位相同的位数。

1
2
3
4
int clrsb64(int64_t x) {
  if (x == 0 || x == -1) return 63;
  return clz64(x ^ x >> 63);
}

很简单,利用有符号整数的右移运算符会用符号填充高位的性质,如果本来是负数那么 x >> 63 就是 -1,如果本来是非负数那就是 0。然后和原数异或一次,如果原数是负数那所有 1 被翻转为 0,如果原数是正数那保持不变。最后求一次 CLZ 即可。

Population Count

__builtin_popcount 统计二进制位中 1 的个数。

比较常用,研究比较多,骚优化不少。

朴素

1
2
3
4
5
6
int popcount64(uint64_t x) {
  int r = 0;
  do r += x & 1;
  while (x >>= 1);
  return r;
}

一个优化

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

x & x - 1 的作用是将最低设置位置 0。很好理解,x - 1 在遇到最低位是 0 的时候会不停向上借位,直到碰到一个 1 也就是最低设置位为止。最后结果最低设置位以及更低位都被翻转,然后和原数按位与一下就把最低设置位置 0 了。

这样迭代次数和参数二进制表达中 1 的个数相同,稍微有些优化。

二分

1
2
3
4
5
6
7
8
9
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;
}

思路是这样,先把二进制按 2 位分组,每个 0x5 即 0b0101 就是分组结果,然后将这每两位中的 0 部分和 1 部分相加。再将加起来的结果按 4 位分组,每个 0x3 即 0b0011 就是分组结果,再将这 4 位中两位 0 部分代表数字和两位 1 部分代表数字相加。依此类推。

或者可以这么理解,先将输入看作 4 进制数字,然后将高位和低位相加;再将相加的结果看作 16 进制数字,然后将高位和低位相加;再将相加的结果看作 256 进制数字,然后将高位和低位相加;依次类推。之所以这么理解,是为了下面的一个优化:

1
2
3
4
5
6
int popcount64(uint64_t x) {
  x = (x & 0x5555555555555555) + (x >> 1 & 0x5555555555555555);
  x = (x & 0x3333333333333333) + (x >> 2 & 0x3333333333333333);
  x = (x & 0x0F0F0F0F0F0F0F0F) + (x >> 4 & 0x0F0F0F0F0F0F0F0F);
  return x % 255;
}

这个优化是在第 4 行这一步过后直接 mod 255 就是结果。因为 \(m^k\equiv1\pmod{m-1}\) ,所以 \(m\) 进制的数字 \(\overline{a_na_{n-1}\cdots a_3a_2a_1}\) 有 \(\sum\limits_{1\le k\le n}a_km^k\equiv\sum\limits_{1\le k\le n}a_k\pmod{m-1}\) 就是按位和。上面已经分析过第三步之后就已经拿到一个 256 进制了,所以直接 mod 255 就是答案。而更上一步则需要 mod 16,这肯定会有问题,因为 64 bits 的 popcount 最高是 64,mod 16 显然会出错。

类似的,在算到第四轮的时候完全不需要维护更高位的部分,因为单个 16 进制已经可以保存下当前完整的 1 计数数了。只需要直接低位部分加上高位部分,最后截取低 8 位即可。该变种仍可继续优化,最后求和的部分可以一个乘法达成:

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

第 2 行因为是按照 2 bits 分组的,所以它的含义,我们只需要枚举 0~4 的映射就可以了,可以看出来和原来是等价的。这么做是为了减少一次运算。

from into
00 0
01 1
10 1
11 2

第三行没变,第四行就是去掉了维护更高位部分的逻辑,上面说过 16 进制已经存得下了。而最后一行乘以 0x0101010101010101 是将每个 0xF 分组部分全部累加到最高的 8 位中,然后直接舍弃低位截取最高字节(右移 56 位)即可。

Parity

__builtin_parity 就是求 __builtin_popcount 的奇偶性。

偷懒

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

二分

1
2
3
4
5
6
7
8
9
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 只会返回 0 或 1,因此可以压位。

参考资料