Please enable JavaScript to view the comments powered by Disqus.

位运算技巧

深坑待填!

本文翻译自这篇文章,原作者Sean Eron Anderson. 如果有翻译错误的地方,欢迎指正。另外,由于年代久远,文中给出的链接可能已经失效。

UPDATE 2016/09/30: 继续翻译,并改进已翻译文本

开头语

本文所有代码片段皆属公有(除非特别指明)——您想怎么用就怎么用。所有的文本皆属Sean Eron Anderson版权所有(1997-2005)。本着希望帮助他人的想法,我将代码及其解释公之于众,但若要这些代码用于特定目的,请记住我并未明显或隐含地保证其可用和适用性。截止于2005年五月5日,所有代码均经过了完整的测试,并经上千人阅览。此外,卡内基·梅隆大学计算机学院院长Randal Bryant教授使用他的轮子uclid代码测试系统进行了测试,测试内容几乎涵盖了所有方面。他没有覆盖到的情况,我已经在32位机器上针对所有可能输入补充测试了一遍。对于第一个找到我代码里正当bug的人,我将给予他10美元的奖金(通过支票或PayPal)。如果将这笔钱转交慈善机构的话,增至20美元。

关于操作统计的说明

统计本文算法中操作数量的时候,所有C语言运算符皆看做是一次操作。不涉及写入内存的中间过程,不在统计范围内。当然,这种统计方式只能是对真实的机器指令和CPU时间的一个近似。在此假定所有的操作耗费相同的时间。这虽然与真实情况不符,但的确是CPU发展努力的方向。一个代码样本在不同系统上运行时间会有很多细微的差别,可能是缓存大小、内存带宽、指令集等等因素造成的。总之,拿去跑一跑是确定方法之间快慢的最好方法,所以在实际机器上测试的时候,我所说的仅供参考。

获取一个整数的正负号

1
2
int v;        //我们想知道v的符号
int sign; //用来存放结果
1
2
3
4
5
6
// CHAR_BIT 代表一字节(byte)是多少位(bit)(一般是8)
sign = -(v < 0); //如果v<0,得到-1,否则是0
// 或者,为了避免在用标记寄存器的CPU (IA32) 上使用转移指令:
sign = -(int)((unsigned int)((int)v) >> (sizeof(int) * CHAR_BIT - 1));
// 或者,这样可以少一个机器指令(但牺牲了可移植性):
sign = v >> (sizeof(int) * CHAR_BIT - 1);

最后一个表达式,对于32位整数求得sign = v >> 31,这比直接求法sign = -(v < 0)少了一个操作。该方法之所以可行,是因为当我们把有符号整数向右移位时,最左边一位的值被复制到了其他位上。对一个有符号整数来说,最左边一位是1表示是负数,其他情况是0;所有位上全是1代表-1. 可惜的是,上述行为是与架构相关的。(译注:带符号数负数的左移在C99标准中属于未定义行为)

另外,如果你希望sign的结果是1-1的话:

1
sign = +1 | (v >> (sizeof(int) * CHAR_BIT - 1));  // v < 0得到 -1, 其余情况是 +1

其次,你希望结果是-1, 01的话:

1
2
3
4
5
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, 使用:

1
sign = 1 ^ ((unsigned int)v >> (sizeof(int) * CHAR_BIT - 1)); // v < 0得到0, 其余是1

作者注:

  • 2003年3月7日,Angus Duggan指出1989 ANSI C标准规定有符号整型的右移运算结果取决于具体实现,因此在某些系统上靠右移的方法可能失效。
  • 为了使代码可移植性更好,2005年9月28日,Toby Speight建议我使用 CHAR_BIT 而不是指定一个字节8位长。2006年3月4日,Angus推荐了上面的使用了类型转换的方法,它可移植性更好。
  • 2009年9月12日,Rohit Garg给出了检测非负数的方法。

检测两个整数是否异号

1
int x, y;               // 待检测的两个数
1
bool f = ((x ^ y) < 0); // 当且仅当x与y符号相反的时候为真

作者注:

  • 2009年11月26日的时候Manfred Weis建议我加上这个条目。

不利用分支来计算整数的绝对值

1
2
3
int v;           // 我们现在要求v的绝对值
unsigned int r; // 输出结果
int const mask = v >> sizeof(int) * CHAR_BIT - 1;
1
r = (v + mask) ^ mask;

另一个专利所有的变体:

1
r = (v ^ mask) - mask;

某些CPU没有绝对值指令(或者编译器未能使用到它)。在分支指令开销很大的机器上,上面的表达式能够比直接求法r = (v < 0) ? -(unsigned)v : v快很多,尽管操作的数量一样。

