位运算三则

· · 个人记录

大概是一篇介绍如何简单在 WRM 下实现 O(1) 找正整数二进制最高位和最低位，同时在 O(\log^*w) 的时间内数正整数二进制中一的数量。随着时代的发展这仨现在基本上已经都有指令能直接调了，所以这个东西基本上 useless。

typedef uint8_t u8;
typedef uint32_t u32;
typedef uint64_t u64;

我们先来解决一个比较简单的问题：给定 2^h（0\le h\lt 64），求 h。这个事情，说白了就是给你 64 个数然后要在其中高效地查找。这个问题我们一般有两种处理。第一种是直接上平衡树——很遗憾，这东西看上去不像是能优化；第二种是用一张哈希表——那么现在的问题变成了如何制作一个优秀的哈希函数。具体而言，我们现在需要找一个函数 f，满足 g(h)=f(2^h) 在 0,1,\cdots,63 上构成双射，且 f 能快速计算。一个比较简单的思路是把 x 乘上一个固定常数然后取其中连续的若干个二进制位，由于 x 是二的幂，我们本质是取了固定常数二进制中的一段。那我们要满足固定常数中能生成 0 到 63 的所有二进制表达，这不禁让我们想起一个题。具体而言，把一个 [0,64) 间的数认为是它前五位和后五位之间连的边，那我们要构造的无非就是这个图的一条欧拉回路。然后我们就做完了，取乘积的第 [58,64) 个二进制位即可。

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

解决了这个问题之后二进制最低位是送的。不会的话左转树状数组题解区。

int ut_lglbt64(u64 x) { return ut_lg64(x & -x); }

接下来我们要解决的问题是 popcount。我们先会有一个比较简单的并行思路。初始时，64 个二进制位被分为 64 块，每块长度为 1，一个块内的二进制值代表了这个块内的一的数量（显然）。每一轮迭代中，我们把相邻两段合并起来，这样我们的块长就翻了一倍。比如这两行可以把块长迭代到 2 再迭代到 4。

  x -= (x & 0xaaaaaaaaaaaaaaaaull) >> 1;
  x = (x & 0x3333333333333333ull) + ((x & 0xccccccccccccccccull) >> 2);

然而这样迭代下去的复杂度是 O(\log w) 的，不是很优。我们考虑去挖掘乘法的潜力。我们设现在的块长是 B，那我们如果制备一个二进制数，第 iB（i\in[0,t)）个二进制位为一，那么我们会发现把这个二进制数和目前的分块数乘起来的结果会导致新数的第 i 块是原数中第 (i-t,i] 的块的数的和，提取一下需要的位置就可以把块长扩大到 tB 了。然后我们考虑 t 的取值，注意到为了让卷积不溢出必须有 tB\le 2^{B}-1。假设 B=2^h 且要求 t 为二的幂，则 t\le 2^{2^h-h-1}，tB\le 2^{2^h-1}。比如我们块长为四，一轮乘法就可以把块长整到 2^{4-1}=8（好吧这里没啥用），然后 8 一轮到 2^{8-1}=128\ge 64，做完了。注意到 2^{2-1}=2，所以前两轮还是得用倍增迭代。

int ut_popcnt64(u64 x) {
  x -= (x & 0xaaaaaaaaaaaaaaaaull) >> 1;
  x = (x & 0x3333333333333333ull) + ((x & 0xccccccccccccccccull) >> 2);
  x = ((x * 0x11) & 0xf0f0f0f0f0f0f0f0ull) >> 4;
  return (x * 0x0101010101010101ull) >> 56;
}

接下来我们要解决一个数二进制最高位。考虑对 64 个二进制位按 B=\sqrt w=8 分块，我们会依次解决两个问题：判断数在哪一块，[0,2^B) 范围内的二进制最高位问题。

我们先来第一个问题，判断数在哪一块。我们把这部分分为两个部分：把所有有值的块化成某个固定形式，根据固定形式找值。我们这里的固定形式是：一个块有值当且仅当其及其之后的所有块最高位为一。

首先我们考虑怎么化形式。第一步我们需要先把一个块有值当且仅当这个块最高位为一这个事情给做出来。

方便起见我们假设所有块的最高位为 0，具体而言就是把最高位是 1 搬到第二高位上去。然后我们考虑给所有块最高位统一改成一，然后给所有块统一减一。不难发现只有整块全空才会在最高位触发退位把它置零，于是我们完成了目标。

u64 ut_blkpstv(u64 x) {
  x = (x & 0x7f7f7f7f7f7f7f7full) | (x & 0x8080808080808080ull) >> 1;
  return (x | 0x8080808080808080ull) - 0x0101010101010101ull;
}

然后我们考虑怎么进一步推到需要的形式。我们先把每一块最高位的一提取出来（其他位置置零）并把它们搬运到对应块的最低位。然后我们制备一个二进制数，满足第 iB（i\in[0,B)）个二进制位为一，那么我们会发现把这个二进制数和目前的分块数乘起来的结果会导致新数的第 i 块是原数中第 i 块及以前的前缀和，而第 i+B-1 块是原数第 i 块及以后的后缀和。注意到这里我们会让数的二进制长度几乎翻倍，于是我们最好先迭代一轮把二进制长度事先减半。另一方面，一块内的数至多为 B\lt 2^B，不会溢出。然后我们就可以用一步位运算完成等价于把每个块的数求后缀和的效果。

我们对这个新数再做一次上一步的化形式。不难发现，此时有值的块构成后缀，也就是说这个更新的的数恰好一个后缀的块最高位为一。我们完成了我们的目标。

接下来我们考虑给定化好形式的数怎么做。不难发现把每一块最高位的一搬到最低位并把其它位置零，然后制备一整块的一乘上去，得到二的幂减一，可以直接处理了。

int ut_blkfst(u64 x) {
  x = x >> 7 & 0x0101010101010101ull;
  return (ut_lg64(x * 0xff + 1) >> 3) - 1;
}

于是我们接下来只需要解决 [0,2^B) 范围内的二进制最高位问题就可以了，首先判掉最高位把值域降到 [0,2^{B-1})。我们充分利用小值域的特性，制备一个每个 B 的位置都是一的二进制数和它乘起来，这样等于是把这个块的内容平铺到整数的每个块上。然后把每个块的最高位置一并生成一个新的二进制数，满足这个数第 i 块从低到高第 i 个二进制位为一。这样而言，假设原数二进制最高位是 i，那么在新数中它第 i 个块以内会不触发退位。不难发现这个数满足我们化好的形式，直接做就行。

完整代码如下：

u64 ut_blkpstv(u64 x) {
  x = (x & 0x7f7f7f7f7f7f7f7full) | (x & 0x8080808080808080ull) >> 1;
  return (x | 0x8080808080808080ull) - 0x0101010101010101ull;
}
int ut_blkfst(u64 x) {
  x = x >> 7 & 0x0101010101010101ull;
  return (ut_lg64(x * 0xff + 1) >> 3) - 1;
}
int ut_bsrsi8(u8 x) {
  if (x >> 7) return 7;
  u64 y = (x | 0x80) * 0x0101010101010101ull;
  return ut_blkfst(y - 0x8040201008040201ull);
}
int ut_bsrsi64(u64 x) {
  u64 y = ut_blkpstv(x) >> 7 & 0x0101010101010101ull;
  int c = (y >> 32) ? (y >>= 32, 8) : 0;
  c += ut_blkfst(ut_blkpstv(y * 0x01010101 >> 24)) << 3;
  return c | ut_bsrsi8(x >> c);
}

以上。