位运算与组合搜索（二）

People who play with bits should expect to get bitten.

-- Jurg Nievergelt

I failed math twice, never fully grasping probability theory.

I mean, first off, who cares if you pick a black ball or a white ball out of the bag?

And second if you’re bent over about the color, don’t leave it to chance.

Look in the damn bag and pick the color you want.

-- Stephanie Plum

这篇文章接着讲怎样高效地遍历所有的组合。同样，假定全集的大小不大于机器字长，计算模型为 word-RAM，即诸如 +, –, *, /, %, &, |, >>, << 等此类操作皆可以在 O(1) 时间内完成。当然，一般 / 和 % 操作明显要慢上一些，因此我们总是希望能够尽量避免使用两个操作。在上一篇 blog 中的子集遍历比较简单，基本上只用到了 +, – , & 三种操作。而组合遍历相对要复杂得多，一些让人不舒服的操作总是难以避免。下面将要介绍两种完全不同的组合遍历算法，其中一种用到了 / 操作，而另一种则使用了三目运算符。尽管不算十分完美，也应该是足够高效啦。

---

2 遍历所有组合

2.1 colex & reverse colex

说到各种位运算技巧，早年从 MIT 流传出来的一份技术报告 HAKMEM 可谓是一本黑暗圣经。在 HAKMEM 的第175条中记录着一个非常巧妙而实用的技巧，被称为 Gosper’s hack，它仅仅使用几个非常简单的算术运算和位运算，即可得到与当前所输入的整数含有相同数目的1的下一个整数：

s = x & (-x);

r = s + x;

n = r | (((x ^ r) >> 2) / s);

在上面这段代码中 x 是输入，n 是输出，为大于 x 且与 x 含1个数相同的最小整数。比如若输入 x = 0b0101, 那么将输出 n = 0b0110。使用这一技巧使得我们可以非常容易地生成所有组合，代码如下：（这是一个成员函数，完整的代码可在位运算与组合搜索（一）所附的压缩包中找到）：

bool next(unsigned long &x) const
{
    if (x == last_) return false;
    unsigned long r, s;
    s = x & (-(long)x);
    r = s + x;
    x = r | (((x ^ r) >> 2) / s);
    return true;
}

上面代码中的 last_ 表示的是最后一个组合。这里遍历组合的序为 colex，最小的组合是所有1都在低位，而最大的组合（即 last_） 是当所有1都在高位。比如若全集为 {a, b, c, d, e}，我们用以上代码遍历其所有大小为2的子集，顺序将如下表所示：

序号	值	位串	子集
1	0b00011	11000	{a, b}
2	0b00101	10100	{a, c}
3	0b00110	01100	{b, c}
4	0b01001	10010	{a, d}
5	0b01010	01010	{b, d}
6	0b01100	00110	{c, d}
7	0b10001	10001	{a, e}
8	0b10010	01001	{b, e}
9	0b10100	00101	{c, e}
10	0b11000	00011	{d, e}

现在稍微来解释 Gosper’ hack 是怎样工作的：

第一条语句：s = x & (-x)， 用于标识出 x 最低位的1（设最低的1右边有 c 个0）。 e.g. 0b10110 –> 0b00010

第二条语句：r = s + x， 将 x 右端的连续一段1清零（红色标识的部分，设这一段有 k 个1），并将前一位设为1。 e.g. 0b10110 –> 0b11000

第三条语句：n = r | (((x ^ r) >> 2) / s)， 这里先用 x 异或 r 得到 k + 1 + c 个连续的1。然后右移 2 位，再除于 s （相当于右移 c 位），得到 k – 1 位连续的1，最后添加到 r 的最右边，打完收工。e.g. 0b11000 | 0b00001 = 0b11001

由于该 hack 中的除法实际上只是用来移位的，因此可以想办法绕过去 （如果你实在看不顺眼那个除号的话）。比如可以使用 bsr 指令计算出 c ，然后直接移位即可。但经过我的测试，发现还是直接除法来得比较快。

// Find last bit set
static inline unsigned long __fls(unsigned long x)
{
    __asm bsr eax, x;
}

现在如果想要反向生成所有的组合那又该如何呢，其实很简单，因为 colex 具有一种某种意义上的对称性：某个组合的前一个组合等于这个组合的补集的下一个组合的补集。如果我们想要得到组合 x 按照 colex 的上一个组合，只需生成 ~x 的下一个组合，再取反即可：

bool prev(unsigned long &x) const
{
    if (x == first_) return false;
    x = ~x;
    next(x);
    x = ~x;
    return true;
}

---

2.2 cool-lex & reverse cool-lex

cool-lex，顾名思义，就是非常 cool 的 lex。cool-lex 是由 Frank Ruskey 和 Aaron Williams 发明的，如果想要详细的了解 cool-lex 的性质，可以看一下参考文献6。另外在这里还有一段 cool-lex 的音乐，感兴趣的可以试听一下。虽然它不怎么好听，也显然不可能给你带来关于 cool-lex 的任何洞见。下面我只简单介绍一下怎样按照 cool-lex 或者反向 cool-lex 进行组合遍历。

