写者按:这篇文章原计划写成Bit Twiddling Hacks的一篇翻译,但因原文除了bit位操作技巧还花费不少笔墨记录了发现者与修订者,本文不含上述内容,只注重bit相关的操作技巧,对这些代码作者感兴趣的请移步观看。

前言

在统计总操作数时,任何一个c语句认为是一个操作。间接赋值语句因其不用写入内存,认为其占0个操作数。当然,操作总数只是近似表示真实的计算机指令数和操作时间。本文假定所有的操作花费的时间相同,虽然实际上是不可能的,但在CPUs速度快速增长的今天,其中的差距已经非常小了。一段样例代码的执行时间还和系统的缓存大小、内存带宽、指令集等有关。实际测试是考较一段代码效率的最终手段,以下代码执行结果以你机子实际测试结果为准。

求整型的符号

int v;      // 输入,求该整型的符号
int sign;   // 整型符号(-1表示负数,0表示正数或0) 

// CHAR_BIT 指每字节bit数 (一般为 8).
sign = -(v < 0);  // if v < 0 then -1, else 0. 
// 或者,避免分支判断的方法如下:
sign = -(int)((unsigned int)((int)v) >> (sizeof(int) * CHAR_BIT - 1));
// 或者,更少的操作数但是可移植行较差的方法:
sign = v >> (sizeof(int) * CHAR_BIT - 1); 

上述代码的最后一个表达式用sign = v >> 31计算一个32位整型的符号。这步操作会比直观的sign = -(v < 0)更快。该方法之所以行的通是因为,符号整型右移时,最左的符号位会向右填充。负整型最左符号位为1, 0和正整型符号为则为0;所有位为1表示-1(补码表示法)。不幸的是,该方法只适用与特定的架构。

或者,你可能希望结果表示成-1或+1,那么可以用:

sign = +1 | (v >> (sizeof(int) * CHAR_BIT - 1));  // if v < 0 then -1, else +1

又或者,你希望结果表示成-1、0或+1,那么用:

sign = (v != 0) | -(int)((unsigned int)((int)v) >> (sizeof(int) * CHAR_BIT - 1));
// 或者,更快但可移植性较差的版本:
sign = (v != 0) | (v >> (sizeof(int) * CHAR_BIT - 1));  // -1, 0, or +1
// 或者,更简短直观且可移植行较好的版本:
sign = (v > 0) - (v < 0); // -1, 0, or +1

如果你想要知道某数是否为非负数,用+1表示非负数,0表示负数,可以用:

sign = 1 ^ ((unsigned int)v >> (sizeof(int) * CHAR_BIT - 1)); // if v < 0 then 0, else 1

警告: 1989 ANSI C specification 不再定义符号右移操作的结果,在某些系统上上述涉及右移操作的方法可能无效。

判断两个整型是否符号相反

int x, y;               // 输入,判断这两个数符号是否相反

bool f = ((x ^ y) < 0); // 如果x,y符号相反则为true

求整型的绝对值(abs, 无分支判断)

int v;           // 输入的整型
unsigned int r;  // 所求绝对值 
int const mask = v >> sizeof(int) * CHAR_BIT - 1;

r = (v + mask) ^ mask;
// 相同效果的变式:
r = (v ^ mask) - mask;

某些CPU没有求整型绝对值的指令(或者编译器没使用它)。在某些机器上分支判断操作可能是比较昂贵的,上述表达式比简单方法r = (v < 0) ? -(unsigned)v : v要快,即使它们的操作数是相同的。

求两个整型的最大最小值(无分支判断)

int x;  // 求x,y 的最小值
int y;   
int r;  // 结果

r = y ^ ((x ^ y) & -(x < y)); // min(x, y)

在某些机器上分支判断操作可能是比较昂贵的,上式比简单方法r = (x < y) ? x : y要快,即使上式多了2个操作。(当然,一般来说最好使用简单方法,易于理解。)上式的原理是:如果x<y,那么-(x-y)全是1,所以r = y ^ (x ^ y) & ~0 = y ^ x ^ y = x。此外,如果x>=y,那么-(x<y)全是0, 所以r = y ^ ((x ^ y) & 0) = y。在某些机子上,x<y可能需要一个分支判断指令,上式将比简单方法差。