作者注:

  • 2003年3月7日,Angus Duggan指出1989 ANSI C标准规定有符号整型的右移运算结果取决于具体实现,因此在某些系统上靠右移的方法可能失效。我拜读了ANSI C, 其中并没有要求将数值表示为补码,所以说这也可能是方法失效的原因(在那些逐渐退出历史舞台的老爷机上仍然采用反码表示)。
  • 2004年3月14日,Keith H. Duggar告诉了我上面的变体,比我刚开始想出来的那个r=(+1|(v>>(sizeof(int)*CHAR_BIT-1)))*v要好,因为少用了个乘法。不幸的是,那个方法于2000年6月6日由Vladimir Yu Volkonsky在美国申请了专利,归属于Sun Microsystems.
  • 2006年8月13日,Yuriy Kaminskiy告诉我专利可能失效因为在专利申请之前它就已经被发布出去了,比如早在1996年11月9日Agner Fog就已经发表的How to Optimize for the Pentium Processor. Yuriy还告诉我那篇文档在1997年翻译成了俄文,而Vladimir可能读过。此外,互联网档案也有个指向它的旧链接
  • 2007年1月30日,Peter Kankowski对我分享了他受Microsoft Visual C++编译器输出启发而想到的绝对值版本。该方法被拿来作为本题的主要解决方案。
  • 2007年12月6日,Hai Jin抱怨说结果是有符号的,所以计算大多数负数绝对值的时候结果还是负的。
  • 2008年4月15日,Andrew Shapira指出那个直接求法可能会溢出,因为缺少一个(unsigned)的类型转换。为了尽可能满足移植性要求,他给我的改进求法是(v < 0) ? (1 + ((unsigned)(-1-v))) : (unsigned)v.
  • 2008年7月9日,Vincent Lefèvre引用ISO C99标准说服了我,让把上面那个删掉,因为即使不是使用补码的机器-(unsigned)v也能得到正确的结果。计算-(unsigned)v的时候,首先加上2**N来将负值v转化成无符号值,得到v的补码表示(记做U)。然后,取U的相反数,得到结果:-U = 0 - U = 2**N - U = 2**N - (v + 2**N) = -v = abs(v).

译注:2**N 表示 $2^N$.

不利用分支来计算两个整数的最大(max)最小(min)值

1
2
3
int x;
int y; // 我们要计算x, y的最大最小值
int r; // 存放结果
1
r = y ^ ((x ^ y) & -(x < y)); // min(x, y)

