《代码之美》（十）－－寻求快速的种群计数

Henry S. Warren, Jr.

种群计数（Population Count）或横列和（Sideways Sums）表示一个计算机字长所包含的“1”的总位数。有些计算机有相应的指令，但仅支持RISC或CISC的计算机并没有，它们所支持的指令为：移位 （shift）、加（add）、与（and）、取数（load）、条件分支（branch）等等。本章就解决了基于支持RISC或CISC的计算机上的种 群计数问题。

基本方法

最容易想到的方法就是下面的C程序：

pop = 0;
for (i = 0; i < 32; i++) {
    if (x & 1) pop = pop + 1;
    x = x >> 1;
}

上述循环体代码会编译成约7条指令，其中2条属于条件分支（其中一条用于循环控制）。7条指令循环执行的次数为32次，但其中if里的指令平均只执行一半的时间，所以循环体总共大约执行了32x6.5=208条指令。

上面的代码可以改进。首先，在大多数计算机上，从31倒计数到0比从0累加到31更高效，因为前者可以少一条比较指令。（博主：没搞懂为什么……）

继续改进：我们并不需要对所有的位都逐个计数，当x为0的时候，循环就可以结束了，这样当高位是0时就可以省去一些迭代。此外，if条件判断也可以去掉：

pop = 0;
while (x) {
    pop = pop + (x & 1);
    x = x >> 1;
}

上述循环体代码仅包含4条或5条RISC指令（取决于x与0是否需要比较），且只有一个分支。当x以1开始时，代码最多，耗费128～160条指令。当x有很多前缀0的话，所耗费的指令数会少很多。

表达式x&(x-1)可以将x的最低位清零，所以我们可以用如下方式来改进代码：

pop = 0;
while (x) {
    pop = pop + 1;
    x = x & (x - 1);
}

和前面那段代码一样，上述循环体的代码也需要花费4条或5条指令，但它的循环次数等于x中包含的1位的数目。

如果1位数目的期望值很大，可以用表达式x = x | (x + 1)来把x最右边的0位翻转为1位，直到所有的比特位都变成1为止（x值为－1）。这时迭代次数n就是0位的数目，32-n就是1位的数目。

最后还可以用查表法，把x中的1个字节转换成该字节所包含的1位的数目：

static char table[256] = {0, 1, 1, 2, 1, 2, 2, 3, ..., 8};
pop = table[x & 0xFF] + table[(x>>8) & 0xFF] + table[(x>>16) & 0xFF] + table[x>>24];

上述代码非常简短，且在很多机器上的运行速度都相当快。在不支持索引取数（indexed loads）的初级RISC计算机上，大约是17条指令。

分治法

我们可以把字长的32位分为两部份，每部分16位，最后的结果为两个部份中1位的数目的总和。同样，对于16位的部分，也可以分成8位的两部分。以此类推，直到每个部分只包含一个比特位。

如果不能对各部分并行计算的话，这样的分治毫无意义，因为额外的一个加操作（把两部分的1位数目相加）会增加开销。（博主：之前的循环方法复杂度为 O(n)，分治法同样为O(n)）多处理器或支持SIMD指令的计算机上可以通过并行可以很好地加速分治算法（单指令多数据源（SIMD）：可以并行操作 计算机字中的多个字段，比如字节或者半字）。

下面给出一个在单处理器的RISC或CISC计算机上并行分治的算法。

给出如下一个32位的字：10111100011000110111111011111111，用分治法来计算它的1位数目。
10    11    11    00    01    10    00    11    01    11    11    10    11    11    11    11 －－ 分治后的最末端
01    10    10    00    01    01    00    10    01    10    10    01    10    10    10    10 －－ 2个比特位中1位数目的和（二进制表示）。
  0011        0010        0010        0010        0011        0011        0100        0100   －－ 4个比特位中1位数目的和（二进制表示）。
      00000101               00000100                0000110                  00001000         －－ 8个比特位中1位数目的和（二进制表示）。
              0000000000001001                               0000000000001110                 －－ 16个比特位中1位数目的和（二进制表示）。
                            00000000000000000000000000010111                                  －－ 32个比特位中1位数目的和（二进制表示）。   