求最大值:

r = x ^ ((x ^ y) & -(x < y)); // max(x, y)

判断一个整型是否是2的幂

unsigned int v; // 输入
bool f;         // 若v是2的冪,f为true

f = (v & (v - 1)) == 0;

上式会错误地将0认为是2的幂,为了修正这点,使用:

f = v && !(v & (v - 1));

定长符号扩展

对于内建类型(chars ints)符号扩展是自动完成的,但假设你有一个有符号使用b位补码表示的数字x,想将之强转成一个int。如果x是正数,简单拷贝即可,但如果是负数,符号位必须调整到最高位。举个例子,我们用4位表示-3是1101。如果用8位表示-3则是11111101。可以看到左边多出的4位使用了符号位填充,即符号扩展。在C中,将定长类型进行符号扩展是很简单的。如下式子将一个5bits扩展成一个整型:

int x; // 输入,该int只使用最低的5 bits
int r; // 符号扩展的结果
struct {signed int x:5;} s;
r = s.x = x;

以下C++模板函数使用相同的方式将B bits的转化成类型T。

template 
inline T signextend(const T x)
{
  struct {T x:B;} s;
  return s.x = x;
}

int r = signextend(x);  // 将5 bits 符号扩展为整型r

条件置位

bool f;         // flag
unsigned int m; // 掩码
unsigned int w; // 待置位的字:  if (f) w |= m; else w &= ~m; 

w ^= (-f ^ w) & m;

// 或者,对于超标量CPUs
w = (w & ~m) | (-f & m);

对于某些架构,无需分支判断的版本能弥补多出的操作数。在非正式速度测试中,在AMD Athlon™ XP 2100+上上式快5-10%。第二个式子在Intel Core 2 Duo 超标量版本中比第一个式子要快16%。

条件取负

如果你想仅当某flag为false是将某数取反,可以使用以下方法避免分支判断:

bool fDontNegate; 
int v;             // 输入,如果fDontNegate为false,v取反
int r;             // r = fDontNegate ? v : -v;

r = (fDontNegate ^ (fDontNegate - 1)) * v;

如果你想仅当flag为true时取反,那么使用:

bool fNegate;  
int v;         // 输入,如果fDontNegate为true,v取反
int r;         // result = fNegate ? -v : v;

r = (v ^ -fNegate) + fNegate;

根据掩码和并位

unsigned int a;    
unsigned int b;    
unsigned int mask; // 掩码,如果某位为0则从a取值,否则从b取
unsigned int r;    // r = (a & ~mask) | (b & mask)

r = a ^ ((a ^ b) & mask); 

上式比注释中的简单方法少了一个操作,但如果mask是常量的话,上式将不具优势。

统计bits数(简单方法)

unsigned int v; // 统计v中1的个数
unsigned int c; // 结果

for (c = 0; v; v >>= 1)
{
  c += v & 1;
}

简单法每个bit循环一次,对于最高位为1的32-bit的字,将循环32次。

统计bits数(查表法)

static const unsigned char BitsSetTable256[256] = 
{
#   define B2(n) n,     n+1,     n+1,     n+2
#   define B4(n) B2(n), B2(n+1), B2(n+1), B2(n+2)
#   define B6(n) B4(n), B4(n+1), B4(n+1), B4(n+2)
    B6(0), B6(1), B6(1), B6(2)
};

unsigned int v; // 统计v中1的个数
unsigned int c; // 结果

// 方法 1:
c = BitsSetTable256[v & 0xff] + 
    BitsSetTable256[(v >> 8) & 0xff] + 
    BitsSetTable256[(v >> 16) & 0xff] + 
    BitsSetTable256[v >> 24]; 

// 方法 2:
unsigned char * p = (unsigned char *) &v;
c = BitsSetTable256[p[0]] + 
    BitsSetTable256[p[1]] + 
    BitsSetTable256[p[2]] +	
    BitsSetTable256[p[3]];