这个世界上还有使用条件分支非常耗费时间的老爷机,而且没有条件移动指令。在这些机器上上面的方法可能会比直接求法r = (x < y) ? x : y要快,尽管位运算解法涉及了两个以上的指令(虽然一般说来,直接解法是最好的)。 解法证明如下:如果x < y, 那么-(x < y)的各二进制位就全是1, 所以r = y ^ ((x ^ y) & -(x < y)) = y ^ ((x ^ y) & ~0 = y ^ x ^ y = x. 另外,如果x >= y, 那么-(x < y)就全是0, 所以r = y ^ ((x ^ y) & 0) = y. 在某些机器上,计算(x < y)要得到0或1的话需要分支指令,因此该解法不具备任何优势。

求最大值的话,使用:

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

快而脏的解法

如果你知道INT_MIN <= x - y <= INT_MAX的话,那你就可以用下面这种解法,由于它只需要计算一次(x - y)所以会快一点。

1
2
r = y + ((x - y) & ((x - y) >> (sizeof(int) * CHAR_BIT - 1))); // min(x, y)
r = x - ((x - y) & ((x - y) >> (sizeof(int) * CHAR_BIT - 1))); // max(x, y)

由于1989 ANSI C未规定有符号数右移的结果,所以该解法不可移植。如果溢出时抛出异常,那么减法时xy的值就是无符号型或者被转换成无符号型来防止不必要的异常抛出。然而这里右移运算需要一个有符号的操作数,使得右移负数时填充1,所以这里转换成有符号整型。

作者注:

  • 2003年3月7日,Angus Duggan指出右移的可移植性问题。
  • 2005年5月3日,Randal E. Bryant警告我说需要一个INT_MIN <= x - y <= INT_MAX的前提条件,并给出了一个慢而脏的修改版本。以上只在这个快而脏的解法里考虑了。
  • 2005年7月6日,Nigel Horspoon观察到在奔腾上,由于gcc对(x < y)的求值问题,它生成了和直接解法一样的代码。
  • 2008年7月9日,Vincent Lefèvre指出r = y + ((x - y) & -(x < y))里的减法有潜在的溢出风险(这是我前一版的代码)。
  • 2009年6月2日,Timothy B. Terriberry建议我使用异或而不是加减来避免溢出的危险。

检测一个数是不是2的次幂

1
2
unsigned int v; // 要检测的数是v
bool f; // 存放结果
1
f = (v & (v - 1)) == 0;

注意0并不是2的多少次幂,因此做如下修改:

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

符号扩展

从给定位宽扩展

符号扩展对于内建类型来说是自动的,比如从charint. 但是,假设你有一个以b位补码表示的有符号整数x, 你想把x转换成int型的,此时存储空间大于b位。只是简单地复制的话,当x是正数的时候还可以,是负数的时候就不对了,必须进行符号扩展。例如,假设我们只用4位来存储数字,那么-3表示成二进制就是1101. 如果用8位来存,那么就是11111101. 当扩展一个负数的时候需要将扩展的高位全赋为1, 这就是符号扩展。在C语言里,符号扩展就是小事一桩,因为可以在structunion里使用位域。例如,要将5位长整数扩展到int型:

1
2
int x; // 将它从5位长转换到int
int r; // 存放结果
1
2
struct {signed int x:5;} s;
r = s.x = x;

下面是一个使用了相同语言特性的C++模板函数,仅用一个操作,就能从b位长扩展(尽管编译器生成的肯定不止这么点):

1
2
3
4
5
6
7
8
template <typename T, unsigned B>
inline T signextend(const T x)
{
struct {T x:B;} s;
return s.x = x;
}

int r = signextend<signed int,5>(x); // 从5位长的x扩展到int型的r

作者注:

  • 2005年5月2日,John Byrd发现了代码里一处拼写错误(都是HTML的锅)。
  • 2006年3月4日,Pat Wood指出ANSI C标准要求位域需要关键字signed来表示有符号数,否则其符号未定义。

从不定长位宽扩展

有时我们需要扩展的数并不知道它的位长b. 或者我们用的是一门没有位域特性的语言,比如 Java.

1
2
3
4
5
6
7
unsigned b; // 表示x的位长
int x; // 扩展这个b位长的数x
int r; // 存放结果
int const m = 1U << (b - 1); // 如果给出b的话,位掩码可以提前算出来

x = x & ((1U << b) - 1); // 如果x的b位以上已经是0的话,跳过这一步
r = (x ^ m) - m;

上面的代码需要四次操作,但如果位宽给定的话,我们就只需要两个高效率的操作(假设高位已经全是0)。

一个牺牲了可移植性的稍快点的方法如下,这种方法并不要求xb位以上的数字是0:

1
2
int const m = CHAR_BIT * sizeof(x) - b;
r = (x << m) >> m;

作者注:

  • 2004年6月13日,Sean A. Irvine希望我能在本页加上这些符号扩展的方法。他给了我最初的版本m = (1 << (b - 1)) - 1; r = -(x & ~m) | x, 后来被我优化成了m = 1U << (b - 1); r = -(x & m) | x.
  • 然而2007年5月11日,Shay Green给了我上面解法里的版本,比我的还少一次操作。
  • 2008年10月15日,Vipin Sharma给我改进建议,考虑到了x的高位可能有1的情况。
  • 2009年12月31日,Chris Pirazzi, 建议我加上下面要讲的那种解法。

从不定长位宽扩展,3次操作

下面讲的这个解法可能在某些机器上较慢,因为它需要乘除法。该版本需要4次操作。如果你知道位宽b比1大,请使用r = (x * multipliers[b]) / multipliers[b]完成这种类型的符号扩展,它仅需3次操作,一次查表。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
unsigned b; // 表示x的位长
int x; // 扩展这个b位长的数x
int r; // 存放结果

#define M(B) (1U << ((sizeof(x) * CHAR_BIT) - B)) // CHAR_BIT=bits/byte
static int const multipliers[] =
{
0, M(1), M(2), M(3), M(4), M(5), M(6), M(7),
M(8), M(9), M(10), M(11), M(12), M(13), M(14), M(15),
M(16), M(17), M(18), M(19), M(20), M(21), M(22), M(23),
M(24), M(25), M(26), M(27), M(28), M(29), M(30), M(31),
M(32)
}; // (如果大于64位,继续添加)
static int const divisors[] =
{
1, ~M(1), M(2), M(3), M(4), M(5), M(6), M(7),
M(8), M(9), M(10), M(11), M(12), M(13), M(14), M(15),
M(16), M(17), M(18), M(19), M(20), M(21), M(22), M(23),
M(24), M(25), M(26), M(27), M(28), M(29), M(30), M(31),
M(32)
}; // (同上)
#undef M

r = (x * multipliers[b]) / divisors[b];

下面的这个变体并不可移植,但在使用算数右移的架构上,能保留符号,应该能快一点。

1
2
const int s = -b; // 或者:  sizeof(x) * CHAR_BIT - b;
r = (x << s) >> s;

作者注:

  • 2005年5月3日Randal E. Bryant发现了早期版本里的一个bug(把multipliers[]用作了divisors[]),导致了当x=1,b=1的时候出错。

选择性地置位清零,不利用分支

1
2
3
4
5
6
7
8
bool f;         // 表示条件
unsigned int m; // 位掩码
unsigned int w; // 待修改的输入:if (f) w |= m; else w &= ~m;

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

// 或,针对超标量CPU架构:
w = (w & ~m) | (-f & m);

作者注:

  • 在某些架构上,没有条件分支要强于使用多两倍的操作数量来计算分支。例如,据非正式测试显示,在AMD Athlon™ XP 2100+上能快上5-10%. Intel Core 2 Duo运行超标量版比普通版快了约16%.
  • 2003年12月11日,Glenn Slayden告诉了我第一个版本。
  • 2007年4月3日,Marco Yu告诉了我超标量版并于两天后指出了两个拼写错误。

选择性地取负,不利用分支

如果你希望仅当条件为假的时候取负,下面是不涉及条件分支的版本:

1
2
3
bool fDontNegate;  // 标识我们什么时候不对v取反
int v; // 待操作的数v
int r; // result = fDontNegate ? v : -v;
1
r = (fDontNegate ^ (fDontNegate - 1)) * v;

如果你希望当条件为真时操作:

1
2
3
bool fNegate;  // 标识我们什么时候对v取反
int v; // 待操作的数v
int r; // result = fNegate ? -v : v;
1
r = (v ^ -fNegate) + fNegate;

作者注:

  • 2009年6月2日,Avraham Plotnitzky建议我加上第一个版本。
  • 我又想避免使用乘法,因此在2009年6月8日想出了第二个版本。
  • 2009年11月26日,Alfonso De Gregorio指出我丢了几个括号,并收到了bug奖励。

根据掩码合并两个数

1
2
3
4
unsigned int a;    // 合并未屏蔽位
unsigned int b; // 合并屏蔽位
unsigned int mask; // 1表示要合并b的对应位; 0针对a.
unsigned int r; // (a & ~mask) | (b & mask) 的结果
1
r = a ^ ((a ^ b) & mask);

相比合并两个数的直接解法,我的比它少了一个操作。如果掩码是常量的话,这种方法不占任何优势。

作者注:

  • 2006年2月9日,Ron Jeffery给了我这个。

统计数位中1的个数

简单粗暴的解法

1
2
3
4
5
6
7
unsigned int v; // 待统计的数
unsigned int c; // 用来累加计数

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

这种解法每一位一次循环,直到没有1. 所以对一个只有最高位是1的32位数来说,它会循环32次。

查表法

1
2
3
4
5
6
7
8
9
10
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; // 待统计的数
unsigned int c; // 统计1的个数
1
2
3
4
c = BitsSetTable256[v & 0xff] + 
BitsSetTable256[(v >> 8) & 0xff] +
BitsSetTable256[(v >> 16) & 0xff] +
BitsSetTable256[v >> 24];

或者:

1
2
3
4
5
unsigned char * p = (unsigned char *) &v;
c = BitsSetTable256[p[0]] +
BitsSetTable256[p[1]] +
BitsSetTable256[p[2]] +
BitsSetTable256[p[3]];
1
2
3
4
5
6
// 用算法来初始化数组
BitsSetTable256[0] = 0;
for (int i = 0; i < 256; i++)
{
BitsSetTable256[i] = (i & 1) + BitsSetTable256[i / 2];
}

作者注:

  • 2009年7月14日Hallvard Furuseth建议我使用宏来生成表。

Brian Kernighan的方法

1
2
3
4
5
6
unsigned int v;
unsigned int c;
for (c = 0; v; c++)
{
v &= v - 1; // 清除最低非0位
}

这种方法的循环次数和1的个数一样多。所以如果有一个只有最高位是1的数,只会循环一次。

作者注:

  • 该方法发表于1988, the C Programming Language 2nd Ed (作者Brian W. Kernighan和Dennis M. Ritchie). 书中练习2-9提到了这种方法。
  • 2006年4月19日,Don Knuth指出这种方法是由Peter Wegner在 CACM 3 (1960), 322第一次发表。(同时也由Derrick Lehmer指出该方法在1964年Beckenbach编过的一本书里发表过。)

针对14, 24, 或32位长数字,使用64位指令集

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
unsigned int v;
unsigned int c;

// 针对至多14位长的v:
c = (v * 0x200040008001ULL & 0x111111111111111ULL) % 0xf;

// 针对至多24位长的v:
c = ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL)
% 0x1f;