这里再对上面的分治算法作补充说明：

1、第一行是32位字本身的值，根据分治法，最后会分为位数仅为1的32个部分。在此基础上，便是开始对各部份的1位数目进行求和了。

2、第二行求出每2个比特位中1位的数目。比如00中没有1位，得到的也是00；01或10中有1个1位，得到的便是01，而11有2个1位，得到的便是10。

3、以此类推，将求得的各个小部分依次合并成更大的部分，直到最后求出32位字长中1位数目的总和（为23）。

对应于第一行到第二行的转换，我们可以这样实现：
x = (x & 0x55555555) + (x & 0xAAAAAAAA) >> 1;

0x55555555对应的二进制表示为
01  01  01  01  01  01  01  01  01  01  01  01  01  01  01  01
x & 0x55555555便求出奇数位（从右往左）的1的数目。

相应的0xAAAAAAAA对应的二进制表示为
10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10
(x & 0xAAAAAAAA ) >> 1便求出了偶数位的1的数目（并把结果放在奇数位）。

之后把奇数位１的数目与偶数位１的数目相加，便得到了第二行的结果。（通过这种移位（shifting）和屏蔽（masking）的方式可以并行地把分治 的最末端的值求出。这里基于一个特性：两个位１位数目的和最多只为２（即二进制的10），只占两位，所以不会出现进位的情况。）

对于表达式(x & 0xAAAAAAAA) >> 1，我们可以改写成(x >> 1) & 0x55555555。之所以这样改写，是为了避免在寄存器中产生两个巨大的常量。（博主：不是很明白为什么）

对于其它行的转换，也可以进行类似的处理，所以前面的转换过程用Ｃ语言描述为：

x = (x & 0x55555555) + ((x >> 1) & 0x55555555);
x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
x = (x & 0x0F0F0F0F) + ((x >> 4) & 0x0F0F0F0F);
x = (x & 0x00FF00FF) + ((x >> 8) & 0x00FF00FF);
x = (x & 0x0000FFFF) + ((x >> 16) & 0x0000FFFF);

以上代码仍然可以继续改进：
１、第一行可以改写为x - ((x>>1) & 0x55555555)（博主：证明如下：仅考虑两个比特位pq，我们只需证明pq & 01 + (pq >> 1) & 01 == pq - (pq >>1) & 01，就可以证明上面的改写是正确的。pq & 01的值为0q，(pq >> 1) & 01的值就为0p，所以上面的等式可以改写为0q + 0p == pq - 0p，等价于0p+0p == pq-0q == p0。显然，p为１或０时等式都成立，证明完毕。），这样可以少用一条指令（少了一条与操作的指令）。

２、第二行无需改进。

３、第三行里，两个4字节的部分合并时，8字节里1位的数目的和最多为0100（二进制），不会超过1111（二进制）即0x0F，所以(x & 0x0F0F0F0F) + ((x>>4) & 0x0F0F0F0F)不会因进位超出掩码0x0F0F0F0F的范围，所以这行可以改写为x = (x + (x >> 4)) & 0x0F0F0F0F。

４、同理，可以改写为(x+(x>>8)) & 0x00FF00FF。

５、最后一行x >> 16后，前面16位就都是０了，所以再进行0x0000FFFF的与操作就是x >> 16本身。所以本行改写为x & 0x0000FFFF + (x>>16)。

６、把第四行中的0x00FF00FF操作推迟到第五行，这样第四行就简化为x = x + (x >> 8)，而第五行变成了((x & 0x00FF00FF) & 0x0000FFFF) + ((x & 0x00FF00FF) >> 16) => x & 0x000000FF + ((x>>16) & 0x000000FF)，同第三行与第四行，可以提取出来成为(x + (x >> 16)) & 0x000000FF。因为32位字长里，１位数最多为不超过10000（二进制），所以这个值也可以为(x + (x>>16)) & 0x0000003F。