// 或者用一下语句初始化查询表:
BitsSetTable256[0] = 0;
for (int i = 0; i < 256; i++)
{
  BitsSetTable256[i] = (i & 1) + BitsSetTable256[i / 2];
}

统计bits数(Brian Kernighan法)

unsigned int v; // 统计v中1的个数
unsigned int c; // 结果
for (c = 0; v; c++)
{
  v &= v - 1; // 清除最低位的1
}

Brian Kernighan法循环次数于bits数相同,对于仅有最高位为1的32-bit字也只需循环一次。该方法最早出现1988年《C Programming Language》第二版的习题2-9。(作者是Brian W. Kernighan and Dennis M. Ritchie)

统计bits数(并行法)

unsigned int v; // 统计v中1的个数
unsigned int c; // 结果
static const int S[] = {1, 2, 4, 8, 16}; // 神奇的2进制数
static const int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF, 0x0000FFFF};

c = v - ((v >> 1) & B[0]);
c = ((c >> S[1]) & B[1]) + (c & B[1]);
c = ((c >> S[2]) + c) & B[2];
c = ((c >> S[3]) + c) & B[3];
c = ((c >> S[4]) + c) & B[4];

数组B中的数字用2进制表示是:

B[0] = 0x55555555 = 01010101 01010101 01010101 01010101
B[1] = 0x33333333 = 00110011 00110011 00110011 00110011
B[2] = 0x0F0F0F0F = 00001111 00001111 00001111 00001111
B[3] = 0x00FF00FF = 00000000 11111111 00000000 11111111
B[4] = 0x0000FFFF = 00000000 00000000 11111111 11111111

对于更大的整型,我们只需扩展B和S数组。对于k bits,S数组与B数组长度为ceil(lg(k))。(此处lg是以2为底的对数函数,下同,写者注)

统计32-bit字bits数最好的方法是下式:

v = v - ((v >> 1) & 0x55555555);                    // 00->00, 01->01, 10->01, 11->10
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);     // 临时变量
c = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24; // 结果

根据类型T可以将上述方法参数化(对128及以内长度类型都适用):

v = v - ((v >> 1) & (T)~(T)0/3);                           
v = (v & (T)~(T)0/15*3) + ((v >> 2) & (T)~(T)0/15*3);      
v = (v + (v >> 4)) & (T)~(T)0/255*15;                      
c = (T)(v * ((T)~(T)0/255)) >> (sizeof(T) - 1) * CHAR_BIT; // 结果

上式是统计bits数最好的方法,操作数和查表法相同,但无需额外的内存开销。见Ian Ashdown’s nice newsgroup post了解更多详情,该法也称为sideways addition。

求总bits数的奇偶性

unsigned int v;       // 输入,求其总bits数的奇偶性
bool parity = false;  // 结果,用true(1)表示奇数,false(0)表示偶数

while (v)
{
  parity = !parity;
  v = v & (v - 1);
}

上式类似Brian Kernigan法统计bits总数,循环次数与总bits数相等。

求总bits数的奇偶性(查表法)

static const bool ParityTable256[256] = 
{
#   define P2(n) n, n^1, n^1, n
#   define P4(n) P2(n), P2(n^1), P2(n^1), P2(n)
#   define P6(n) P4(n), P4(n^1), P4(n^1), P4(n)
    P6(0), P6(1), P6(1), P6(0)
};

unsigned char b;  // 输入,求其总bits数的奇偶性
bool parity = ParityTable256[b];

// 对于32-bit字:
unsigned int v;
v ^= v >> 16;
v ^= v >> 8;
bool parity = ParityTable256[v & 0xff];

// 变式:
unsigned char * p = (unsigned char *) &v;
parity = ParityTable256[p[0] ^ p[1] ^ p[2] ^ p[3]];

求总bits数的奇偶性(取余法)