cool-lex 的生成算法是基于后缀旋转的（如果是针对位串表示则是前缀旋转，但下面我们都是针对二进制整数表示，也就是低位在右边）：找到最短的以010或者110开始的后缀（如果不存在则选定全部位），然后向左旋转1位。比如组合0b01101,  首先找出最短的以010或者110开始的后缀（用红色表示）：0b01101，然后将这个后缀向左旋转1位（即循环左移1位）即得到下一个组合：0b01011。

如何借助于位运算高效的完成后缀旋转呢，Donald 在 TAoCP 中7.2.1.3节习题55的答案中给出了一个 MMIX 实现。下面的代码是我写的一个C++版：

bool next(unsigned long &x) const
{
    if (x == last_) return false;
 
    unsigned long r, s;
    r = x & (x + 1);
    s = r ^ (r - 1);
    r = ((s + 1) & x) ? s : 0;
    x = x + (x & s) - r;
    return true;
}

上面代码中的 last_ 当然也是指最后一个组合。cool-lex 中的第一个组合也是所有1在低位，即类似于这样：0b0…01…1。最后一个组合是1个1在最高位，而其余的1在低位，即形如 0b10…01…1。这段代码到底是怎么起作用的？你猜！我就不分析了，不过我等下会详细解释生成 reverse cool-lex 的代码。下表是 cool-lex 序的一个例子（同样，全集为 {a, b, c, d, e}，子集大小为 2）：

序号	值	位串	子集
1	0b00011	11000	{a, b}
2	0b00110	01100	{b, c}
3	0b00101	10100	{a, c}
4	0b01010	01010	{b, d}
5	0b01100	00110	{c, d}
6	0b01001	10010	{a, d}
7	0b10010	01001	{b, e}
8	0b10100	00101	{c, e}
9	0b11000	00011	{d, e}
10	0b10001	10001	{a, e}

现在来讲怎样反向遍历 cool-lex。reverse cool-lex 被提到的不多，网上以及各种文献上也并没有生成 reverse cool-lex 的代码，因此我只好自己写了一个。想要得到高效的 cool-lex 反向遍历代码，首先需要一个简单的生成规则。这个规则其实根据正向 cool-lex 的规则可以很容易地yy出来：找到最短的以100或者101开始的后缀（如果不存在则选定全部位），然后向右旋转1位。（后来我向 Frank 请教了一下，他说这个规则的确是正确的，另外还告诉我 Aaron 的另一篇文章 “loopless generation of multiset permutations by prefix shifts” 对 reverse cool-lex 作了介绍。）

规则有了，还剩下最后一个问题，那就是怎样借助于位运算高效的实现这个规则。下面是我的实现：

bool prev(unsigned long &x) const
{
    if (x == first_) return false;
    unsigned long r, s, v;
    v = x | 1;
    r = v & (v + 1);
    s = r ^ (r - 1);
    v = s & x;
    r = (v & 1) ? s - (s >> 1) : 0;
    x = x & (~s) | r | (v >> 1);
    return true;
}

上面的代码中，基本上都是非常基础的位运算技巧，如果你对此并不熟悉，不妨看一下参考文献1或3。首先，我们需要找到最短的以 100 或者 101 开始的后缀，这将通过下面四条语句来完成：

第一条语句：v = x | 1，将最低位置1。e.g. 0b01010 –> 0b01011

第二条语句：r = v & (v + 1)，清除右边连续的1。 e.g. 0b01011 –> 0b01000

第三条语句：s = r ^ (r – 1)，标记最低位的1以及其后的0。e.g. 0b01000 –> 0b01111

第四条语句：v = s & x，得到后缀。e.g. 0b01111 & 0b01010 –> 0b01010

至此，满足条件的后缀已经找出来了，下一步的工作就是将它右旋一位：

第五条语句：r = (v & 1) ? s - (s >> 1) : 0， 得到旋转后的后缀的最高位。 e.g. 0b01111 - 0b00111 –> 0b01000

第六条语句：x = x & (~s) | r | (v >> 1)，将后缀右移一位，与最高位相或，再与其余不相干的位合并，即得到最终结果。 e.g. 0b00101

在第五条语句用到了三目运算符，这里其实也可以借助 bsr 指令绕过去。我并没有比较哪种更快一些。

完。（做人要厚道，转载请注明出处：http://www.cnblogs.com/atyuwen/）

---

3 参考文献

1. Henry S. W., Hacker’s Delight.
2. Jörg A., Matters Computational.
3. Donald E. K., The Art of Computer Programming: Bitwise Tricks and Techniques. Volume 4, Pre-Fascicle 1A.
4. Donald E. K., The Art of Computer Programming: Generating all Combinations and Partitions. Volume 4, Fascicle 3.
5. Beeler M., Gosper R. W., and Schroeppel R., HAKMEM
6. Frank R., Aaron W., The Coolest Way To Generate Combinations.

P.S. 对于我这种从小就怕写作文的人来说，写篇稍正式一点的技术文章实在是太辛苦了，因此关于位运算与组合搜索就先写到这里，虽然有很多想说的还没有谈到。以后有心情了再来讨论怎样高效实现（一）中提到的映射操作及其逆操作。