改进后的代码如下：

int pop(unsigned x) {
    x = x - ((x >> 1) & 0x55555555);
    x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
    x = (x + (x >> 4) & 0x0F0F0F0F;
    x = x + (x >> 8);
    x = x + (x >> 16);
    return x & 0x0000003F;
}

上面的代码只执行了21条指令，并且没有分支（branches）和内存引用（memory references）。不过高效的代价就是把它扩展到６４位机器时，并不是直接明了的事情。仅管如此，这种节约指令的收益还是很大的。

种群计数的公式如下：
pop(x) = x - floor(x/2) - floor(x/4) - floor(x/8) - ... - floor(x/2^31)

对于上面这个公式，可以这样理解：floor(x/2)是x进行右移1位的结果，floor(x/2^n)是x右移n位的结果。上式相当于循环地将x减去 x右移的结果，总共31次。在x对应的二进制序列中的第m位（从低到高），如果它是"1"，它对应的十进制值为2^m，第一次右移，对应该位被减掉的值为 2^(m-1)，...，第m-1次右移后，对应该位被减掉的值为2，之后该位就不再作为减数中的一部分了。如此，第m位的“１”最终结果为 2^m-(2^(m-1) + 2^(m-2) + ... + 2 + 1) = 1。也就是说，每多一个值为“１”的位，最终结果就能加１。这样产生的最后结果就是我们要求的种群计数。

其实上面的代码的第一行x = x - ((x>>1) & 0x5555555);便是并行地对16个2比特位使用公式的前两项：2个比特位b的1位数目和为b-floor(b/2) = b - (b>>1)。

其他方法

另一种类似于分治法的方法是求出相邻的3位字段中1位的数目，再把相邻3位字段合并成6位字段，之后对字x模63便能得到所有6位字段的位数总和。这个算法原本是为36位字长的计算机设计的，但是很容易适用于32位字。算法的C代码如下（长串的常量是8进制的）：

int pop(unsigned x) {
   unsigned n;

   n = (x >> 1) & 033333333333; //对每个
   x = x - n;                   //3位字段
   n = (n >> 1) & 033333333333; //进行1位计数。
   x = x - n;
   x = (x + (x >> 3)) & 030707070707; //6位字段的1位数目和。
   return x%63; // 将所有6位字段和相加。
}

对上述代码进行补充说明：
1、前面4行并行地求出了相邻3位的1位数目，原理和分治法一样，使用了前面计数公式的前三项：x-floor(x/2)-floor(x/4)，所以代码前4行做了两个移位操作，以及两个减法操作；

2、第5行将相邻两个3位字段合并，原理与分治法相同；

3、最后的模63的由来：可以把一个字分为每组6位的各组，假设（从低位到高位）第i组（最低位的那组为第0组）中6个比特位表示的值为bi，那么（32 位）字x的值为b0 + b1*64 + b2*64^2 + .. + bi*64^i +... + b6*64^6。因为(a+b) % c = a%c + b%c - k*c，而b0+b1+b2+..+b6的值不会超过32<63，所以x%63的值可以表示为b0%63 + b1*64%63 + b2*64^2%63 +... + bi*64^i%63 + ... + b6*64^6%63。又因为a*(b+1) % b = a，从而推出a*(b+1)^i % b = a，所以x%63 = b0 + b1 + b2 + ... + b6，即字x中所有1位的数目总和。

这种算法在DEC PDP-10计算机上只需要10条指令，因为这种计算机支持这种指令：通过直接引用内存全字得到的值，将其余为余数计算的第二个操作数使用。在一台 RISC计算机上，大约需要15条指令，并且假设这台计算机提供无符号模指令。这条指令可能运行得不是很快，因为除法总是较慢。另外，这种算法最大适用于 62位字长，较难扩展到64位字长。

另一种令人惊奇的方法是把x逐次循环左移１个位置，总共３１次，然后求３２项的和：
pop(x) = -∑<0→31>(x rot<< i) （rot<<表示循环左移，即左移时最高位的值填到最低位上来）

举个4比特位1001的例子：循环左移三次，得到的值分别为0011, 0110, 1100。这样各项和为1001+0011+0110+1100 = 1110即-2的补码，也即表示原来有2个１位。其实这种方式很好理解，对于每个１位，都为最后的和贡献出1111的值，即-1。所以字长中有多少个１ 位，最后的和就有多少个-1。

这种算法只是比较新奇而已，在大多数计算机上并没有什么实用价值，因为循环要执行３１次，从而需要６３条指令外加循环控制的开销。

两个种群的和与差

计算pop(x)+pop(y)的值，有一种方法可以节省部分时间：首先用分治法中的前两行分别作用于x和y（分别把x，y中的相邻１位合并成２位字段，再合并为４位字段），之后把x与y的值相加，再用分治法的后几行把它们的和从4位字段合并成32位字。

对于pop(x) - pop(y)，我们可以把它转换为pop(x) - (32 - pop(~y)) = pop(x) + pop(~y) - 32。~y是y的反码，pop(~y)是y中0的个数，所以32-pop(~y)是y中１的个数。这样就可以和pop(x)+pop(y)使用一样的方法 来计算结果。下面是pop(x) - pop(y)的Ｃ代码：

int popDiff(unsigned x, unsigned y) {
    x = x - ((x >> 1) & 0x55555555);
    x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
    y = ~y;
    y = y - ((y >> 1) & 0x55555555);
    y = (y & 0x33333333) + ((y >> 2) & 0x33333333);
    x = x + y;
    x = (x & 0x0F0F0F0F) + ((x >> 4) & 0x0F0F0F0F);
    x = x + (x >> 8);
    x = x + (x >> 16);
    return (x & 0x00000007F) - 32;
}

两个字的种群计数比较

有时我们只想知道两个字中哪一个的种群计数更大，而不考虑实际的计数值是多少。基本思想是对每个字中的单个比特位轮流清零，直到其中某个字的比特位全为０，那剩下的非零种群计数值的字就是较大者：

int popComp(unsigned xp, unsigned yp) {
    unsigned x, y;
    x = xp & ~yp;  // 清除同
    y = yp & ~xp;  // 是１的比特位
    while (1) {
        if (x == 0) return y | -y;
        if (y == 0) return 1;
        x = x & (x - 1);  // 每个字
        y = y & (y - 1);  // 清除一个１位
    }
}

xp & ~yp和yp & ~xp可以清掉xp和yp中同时为１的比特位，这样可以减少后面的循环次数。这种清理后，x和y中１位的数目和不会超３２，而其中较少１位数目的字最多只 有１６位，所以后面的循环最多为１６次。最坏情况下，初级RISC计算机需要耗费119条指令（16 x 7 + 7）。y | -y在y为0时的值为０，在y不为０时的值为-1。

数组中的１位种群计数

在没有种群计数指令的情况下，计算一个全字数组的１位数目最简单的方法是：求出每个字的１位数目，再把结果相加。这种方法被称为“幼稚的方法”。

另一种方法，是先对数组中的三字组分别进行１位到２位，２位到４位的合并，之后再把这三个字相加。因为３个４位字段的和不会超过１２，所以不会溢出到左边的相邻字段。这种思想同样适用于８位和１６位字段。这种方法大概比幼稚方法减少20%的指令。

一种更优的方法由Robert Harley和David Seal在1996年左右发明的。它是一种基于进位保留器（CSA）或者3-2编码器（3:2 compressor）的电路。简单来说，CSA就是一组独立的全加法器（由３个１位输入（用于相加）和２个１位输出（和与进位）所组成的电路）。

下面是使用布尔代数符号（并置表示与，＋表示或，⊕表示异或）所表示的每个全加法器的逻辑：
h ← ab + ac + bc = ab + (a + b)c = ab + (a ⊕ b)c
l ← (a ⊕ b) ⊕ c
其中a，b，c都是１位的输入，l是低位输出（和），h是高位输出（进位）。第一行a＋b之所以可以改成a⊕b是经过论证的：a和b都为１时，ab值为１，所以两种形式下值都为１；当a和b都为０时，两种形式的值都为０；当a和b中只有一个１时，ab值为０，a＋b和a⊕b的值都为１。

一种利用CSA的方法是将三个字减为两个字，下面是１６位字的处理过程：
a = 0110 1001 1110 0101 -- 9
b = 1000 1000 0000 0111 -- 6
c = 1100 1010 0011 0101 -- 8
----------------------------------------------
l = 0010 1011 1001 0111 -- 9
h = 1100 1000 0111 0101 -- 7*2 = 14
高位的值是进位，所以最后算１位数总和时，要乘以２。

使用上述的方法，可以把数组中的三个一组地处理，把每组中的三个字减为两个字，再对这两个字应用种群计数 操作（高位中的总数最后要乘以２）。得到的结果在循环体中累计。在典型的RISC计算机上，CSA步骤的指令数nc = 5，为某个字进行种群计数的指令数np = 15。如果忽略数组取数和循环控制，上述方法处理数组中的每个字大约花费(nc + 2*np + 2) / 3约为12.33条指令（式中的+2表示循环体中的２次加法）。

另一种运用CSA的操作方法如下：

#define CSA(h, l, a, b, c) \
    {unsigned u = a ^ b; unsigned v = c; \
        h = (a & b) | (u & v); l = u ^ v; }
int popArray(unsigned A[], int n) {
    int tot, i;
    unsigned ones, twos;
    tot = 0; // 初始化
    ones = 0;
    for (i = 0; i <= n - 2; i = i + 2) {
        CSA(twos, ones, ones, A[i], A[i+1])
        tot = tot + pop(twos);
    }
    tot = 2 * tot + pop(ones);
    if (n & 1)                  // 如果还剩最后一个字，
        tot = tot + pop(A[i]);  // 则把它的种群计数也加上。
    return tot;
}

上述代码持续地把相邻两个字的和的低位累加到ones上，而把和的高位twos用计数器tot记录下来。这种方法平均每字耗费(nc + np + 1) / 2 = 10.5条指令。
上面的CSA宏展开后的代码如下：

u = ones ^ A[i];
v = A[i+1];
tows = (ones & A[i]) | (u & v);
ones = u ^ v;

我们可以进一步减少指令数。这种方法可以通过下方的电路图来解释：

我们之前只是用ones来记录低位的值，把twos的值累加到计数器tot上，从而一次数理两个字。现在，我们增加几个CSA单元，这样可以一次处理８个 字。整个数组处理完毕后，最后的结果为8 x pop(eights) + 4 x pop(fours) + 2 x pop(twos) + pop(ones)。忽略数组取数和循环控制，每个字的循环体的执行时间为(7*nc+np+1)/8 = 6.375条指令。下面是对应的C代码：

int popArray(unsigned A[], int n) {
    int tot, i;
    unsigned ones, twos, twosA, towsB, fours, foursA, foursB, eights;
    tot = 0;
    fours = tows = ones = 0;
    for (i = 0; i <= n - 8; i = i + 8) {
        CSA(twosA, ones, ones, A[i], A[i+1])
        CSA(twosB, ones, ones, A[i+2], A[i+3])
        CSA(foursA, twos, twos, twosA, twosB)
        CSA(twosA, ones, ones, A[i+4], A[i+5])
        CSA(twosB, ones, ones, A[i+6], A[i+7])
        CSA(foursB, twos, twos, twosA, twosB)
        CSA(eights, fours, fours, foursA, foursB)
        tot = tot + pop(eights);
    }
    tot = 8 * tot + 4 * pop(fours) + 2*pop(twos) + pop(ones);
    for (i = i; i < n; i++)
        tot = tot + pop(A[i]);
    return tot;
}

下面是各类方法下数组种群计数的平均每次花费的指令数：

程序	取数和循环控制除外的指令		循环中的所有指令（编译器输出）
	公式	令nc = 5，np = 15
幼稚方法	np+1	16	21
分组大小为2	(nc + np + 1) / 2	10.5	14
分组大小为4	(3*nc + np + 1) / 4	7.75	10
分组大小为8	(7*nc + np + 1) / 8	6.38	8
分组大小为16	(15*nc + np + 1) / 16	5.69	7
分组大小为32	(31*nc + np + 1) / 32	5.34	6.5
分组大小为2^n	nc + (np - nc + 1) / 2^n	5 + 11 / 2^n

应用

种群计数指令的用途很广泛，比如：
１、计算位串（bit strings）所表示的集合大小；

２、计算两个位向量之前的“Hamming距离”，这一应用来自于错误校验码理论。Hamming距离指的是向量间对应位取不同值的位数，即：dist(x, y) = pop(x⊕ y)

３、计算字中后缀0的数目，使用关系式：ntz(x) = pop(┐x & (x - 1)) = 32 - pop(x | -x)可以得到。

这里x作为无符号数。x-1可以把离最低位最近的“１”位变为“011111...”序列，而┐x与x-1的交集正是这些最后的“11111...”序 列，也即x中后缀0的数目。另一方面，根据补码的特性，-x = ┐x + 1 ，这个+1操作将使得从低位开始的连续1序列（对应着x中的后缀0序列）变为0序列，x | -x的值除了+1操作产生的低位连续的0序列外，其余都为１，所以pop(x | -x)就是除了后缀0序列的位数，32 - pop(x | -x)便是x中后缀0的个数。

ntz(x)函数用途同样非常广泛，比如中断处理。早期计算机一旦发生中断，就会在一个特殊的寄存器中存储一个表示“中断原因”的位。优先级越高的中断位离最低位越近。可以通过计算后缀０个数的方式找到最高优先级的中断并处理。

４、快速访问使用特定压缩方式表示的中等稀疏数组。

稀疏数组可以被压缩，在压缩的数组里仅存储已定义的或者非零的数组元素。为了维护原数组的信息，需要额外一个位串。如果该位串的第i位是１，就说明原数组 A[i]有定义。这个位串通常很长，所以它会被分成一组32位的字，这样就会有一个bits数组，数组中每个元素是一个３２位的字。另外维护一个 bitsum的数组，其中元素bitsum[j]表示bits中前j个字里的1位数目的总和，它用来加速索引。下表描述了一个将95个元素压缩成6个元素 的结果：

bits	bitsum	data
0x00000005	0	A[0]
0x00018001	2	A[2]
0x80000000	5	A[32]
		A[47]
		A[48]
		A[95]

在进行数组索引时，需要将“逻辑”下标i转换为“物理”索引sparse_i，比如将A[48]转换为data[5]。48是逻辑下标，5是物理下标。下面是转换的代码：

j = i >> 5; // j = i/32
k = i & 31; // k = rem(i, 32),取余（取模）
mask = 1 << k; // "1"位于第k位。
if ((bits[j] & mask) == 0) goto no_such_element;
mask = mask - 1; // 第k位右边全为1。
sparse_i = bitsum[j] + pop(bits[j] & mask);

j*32 + k = i，通过j，k找到位串中的对应于A[i]的位。最后通过该位之前的1的个数，来确定它真正对应的物理存储的位置。

５、用来生成符合二项式分布的随机整数。

６、计算机界传言，种群计数对于美国国家安全局（National Security Agency）而言非常重要。除了NSA自己的人，没有人知道他们用种群计数做什么，但很有可能是关于加密的工作或者海量材料的搜索。

：《高效程序的奥秘》（Hacker's Delight）的作者，他在IBM有过40年的工作经验。