unsigned char b;  // 输入,求其总bits数的奇偶性
bool parity = 
  (((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;

上式仅有4个操作,仅对bytes有效。

求总bits数的奇偶性(并行法)

unsigned int v;  // 输入,求其总bits数的奇偶性
v ^= v >> 16;
v ^= v >> 8;
v ^= v >> 4;
v &= 0xf;
return (0x6996 >> v) & 1;

上式仅约9个操作,适用于32字。对一个字节可以进行一些优化(删除2,3行)。上式中的0x6996是4-bits奇偶性查找表,移位后第0位即是奇偶性的值。

交换值(异或方式)

#define SWAP(a, b) (((a) ^= (b)), ((b) ^= (a)), ((a) ^= (b)))

这是不用额外临时变量交换值的方法,相信很多人都已经知道。值得注意的是对于地址相同的a,b,比如SWAP(a[i], a[j]) i==j,上式结果将不正确。为了修正这点,考虑使用(((a) == (b)) || (((a) ^= (b)), ((b) ^= (a)), ((a) ^= (b))))表达式作为替代,其更安全。

交换bits(异或方式)方法

unsigned int i, j; // 待交换的位置
unsigned int n;    // 交换的bits个数
unsigned int b;    // 待交换的数
unsigned int r;    // 交换后的结果

unsigned int x = ((b >> i) ^ (b >> j)) & ((1U << n) - 1);
r = b ^ ((x << i) | (x << j));

例如,如果b = 00101111,我们想交换n = 3个连续的bits,分别位于i = 1(从右边往左数),j = 5的位置,结果r = 11100011。

反转bits(简单方法)

unsigned int v;     // 待反转的数
unsigned int r = v; // 反转后的数
int s = sizeof(v) * CHAR_BIT - 1; // 最后需要左移s位

for (v >>= 1; v; v >>= 1)
{   
  r <<= 1;="" r="" |="v" &="" s--;="" }="" <<="s;" 当v的左边为有x个连续的0时,需要左移x-1位<="" code="">

反转bits(查表法)

static const unsigned char BitReverseTable256[256] = 
{
#   define R2(n)     n,     n + 2*64,     n + 1*64,     n + 3*64
#   define R4(n) R2(n), R2(n + 2*16), R2(n + 1*16), R2(n + 3*16)
#   define R6(n) R4(n), R4(n + 2*4 ), R4(n + 1*4 ), R4(n + 3*4 )
    R6(0), R6(2), R6(1), R6(3)
};

unsigned int v; // 待反转bits的输入
unsigned int c; // 结果

// 方法 1:
c = (BitReverseTable256[v & 0xff] << 24) | 
    (BitReverseTable256[(v >> 8) & 0xff] << 16) | 
    (BitReverseTable256[(v >> 16) & 0xff] << 8) |
    (BitReverseTable256[(v >> 24) & 0xff]);

// 方法 2:
unsigned char * p = (unsigned char *) &v;
unsigned char * q = (unsigned char *) &c;
q[3] = BitReverseTable256[p[0]]; 
q[2] = BitReverseTable256[p[1]]; 
q[1] = BitReverseTable256[p[2]]; 
q[0] = BitReverseTable256[p[3]];

方法1需要17个操作,方法2只需12个。(假定你的CPU存取byte无须额外操作)

反转bits(乘法)

unsigned char b; // 待反转数
 
b = (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;

下边展示了该方法是如何作用的,用a,b,c,d,e,f,g,h来表示一个字节的8个位。第一步乘法展开了b的多个副本,后一个乘法将它们整合到右数第5个字节。

                                                                                        abcd efgh
*                                         (0x80200802) 1000 0000  0010 0000  0000 1000  0000 0010 
-------------------------------------------------------------------------------------------------
                                            0abc defg  h00a bcde  fgh0 0abc  defg h00a  bcde fgh0
&                            (0x0884422110) 0000 1000  1000 0100  0100 0010  0010 0001  0001 0000 
-------------------------------------------------------------------------------------------------
                                            0000 d000  h000 0c00  0g00 00b0  00f0 000a  000e 0000
*                            (0x0101010101) 0000 0001  0000 0001  0000 0001  0000 0001  0000 0001 
-------------------------------------------------------------------------------------------------
                                            0000 d000  h000 0c00  0g00 00b0  00f0 000a  000e 0000
                                 0000 d000  h000 0c00  0g00 00b0  00f0 000a  000e 0000
                      0000 d000  h000 0c00  0g00 00b0  00f0 000a  000e 0000
           0000 d000  h000 0c00  0g00 00b0  00f0 000a  000e 0000
0000 d000  h000 0c00  0g00 00b0  00f0 000a  000e 0000
-------------------------------------------------------------------------------------------------
0000 d000  h000 dc00  hg00 dcb0  hgf0 dcba  hgfe dcba  hgfe 0cba  0gfe 00ba  00fe 000a  000e 0000
>> 32
-------------------------------------------------------------------------------------------------
                                            0000 d000  h000 dc00  hg00 dcb0  hgf0 dcba  hgfe dcba  
&                                                                                       1111 1111
-------------------------------------------------------------------------------------------------
                                                                                        hgfe dcba

求余数(不用取余操作)

const unsigned int n;          // 分子
const unsigned int s;
const unsigned int d = 1U << s; // 除数,d 的可能值:1 2 4 8 16 32 ... 
unsigned int m;                // m = n % d
m = n & (d - 1); 

大部分的程序员都知道这个技巧,从完备性角度考虑还是加上它。

求以2为底的对数(简单方法)

unsigned int v;     // 输入,求v以2为底的对数
unsigned int r = 0; // r = lg(v)

while (v >>= 1)
{
  r++;
}

求一个整型以2为底的对数等价与求最高位的1(MSB)。

求以2为底的对数(查表法)

static const char LogTable256[256] = 
{
#define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
    -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
    LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
    LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
};

unsigned int v; // 输入,求v以2为底的对数
unsigned r;     // r = lg(v)
register unsigned int t, tt; // 临时变量

if (tt = v >> 16)
{
  r = (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
}
else 
{
  r = (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
}

对于32-bit数,查表法仅需7个操作。扩展到64-bit,将要9个操作。如果使用4个查找表可以节省1个操作。

以下方法适用于均匀分布的v,如果输入的v值中1是均匀分布的,考虑使用一下方法替代:

if (tt = v >> 24) 
{
  r = 24 + LogTable256[tt];
} 
else if (tt = v >> 16) 
{
  r = 16 + LogTable256[tt];
} 
else if (tt = v >> 8) 
{
  r = 8 + LogTable256[tt];
} 
else 
{
  r = LogTable256[v];
}

可以用以下方式生成查找表:

LogTable256[0] = LogTable256[1] = 0;
for (int i = 2; i < 256; i++) 
{
  LogTable256[i] = 1 + LogTable256[i / 2];
}
LogTable256[0] = -1; // 如果你想让log(0)返回-1

求尾部连续0的个数(并行法)

unsigned int v;      // 求v尾部连续0的个数 如(b00100 -> 2)
unsigned int c = 32; // 输出
v &= -signed(v);
if (v) c--;
if (v & 0x0000FFFF) c -= 16;
if (v & 0x00FF00FF) c -= 8;
if (v & 0x0F0F0F0F) c -= 4;
if (v & 0x33333333) c -= 2;
if (v & 0x55555555) c -= 1;

我们先分离最低位的1,c初始化为最大的32,然后逐步递减c的值。对于N bit数,平均操作数约等于3 * lg(N) + 4。

求尾部连续0的个数(2分查找法)

// 注意: 如果 0 == v, 那么 c = 31.
if (v & 0x1) 
{
  // 对奇数的快速处理,假定奇数占输入的1/2
  c = 0;
}
else
{
  c = 1;
  if ((v & 0xffff) == 0) 
  {  
    v >>= 16;  
    c += 16;
  }
  if ((v & 0xff) == 0) 
  {  
    v >>= 8;  
    c += 8;
  }
  if ((v & 0xf) == 0) 
  {  
    v >>= 4;
    c += 4;
  }
  if ((v & 0x3) == 0) 
  {  
    v >>= 2;
    c += 2;
  }
  c -= v & 0x1;
}	

该方法类似与上一个方法,但是是以积累的方式计算c。它先检测v最低的16bits是否为0, 如果都为0,将v右移16位,c加上16。后续操作每次减半直到为1为止。这个方法比上个方法大约要快33%,因为有些分支不会被执行。

向上取2的幂(浮点法)

unsigned int const v; // 求大于或等于v的2幂数的最小值
unsigned int r;       // 结果 ( v=3 -> r=4; v=8 -> r=8)

if (v > 1) 
{
  float f = (float)v;
  unsigned int const t = 1U << ((*(unsigned int *)&f >> 23) - 0x7f);
  r = t << (t < v);
}
else 
{
  r = 1;
}

上述方法需要8个操作,适用与所有v <= (1<<31)

向上取2的幂

unsigned int v; // 求大于或等于v的2幂数的最小值

v--;
v |= v >> 1;
v |= v >> 2;
v |= v >> 4;
v |= v >> 8;
v |= v >> 16;
v++;

上式等价于1U << (lg(v-1) + 1)。注意上式在v为0时返回0,该值不是2的幂。你可以在最后加上v += (v==0)来修正这结果。该方法比用查表法求v以2为底的对数法要少两个操作数,且可移植行更好。该方法将最高位的1拷贝到所有较低的bits,然后加1,进位后所有的低位变成了0,最高位上一位则变成1。如果v原本是2的幂,会先减1再向上取到原值。

交错和并bits(简单方法)

unsigned short x;   // x贡献z的偶数上的bits
unsigned short y;   // y贡献z的奇数上的bits
unsigned int z = 0; // z 为输出的结果,也叫Morton Number

for (int i = 0; i < sizeof(x) * CHAR_BIT; i++) // 非循环方式将更快
{
  z |= (x & 1U << i) << i | (y & 1U << i) << (i + 1);
}

交错bits(也叫 Morton numbers)在线性2D整型坐标中非常有用,x和y的值被和并到一个数中,Morton Number之差越小表示两点间的x和y越接近。

交错和并bits(查表法)

static const unsigned short MortonTable256[256] = 
{
  0x0000, 0x0001, 0x0004, 0x0005, 0x0010, 0x0011, 0x0014, 0x0015, 
  0x0040, 0x0041, 0x0044, 0x0045, 0x0050, 0x0051, 0x0054, 0x0055, 
  0x0100, 0x0101, 0x0104, 0x0105, 0x0110, 0x0111, 0x0114, 0x0115, 
  0x0140, 0x0141, 0x0144, 0x0145, 0x0150, 0x0151, 0x0154, 0x0155, 
  0x0400, 0x0401, 0x0404, 0x0405, 0x0410, 0x0411, 0x0414, 0x0415, 
  0x0440, 0x0441, 0x0444, 0x0445, 0x0450, 0x0451, 0x0454, 0x0455, 
  0x0500, 0x0501, 0x0504, 0x0505, 0x0510, 0x0511, 0x0514, 0x0515, 
  0x0540, 0x0541, 0x0544, 0x0545, 0x0550, 0x0551, 0x0554, 0x0555, 
  0x1000, 0x1001, 0x1004, 0x1005, 0x1010, 0x1011, 0x1014, 0x1015, 
  0x1040, 0x1041, 0x1044, 0x1045, 0x1050, 0x1051, 0x1054, 0x1055, 
  0x1100, 0x1101, 0x1104, 0x1105, 0x1110, 0x1111, 0x1114, 0x1115, 
  0x1140, 0x1141, 0x1144, 0x1145, 0x1150, 0x1151, 0x1154, 0x1155, 
  0x1400, 0x1401, 0x1404, 0x1405, 0x1410, 0x1411, 0x1414, 0x1415, 
  0x1440, 0x1441, 0x1444, 0x1445, 0x1450, 0x1451, 0x1454, 0x1455, 
  0x1500, 0x1501, 0x1504, 0x1505, 0x1510, 0x1511, 0x1514, 0x1515, 
  0x1540, 0x1541, 0x1544, 0x1545, 0x1550, 0x1551, 0x1554, 0x1555, 
  0x4000, 0x4001, 0x4004, 0x4005, 0x4010, 0x4011, 0x4014, 0x4015, 
  0x4040, 0x4041, 0x4044, 0x4045, 0x4050, 0x4051, 0x4054, 0x4055, 
  0x4100, 0x4101, 0x4104, 0x4105, 0x4110, 0x4111, 0x4114, 0x4115, 
  0x4140, 0x4141, 0x4144, 0x4145, 0x4150, 0x4151, 0x4154, 0x4155, 
  0x4400, 0x4401, 0x4404, 0x4405, 0x4410, 0x4411, 0x4414, 0x4415, 
  0x4440, 0x4441, 0x4444, 0x4445, 0x4450, 0x4451, 0x4454, 0x4455, 
  0x4500, 0x4501, 0x4504, 0x4505, 0x4510, 0x4511, 0x4514, 0x4515, 
  0x4540, 0x4541, 0x4544, 0x4545, 0x4550, 0x4551, 0x4554, 0x4555, 
  0x5000, 0x5001, 0x5004, 0x5005, 0x5010, 0x5011, 0x5014, 0x5015, 
  0x5040, 0x5041, 0x5044, 0x5045, 0x5050, 0x5051, 0x5054, 0x5055, 
  0x5100, 0x5101, 0x5104, 0x5105, 0x5110, 0x5111, 0x5114, 0x5115, 
  0x5140, 0x5141, 0x5144, 0x5145, 0x5150, 0x5151, 0x5154, 0x5155, 
  0x5400, 0x5401, 0x5404, 0x5405, 0x5410, 0x5411, 0x5414, 0x5415, 
  0x5440, 0x5441, 0x5444, 0x5445, 0x5450, 0x5451, 0x5454, 0x5455, 
  0x5500, 0x5501, 0x5504, 0x5505, 0x5510, 0x5511, 0x5514, 0x5515, 
  0x5540, 0x5541, 0x5544, 0x5545, 0x5550, 0x5551, 0x5554, 0x5555
};

unsigned short x; // x贡献z的偶数上的bits
unsigned short y; // y贡献z的奇数上的bits
unsigned int z;   // z 为输出的结果,也叫Morton Number

z = MortonTable256[y >> 8]   << 17 | 
    MortonTable256[x >> 8]   << 16 |
    MortonTable256[y & 0xFF] <<  1 | 
    MortonTable256[x & 0xFF];

如果用一个额外的表存储y的查找表可以减少一个移位操作,但是会需要额外的内存。同理,如果使用4个查找表表示上述4行的结果可以将操作数减少至11个。

交错和并bits(使用神奇的2进制数)

static const unsigned int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF};
static const unsigned int S[] = {1, 2, 4, 8};

unsigned int x; // x贡献z的偶数上的bits
unsigned int y; // y贡献z的奇数上的bits
unsigned int z; // z 为输出的结果,也叫Morton Number
                // x 和 y 必须小于65536

x = (x | (x << S[3])) & B[3];
x = (x | (x << S[2])) & B[2];
x = (x | (x << S[1])) & B[1];
x = (x | (x << S[0])) & B[0];

y = (y | (y << S[3])) & B[3];
y = (y | (y << S[2])) & B[2];
y = (y | (y << S[1])) & B[1];
y = (y | (y << S[0])) & B[0];

z = x | (y << 1);

递推求下一个bit排列

假设我们有N bits为 1的一个整型,我们想求下一个(大于该数)满足有N bits 为 1 的整型排列。例如,如果N为3,当前整型为00010011,下一个则是00010101, 00010110, 00011001,00011010, 00011100, 00100011等等。下面是求该排列的快速方法:

unsigned int v; // 当前排列
unsigned int w; // 下一个排列

unsigned int t = v | (v - 1); // t 等于 v 将低位的0全部置1

w = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));  

内建的__builtin_ctz(v) 函数在GNU C编译器(x86)中返回该数尾部的0 bits 的个数,如果你用微软的编译器,可以使用_BitScanForward方法。如果没有的话,可以使用前面提到求尾部连续0 bits的个数的方法。
下面另有一个使用了除法稍慢的方法,但是不用计算尾部0 bits的个数。

unsigned int t = (v | (v - 1)) + 1;  
w = t | ((((t & -t) / (v & -v)) >> 1) - 1);