// 针对至多32位长的v:
c = ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) %
0x1f;
c += ((v >> 24) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;

该方法需要能快速求余的64位CPU才能高效率。14位长的只需要3次操作;24位的需要10次;32位15次。

作者注:

  • Rich Schroeppel最初原创了针对9位的求法,类似于我这里的14位的方法;详见Beeler, M., Gosper, R. W., and Schroeppel, R. HAKMEM. MIT AI Memo 239, Feb. 29, 1972.的Programming Hacks部分。它启发了Sean Anderson(译注:作者自己), 并想出了上面的几种方法。
  • 2005年5月3日,Randal E. Bryant找出了一堆bug.
  • 2007年2月1日,Bruce Dawson对12位方法做了微调,使之变成了现在的14位版本。

并行方法

1
2
3
4
5
6
7
8
9
10
unsigned int v; // 要统计v(32位数值)
unsigned int c; // 存放结果
static const int S[] = {1, 2, 4, 8, 16}; // 二进制魔法数字
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表示为二进制,就是:

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

对于大点的数据,我们可以延长魔法数字BS来适配。如果数字有k位,数组SB就应该有$\lceil \lg(k) \rceil$个元素,而且我们要写出关于c的同样多个数的表达式。对于一个32位长的v, 进行了16次操作。

对于一个32位整数v的最佳解法是:

1
2
3
v = v - ((v >> 1) & 0x55555555);                    // 用v来存储中间值
v = (v & 0x33333333) + ((v >> 2) & 0x33333333); // 同上
c = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24; // 计数

该解法只用了12次操作,与查表法的相同,并避免使用额外内存和缓存未命中表的情况。该解法混合了上面并行法和更前面的使用乘法的解法(使用64位指令集的那个解法),然而它并没有用到64位指令集。统计在并行部分完成,计数求和由乘以0x1010101并右移24位完成。

一个最佳通用解法如下,它最多支持到128位(参数化表示为T):

1
2
3
4
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; // 计数

作者注:

  • 更多信息见Ian Ashdown的贴子
  • 最佳解法是2005年10月5日由Andrew Shapira告诉我的;他在Software Optimization Guide for AMD Athlon™ 64 and Opteron™ Processors的187-188页找到了它。
  • 2005年12月14日Charlie Gordon给了我一个能优化掉一次操作的建议。
  • 2005年12月30日Don Clugston从中整理出了三种方法。
  • 我接受了建议后又犯了几个拼写错误,于2006年1月8日由Eric Cole指出。此后,他给了上面那个通用解法的建议。
  • 2007年4月5日,Al Williams告知我在最开始的方法里有一句冗余代码。

从最高位到指定位

下面给出求从最高位到指定位1的个数(rank)的解法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
uint64_t v;       // 待计算的数v
unsigned int pos; // 指定的位置
uint64_t r; // 存放结果

// 将给定位置以下的数字删除.
r = v >> (sizeof(v) * CHAR_BIT - pos);
// 并行计算
// r = (r & 0x5555...) + ((r >> 1) & 0x5555...);
r = r - ((r >> 1) & ~0UL/3);
// r = (r & 0x3333...) + ((r >> 2) & 0x3333...);
r = (r & ~0UL/5) + ((r >> 2) & ~0UL/5);
// r = (r & 0x0f0f...) + ((r >> 4) & 0x0f0f...);
r = (r + (r >> 4)) & ~0UL/17;
// r = r % 255;
r = (r * (~0UL/255)) >> ((sizeof(v) - 1) * CHAR_BIT);

作者注:

  • 2009年11月21日,Juha Järvi给了我这个下面问题的镜像问题。

已知rank求出位数

下面的64位代码计算出从左至右第r1. 也就是说,从最高位向右统计1的个数,直到得到rankr, 返回对应的这个位。如果给出的rank超出了范围,返回64. 代码可以针对32位或者从左往右计算做出修改。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
uint64_t v;          // 待计算的数v
unsigned int r; // 输入rank, 范围是1~64.
unsigned int s; // 输出rank对应的位数
uint64_t a, b, c, d; // 中间变量
unsigned int t; // 位数统计

// 对64位整数做并行处理,中间结果都存储起来
// a = (v & 0x5555...) + ((v >> 1) & 0x5555...);
a = v - ((v >> 1) & ~0UL/3);
// b = (a & 0x3333...) + ((a >> 2) & 0x3333...);
b = (a & ~0UL/5) + ((a >> 2) & ~0UL/5);
// c = (b & 0x0f0f...) + ((b >> 4) & 0x0f0f...);
c = (b + (b >> 4)) & ~0UL/0x11;
// d = (c & 0x00ff...) + ((c >> 8) & 0x00ff...);
d = (c + (c >> 8)) & ~0UL/0x101;
t = (d >> 32) + (d >> 48);
// 开始不使用条件分支的计算
s = 64;
// if (r > t) {s -= 32; r -= t;}
s -= ((t - r) & 256) >> 3; r -= (t & ((t - r) >> 8));
t = (d >> (s - 16)) & 0xff;
// if (r > t) {s -= 16; r -= t;}
s -= ((t - r) & 256) >> 4; r -= (t & ((t - r) >> 8));
t = (c >> (s - 8)) & 0xf;
// if (r > t) {s -= 8; r -= t;}
s -= ((t - r) & 256) >> 5; r -= (t & ((t - r) >> 8));
t = (b >> (s - 4)) & 0x7;
// if (r > t) {s -= 4; r -= t;}
s -= ((t - r) & 256) >> 6; r -= (t & ((t - r) >> 8));
t = (a >> (s - 2)) & 0x3;
// if (r > t) {s -= 2; r -= t;}
s -= ((t - r) & 256) >> 7; r -= (t & ((t - r) >> 8));
t = (v >> (s - 1)) & 0x1;
// if (r > t) s--;
s -= ((t - r) & 256) >> 8;
s = 65 - s;

如果分支在你的CPU上速度很快的话,取消使用注释里的条件语句进行替换。

奇偶校验

译注:信息是以比特流的方式传输的,类似01000001. 在传输过程中,有可能会发生错误,比如,我们存储了01000001,但是取出来却是01000000,即低位由0变成了1。为了检测到这种错误,我们可以通过“奇偶校验”来实现。假如,我们存储的数据是一个字节,8个比特位,那我们就可以计算每个字节比特位是1的个数,如果是偶数个1,那么,我们就把第九个位设为1,如果是奇数个1,那么就把第九个位设为0,这样连续9个字节比特位为1的位数肯定是奇数。这种方法叫做“奇校验”,“偶校验”和此类似。当然,在实际应用中,也可以把一个字节的前7位作为数据位,最后一个为作为校验位。

简单粗暴的解法

1
2
3
4
5
6
7
8
unsigned int v;       // 待检验的数v
bool parity = false; // 存放结果,奇数个为true

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

该解法使用了类似于前述Brian Kernigan解法的思路,消耗时间与二进制1的数量成正比。

查表法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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; // 待计算的数b
bool parity = ParityTable256[b];

// 针对32位长的数:
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]];

作者注:

  • 2005年5月3日Randal E. Bryant支持我加上那个很直接的解法(最后那个变体)。
  • 2005年9月27日,Bruce Rawles找到了几个表名的拼写错误,并因此收到了$10的找bug奖金。
  • 2006年10月9日,Fabrice Bellard给了我上面的32位长的版本,它只需一次查表,原来的版本需要4次查表,更慢。
  • 2009年7月14日Hallvard Furuseth建议我使用宏来生成表。

使用64位指令集的乘法和求余

1
2
unsigned char b;  // 待计算的数b
bool parity = (((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;

上面的方法大概需要4次操作,但只对整字节有效。

使用一次乘法

以下方法能仅用8次操作处理一个32位长的数值,而且只用了一次乘法。

1
2
3
4
5
unsigned int v; // 32-bit word
v ^= v >> 1;
v ^= v >> 2;
v = (v & 0x11111111U) * 0x11111111U;
return (v >> 28) & 1;

对于64位长的,8次操作也足够:

1
2
3
4
5
unsigned long long v; // 64-bit word
v ^= v >> 1;
v ^= v >> 2;
v = (v & 0x1111111111111111UL) * 0x1111111111111111UL;
return (v >> 60) & 1;

作者注:

  • Andrew Shapira想出了这个方法并于2007年9月2日给了我。

并行法

1
2
3
4
5
6
unsigned int v;
v ^= v >> 16;
v ^= v >> 8;
v ^= v >> 4;
v &= 0xf;
return (0x6996 >> v) & 1;

以上方法需要大概9次操作,对32位的有效。对于整字节也可以删掉unsigned int v后面的两行来优化到只用5次操作。该方法首先使用移位和异或对v的每段(每段8bit)进行处理,异或后只取后4个bit, 然后二进制数0110 1001 1001 0110(16进制是0x6996)右移这么多位。这就好像是拿v的后四位去小型的16位奇偶表里去查。v的奇偶结果以二进制1的形式保存并返回。

作者注:

  • 感谢Mathew Hendry于2002年12月15日告诉我的移位-查表方法。优化后能少两次操作并且只用了移位和异或。

交换两个数的值

使用加减法

译注:不要把右值往a, b里代入。

1
2
#define SWAP(a, b) ((&(a) == &(b)) || \
(((a) -= (b)), ((b) += (a)), ((a) = (b) - (a))))

该方法能不使用额外的临时变量来交换ab的值。开始的那句是检查ab是否使用了同一块内存空间,是的话那么该方法结果就不对了(编译器可能会把这个检查优化掉)。如果你开启了溢出异常检测,那么要传递无符号值以避免触发异常。下面的异或方法在某些机器上可能会快一点。不要对浮点数也这么操作,除非你操作的是它们的二进制值。

作者注:

  • 2007年6月12日Sanjeev Sivasankaran建议我加上这一条。
  • 2008年7月9日Vincent Lefèvre指出它有潜在的溢出风险。

使用异或

译注:同上。

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

这是一种很古老的,不使用额外临时变量的交换技巧。

作者注:

  • 2005年1月20日,Iain A. Fleming指出该方法当作用于同一内存空间的时候失效,比如我当i==j的时候SWAP(a[i], a[j]). 所以为避免这种情况的发生,把宏定义改成(((a) == (b)) || (((a) ^= (b)), ((b) ^= (a)), ((a) ^= (b)))).
  • 2009年7月14日,Hallvard Furuseth告诉我说,在某些机器上(((a) ^ (b)) && ((b) ^= (a) ^= (b), (a) ^= (b)))会更快一点,因为(a) ^ (b)的值能重复利用。

交换指定的某几个二进制位

1
2
3
4
5
6
7
unsigned int i, j; // 给定位置
unsigned int n; // 想交换多少位
unsigned int b; // 被交换数b
unsigned int r; // 交换结果

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

举个例子,假如我们有个数b = 00101111(二进制表示),我们想把它里面从i = 1开始的连续n = 3位与从j = 5开始的连续n位交换,其结果应该是r = 11100011(二进制表示)。

该方法的思路与异或法相似。x用于保存我们想要交换的那几位,然后这几位的值被赋值为它们自己与x的异或结果。当然,如果序列超出范围其结果未定义。

作者注:

  • 2009年7月14日,Hallvard Furuseth建议我把1 << n改成1U << n来强行赋值为无符号型来避免对意外地有符号数移位。

倒置二进制位

直接解法

1
2
3
4
5
6
7
8
9
10
11
unsigned int v;     // 待处理的数v
unsigned int r = v; // r存放结果
int s = sizeof(v) * CHAR_BIT - 1; // 最后一步的额外处理

for (v >>= 1; v; v >>= 1)
{
r <<= 1;
r |= v & 1;
s--;
}
r <<= s; // 当v的最高位是0的时候移位

作者注:

  • 2004年10月15日,Michael Hoisie指出最初版本有一个bug.
  • 2005年5月3日Randal E. Bryant建议我删掉一个冗余操作。2005年5月18日Behdad Esfabod给了我一个小建议,它能减少一次循环。
  • 2007年2月6日,Liyong Zhou给了我点建议,使得循环只在v不是0的时候运行,而不是对每一位都循环一次。

查表法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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; // 32位值,每次倒置8位
unsigned int c; // v倒置的结果

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

// Option 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]];

第一个方法要用17次操作,第二种需要12次,如果你的CPU能够快速读取和存储数据的话。

作者注:

  • 2009年7月14日Hallvard Furuseth建议我使用宏来生成表。

倒置一比特,仅用3次操作(使用64位的乘法和求余)

1
2
3
unsigned char b; // reverse this (8-bit) byte

b = (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;

乘法操作生成了五分相互独立的拷贝,每份8位,并合并成一个64位值。将该值分为每十位一组,与操作选出哪些位在组中正确的(倒置)位上。乘法和与运算从原数据拷贝到这里来,所以每个组中一定会有一个原数据。比特位的倒置位置与它们在组中的相对位置相同。最后,求得模2^10 - 1的余数,即能把每组中的选中的比特位聚合到一起。它们相互不会覆盖,所以求余操作的底层操作就像逻辑或操作一样。

该方法见于Beeler, M., Gosper, R. W., and Schroeppel, R. HAKMEM. MIT AI Memo 239, Feb. 29, 1972.的Programming Hacks部分。

倒置一比特,仅用4次操作(使用64位的乘法,不使用除法)

1
2
3
unsigned char b; // reverse this byte

b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;

下面就是该方法的原理展示,使用a,b,c,d,e,f,g,h来表示该比特的各个位。注意第一次乘法是如何将这个数拷贝多份的。以及最后一次乘法是如何将这几个位合并到从右往左数第五个比特中去的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
                                                                                        abcd efgh (-> hgfe dcba)
* 1000 0000 0010 0000 0000 1000 0000 0010 (0x80200802)
-------------------------------------------------------------------------------------------------
0abc defg h00a bcde fgh0 0abc defg h00a bcde fgh0
& 0000 1000 1000 0100 0100 0010 0010 0001 0001 0000 (0x0884422110)
-------------------------------------------------------------------------------------------------
0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
* 0000 0001 0000 0001 0000 0001 0000 0001 0000 0001 (0x0101010101)
-------------------------------------------------------------------------------------------------
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

最后两步可以合并,因为某些处理器的寄存器可以按照字节来读写;只用乘法,然后寄存器存储高位的32位,再取低八位,这样就只要6次操作。

该方法于2001年7月13日由Sean Anderson创造而出。

倒置一比特,仅用7次操作(不使用64位)

1
b = ((b * 0x0802LU & 0x22110LU) | (b * 0x8020LU & 0x88440LU)) * 0x10101LU >> 16;

注意要把结果赋值给或转换到unsigned char以去除不需要的垃圾位。

作者注:

  • 2001年7月13日,该方法由Sean Anderson创造而出。
  • 2002年1月3日,Mike Keith指出并修正了拼写错误。

倒置N位,并行法,使用 5*lg(N) 次操作

1
2
3
4
5
6
7
8
9
10
11
12
unsigned int v; // 32-bit word to reverse bit order

// swap odd and even bits
v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
// swap consecutive pairs
v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
// swap nibbles ...
v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
// swap bytes
v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
// swap 2-byte long pairs
v = ( v >> 16 ) | ( v << 16);

下面的变体依然是$O(\lg N)$的,不过它需要更多次操作来倒置v. 它好在即时计算常量来减少内存使用。

1
2
3
4
5
6
7
unsigned int s = sizeof(v) * CHAR_BIT; // bit位数,必须是2的次幂
unsigned int mask = ~0;
while ((s >>= 1) > 0)
{
mask ^= (mask << s);
v = ((v >> s) & mask) | ((v << s) & ~mask);
}

以上方法非常适合于N很大的情形。对于64位或以上的整数,还得按照套路加几行,否则只有低32位被倒置。

作者注:

  • 更多内容详见Dr. Dobb’s Journal 1983, Edwin Freed关于Binary Magic Numbers的文章。
  • 2005年9月13日,Ken Raeburn推荐给我第二种变体。
  • 2006年3月19日,Veldmeijer提出第一种方法的最后一行可以不使用与运算。

求余除法

不使用除法,计算关于 1<<s 的余数

1
2
3
4
5
const unsigned int n;          // 被除数
const unsigned int s;
const unsigned int d = 1U << s; // So d will be one of: 1, 2, 4, 8, 16, 32, ...
unsigned int m; // m will be n % d
m = n & (d - 1);

大部分程序员很早就知道这个技巧了,为了完整性我将其添加进来。

不使用除法,计算关于 (1<<s) - 1 的余数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
unsigned int n;                      // 被除数
const unsigned int s; // s > 0
const unsigned int d = (1 << s) - 1; // so d is either 1, 3, 7, 15, 31, ...).
unsigned int m; // n % d goes here.

for (m = n; n > d; n = m)
{
for (m = 0; n; n >>= s)
{
m += n & d;
}
}
// m值域为[0,d], 但m是余数所以m为d时赋值为0
m = m == d ? 0 : m;

求关于 $2^N-1$ 的余数需要 $5 + (4 + 5 \lceil {N \over s} \rceil) \cdot \lceil \lg{N \over s} \rceil$ 次操作,其中N是被除数的位数。也就是说,时间复杂度为$O(\lg N)$.

作者注:

  • 该方法于2001年8月15日由Sean Anderson想出。
  • 在2004年6月17日Sean A. Irvine指正之前,我错误地评论说可以在最末尾使用m = ((m + 1) & d) - 1;.
  • 2005年4月25日,Michael Miller指出代码中一处拼写错误。

不使用除法,并行计算关于 (1<<s) - 1 的余数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
// The following is for a word size of 32 bits!

static const unsigned int M[] =
{
0x00000000, 0x55555555, 0x33333333, 0xc71c71c7,
0x0f0f0f0f, 0xc1f07c1f, 0x3f03f03f, 0xf01fc07f,
0x00ff00ff, 0x07fc01ff, 0x3ff003ff, 0xffc007ff,
0xff000fff, 0xfc001fff, 0xf0003fff, 0xc0007fff,
0x0000ffff, 0x0001ffff, 0x0003ffff, 0x0007ffff,
0x000fffff, 0x001fffff, 0x003fffff, 0x007fffff,
0x00ffffff, 0x01ffffff, 0x03ffffff, 0x07ffffff,
0x0fffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff
};

static const unsigned int Q[][6] =
{
{ 0, 0, 0, 0, 0, 0}, {16, 8, 4, 2, 1, 1}, {16, 8, 4, 2, 2, 2},
{15, 6, 3, 3, 3, 3}, {16, 8, 4, 4, 4, 4}, {15, 5, 5, 5, 5, 5},
{12, 6, 6, 6 , 6, 6}, {14, 7, 7, 7, 7, 7}, {16, 8, 8, 8, 8, 8},
{ 9, 9, 9, 9, 9, 9}, {10, 10, 10, 10, 10, 10}, {11, 11, 11, 11, 11, 11},
{12, 12, 12, 12, 12, 12}, {13, 13, 13, 13, 13, 13}, {14, 14, 14, 14, 14, 14},
{15, 15, 15, 15, 15, 15}, {16, 16, 16, 16, 16, 16}, {17, 17, 17, 17, 17, 17},
{18, 18, 18, 18, 18, 18}, {19, 19, 19, 19, 19, 19}, {20, 20, 20, 20, 20, 20},
{21, 21, 21, 21, 21, 21}, {22, 22, 22, 22, 22, 22}, {23, 23, 23, 23, 23, 23},
{24, 24, 24, 24, 24, 24}, {25, 25, 25, 25, 25, 25}, {26, 26, 26, 26, 26, 26},
{27, 27, 27, 27, 27, 27}, {28, 28, 28, 28, 28, 28}, {29, 29, 29, 29, 29, 29},
{30, 30, 30, 30, 30, 30}, {31, 31, 31, 31, 31, 31}
};

static const unsigned int R[][6] =
{
{0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000},
{0x0000ffff, 0x000000ff, 0x0000000f, 0x00000003, 0x00000001, 0x00000001},
{0x0000ffff, 0x000000ff, 0x0000000f, 0x00000003, 0x00000003, 0x00000003},
{0x00007fff, 0x0000003f, 0x00000007, 0x00000007, 0x00000007, 0x00000007},
{0x0000ffff, 0x000000ff, 0x0000000f, 0x0000000f, 0x0000000f, 0x0000000f},
{0x00007fff, 0x0000001f, 0x0000001f, 0x0000001f, 0x0000001f, 0x0000001f},
{0x00000fff, 0x0000003f, 0x0000003f, 0x0000003f, 0x0000003f, 0x0000003f},
{0x00003fff, 0x0000007f, 0x0000007f, 0x0000007f, 0x0000007f, 0x0000007f},
{0x0000ffff, 0x000000ff, 0x000000ff, 0x000000ff, 0x000000ff, 0x000000ff},
{0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff},
{0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff},
{0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff},
{0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff},
{0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff},
{0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff},
{0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff},
{0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff},
{0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff},
{0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff},
{0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff},
{0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff},
{0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff},
{0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff},
{0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff},
{0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff},
{0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff},
{0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff},
{0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff},
{0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff},
{0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff},
{0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff},
{0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff}
};

unsigned int n; // numerator
const unsigned int s; // s > 0
const unsigned int d = (1 << s) - 1; // so d is either 1, 3, 7, 15, 31, ...).
unsigned int m; // n % d goes here.

m = (n & M[s]) + ((n >> s) & M[s]);

for (const unsigned int * q = &Q[s][0], * r = &R[s][0]; m > d; q++, r++)
{
m = (m >> *q) + (m & *r);
}
m = m == d ? 0 : m; // OR, less portably: m = m & -((signed)(m - d) >> s);

该方法求关于 $2^N-1$ 的余数至多需要 $O(\lg N)$ 时间,其中N是被除数位长(以上代码是针对32位的)。所需操作次数至多是 $12 + 9 \lceil \lg N \rceil$. 如果你能在编译器知道除数的话,就能去掉表,只需提取相关条目并做循环展开。该方法很容易能扩展到更长位宽。

作者:Dr. A. Clef
发布日期:2015-12-04
修改日期:2017-06-10
发布协议: BY-SA