位运算整理(2)
从一个固定的位宽延伸出来的符号
对于内置类型,如chars和ints,符号扩展是自动的。但是,假设你有一个有符号的二补数x,只用b位存储。此外,假设你想将x转换为一个int,它的位数超过b位。如果x是正数,简单的复制就可以了,但如果是负数,则必须扩展符号。例如,如果我们只有4位来存储一个数字,那么-3在二进制中表示为1101。如果我们有8位,那么-3就是1111101。当我们转换到有更多比特的表示法时,4比特表示法中最重要的比特会被正向复制以填入目标位,这就是符号扩展。在C语言中,从一个恒定的位宽进行符号扩展是很简单的,因为位域可以在结构或联合中指定。例如,要从5位转换为一个完整的整数。
int x; // convert this from using 5 bits to a full int
int r; // resulting sign extended number goes here
struct {signed int x:5;} s;
r = s.x = x;
下面是一个C++模板函数,它使用同样的语言特征,在一次操作中从B位进行转换(当然编译器生成的更多)
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); // sign extend 5 bit number x to r
2006年3月4日,Pat Wood指出,ANSI C标准要求位域必须有关键字 "signed "才能被签名,否则,签名是未定义的。
从一个可变位宽的标志延伸出来
有时我们需要扩展一个数字的符号,但我们不知道先验地知道它的位数,b,用它来表示。(或者我们可能是在像Java这样的语言中进行编程,而这种语言缺乏位数字段。)
unsigned b; // number of bits representing the number in x
int x; // sign extend this b-bit number to r
int r; // resulting sign-extended number
int const m = 1U << (b - 1); // mask can be pre-computed if b is fixed
x = x & ((1U << b) - 1); // (Skip this if bits in x above position b are already zero.)
r = (x ^ m) - m;
上面的代码需要四次操作,但当位宽是常数而不是变量时,只需要两次快速操作,假设上位已经是零。
一种稍快但不太方便携带的方法,它不依赖于位置b以上的x中的位数为零
int const m = CHAR_BIT * sizeof(x) - b;
r = (x << m) >> m;
在3次操作中从可变位宽延伸出的符号
下面的内容在某些机器上可能会很慢,因为乘法和除法需要费力。这个版本是4次操作。如果你知道你的初始位宽b大于1,你可以通过使用r = (x * multipliers[b]) / multipliers[b]在3次操作中完成这种符号扩展,这只需要一次数组查找。
1 unsigned b; // number of bits representing the number in x
2 int x; // sign extend this b-bit number to r
3 int r; // resulting sign-extended number
4 #define M(B) (1U << ((sizeof(x) * CHAR_BIT) - B)) // CHAR_BIT=bits/byte
5 static int const multipliers[] =
6 {
7 0, M(1), M(2), M(3), M(4), M(5), M(6), M(7),
8 M(8), M(9), M(10), M(11), M(12), M(13), M(14), M(15),
9 M(16), M(17), M(18), M(19), M(20), M(21), M(22), M(23),
10 M(24), M(25), M(26), M(27), M(28), M(29), M(30), M(31),
11 M(32)
12 }; // (add more if using more than 64 bits)
13 static int const divisors[] =
14 {
15 1, ~M(1), M(2), M(3), M(4), M(5), M(6), M(7),
16 M(8), M(9), M(10), M(11), M(12), M(13), M(14), M(15),
17 M(16), M(17), M(18), M(19), M(20), M(21), M(22), M(23),
18 M(24), M(25), M(26), M(27), M(28), M(29), M(30), M(31),
19 M(32)
20 }; // (add more for 64 bits)
21 #undef M
22 r = (x * multipliers[b]) / divisors[b];
下面的变化是不可移植的,但在采用算术右移的架构上,保持符号,应该是很快的。
const int s = -b; // OR: sizeof(x) * CHAR_BIT - b; r = (x << s) >> s;
Randal E. Bryant 在 2005 年 5 月 3 日指出了早期版本中的一个错误 (该版本使用乘数[]来表示除数[]),它在 x=1 和 b=1 的情况下失败。
有条件地设置或清除没有分支的位
1 bool f; // conditional flag
2 unsigned int m; // the bit mask
3 unsigned int w; // the word to modify: if (f) w |= m; else w &= ~m;
4
5 w ^= (-f ^ w) & m;
6
7 // OR, for superscalar CPUs:
8 w = (w & ~m) | (-f & m);
有条件地否定一个没有分支的值
如果你只需要在一个标志为假的时候进行否定,那么使用下面的方法来避免分支。
1 bool fDontNegate; // Flag indicating we should not negate v.
2 int v; // Input value to negate if fDontNegate is false.
3 int r; // result = fDontNegate ? v : -v;
4
5 r = (fDontNegate ^ (fDontNegate - 1)) * v;
6 If you need to negate only when a flag is true, then use this:
7 bool fNegate; // Flag indicating if we should negate v.
8 int v; // Input value to negate if fNegate is true.
9 int r; // result = fNegate ? -v : v;
10
11 r = (v ^ -fNegate) + fNegate;
根据掩码合并两个值的位数
1 unsigned int a; // value to merge in non-masked bits 2 unsigned int b; // value to merge in masked bits 3 unsigned int mask; // 1 where bits from b should be selected; 0 where from a. 4 unsigned int r; // result of (a & ~mask) | (b & mask) goes here 5 6 r = a ^ ((a ^ b) & mask);
这就从根据位掩码组合两组比特的明显方式中减去了一个操作。如果掩码是一个常数,那么可能没有任何优势。
计数位集
unsigned int v; // count the number of bits set in v
unsigned int c; // c accumulates the total bits set in v
for (c = 0; v; v >>= 1)
{
c += v & 1;
}
方法要求每一个位进行一次迭代,直到没有更多的位被设置。所以在一个32位的字上,只有高位被设置,它将经历32次迭代。
由查找表设置的计数位
1 static const unsigned char BitsSetTable256[256] =
2 {
3 # define B2(n) n, n+1, n+1, n+2
4 # define B4(n) B2(n), B2(n+1), B2(n+1), B2(n+2)
5 # define B6(n) B4(n), B4(n+1), B4(n+1), B4(n+2)
6 B6(0), B6(1), B6(1), B6(2)
7 };
8
9 unsigned int v; // count the number of bits set in 32-bit value v
10 unsigned int c; // c is the total bits set in v
11
12 // Option 1:
13 c = BitsSetTable256[v & 0xff] +
14 BitsSetTable256[(v >> 8) & 0xff] +
15 BitsSetTable256[(v >> 16) & 0xff] +
16 BitsSetTable256[v >> 24];
17
18 // Option 2:
19 unsigned char * p = (unsigned char *) &v;
20 c = BitsSetTable256[p[0]] +
21 BitsSetTable256[p[1]] +
22 BitsSetTable256[p[2]] +
23 BitsSetTable256[p[3]];
24
25
26 // To initially generate the table algorithmically:
27 BitsSetTable256[0] = 0;
28 for (int i = 0; i < 256; i++)
29 {
30 BitsSetTable256[i] = (i & 1) + BitsSetTable256[i / 2];
31 }
计数位集,布莱恩-克尼汉的方式
unsigned int v; // count the number of bits set in v
unsigned int c; // c accumulates the total bits set in v
for (c = 0; v; c++)
{
v &= v - 1; // clear the least significant bit set
}
Brian Kernighan 的方法是有多少个设置位就有多少次迭代。因此,如果我们有一个32位的字,只有高位被设置,那么它将只通过一次循环。
使用64位指令对14、24或32位字中设置的位进行计数
1 unsigned int v; // count the number of bits set in v
2 unsigned int c; // c accumulates the total bits set in v
3
4 // option 1, for at most 14-bit values in v:
5 c = (v * 0x200040008001ULL & 0x111111111111111ULL) % 0xf;
6
7 // option 2, for at most 24-bit values in v:
8 c = ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
9 c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL)
10 % 0x1f;
11
12 // option 3, for at most 32-bit values in v:
13 c = ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
14 c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) %
15 0x1f;
16 c += ((v >> 24) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
这种方法需要64位的CPU和快速的模数划分来提高效率。第一个选项只需要3次操作;第二个选项需要10次;第三个选项需要15次。
计数位设置,并行
1 unsigned int v; // count bits set in this (32-bit value)
2 unsigned int c; // store the total here
3 static const int S[] = {1, 2, 4, 8, 16}; // Magic Binary Numbers
4 static const int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF, 0x0000FFFF};
5
6 c = v - ((v >> 1) & B[0]);
7 c = ((c >> S[1]) & B[1]) + (c & B[1]);
8 c = ((c >> S[2]) + c) & B[2];
9 c = ((c >> S[3]) + c) & B[3];
10 c = ((c >> S[4]) + c) & B[4];
B数组,用二进制表示,是:
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位,那么我们需要数组S和B的元素长为ceil(lg(k)),我们必须计算与S或B长相同的c的表达式数量。对于32位的v,要用到16次运算。
32位整数v的最佳计数方法是以下:
v = v - ((v >> 1) & 0x55555555); // reuse input as temporary v = (v & 0x33333333) + ((v >> 2) & 0x33333333); // temp c = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24; // count
最好的位数计算方法只需要12次操作,这与查找表法相同,但避免了表的内存和潜在的缓存遗漏。它是上述纯并行方法和前面使用乘法的方法的混合体(在用64位指令计数位的部分),不过它没有使用64位指令。对字节中设置的位的计数是并行完成的,通过乘以0x1010101并右移24位,计算出字节中设置的位的总和。
最好的位数计算方法对位宽不超过128的整数(参数化为T型)的概括是这样的:
v = v - ((v >> 1) & (T)~(T)0/3); // temp v = (v & (T)~(T)0/15*3) + ((v >> 2) & (T)~(T)0/15*3); // temp v = (v + (v >> 4)) & (T)~(T)0/255*15; // temp c = (T)(v * ((T)~(T)0/255)) >> (sizeof(T) - 1) * CHAR_BIT; // count
从最重要的位到给定位置的计数位设置(等级)
下面查找一个位的等级,即返回从最重要的位到给定位置的位的位数之和。
1 uint64_t v; // Compute the rank (bits set) in v from the MSB to pos.
2 unsigned int pos; // Bit position to count bits upto.
3 uint64_t r; // Resulting rank of bit at pos goes here.
4
5 // Shift out bits after given position.
6 r = v >> (sizeof(v) * CHAR_BIT - pos);
7 // Count set bits in parallel.
8 // r = (r & 0x5555...) + ((r >> 1) & 0x5555...);
9 r = r - ((r >> 1) & ~0UL/3);
10 // r = (r & 0x3333...) + ((r >> 2) & 0x3333...);
11 r = (r & ~0UL/5) + ((r >> 2) & ~0UL/5);
12 // r = (r & 0x0f0f...) + ((r >> 4) & 0x0f0f...);
13 r = (r + (r >> 4)) & ~0UL/17;
14 // r = r % 255;
15 r = (r * (~0UL/255)) >> ((sizeof(v) - 1) * CHAR_BIT);
选择给定的位数(等级)的位数位置(从最重要的位)
下面的64位代码在从左数时选择第r 1位的位置。换句话说,如果我们从最有意义的位开始,然后向右进行,计数设置为1的位数,直到达到所需的位数r,那么我们停止的位置将被返回。如果要求的等级超过了设置的位数,则返回64。该代码可以修改为32位或从右开始计数。
1 uint64_t v; // Input value to find position with rank r.
2 unsigned int r; // Input: bit's desired rank [1-64].
3 unsigned int s; // Output: Resulting position of bit with rank r [1-64]
4 uint64_t a, b, c, d; // Intermediate temporaries for bit count.
5 unsigned int t; // Bit count temporary.
6
7 // Do a normal parallel bit count for a 64-bit integer,
8 // but store all intermediate steps.
9 // a = (v & 0x5555...) + ((v >> 1) & 0x5555...);
10 a = v - ((v >> 1) & ~0UL/3);
11 // b = (a & 0x3333...) + ((a >> 2) & 0x3333...);
12 b = (a & ~0UL/5) + ((a >> 2) & ~0UL/5);
13 // c = (b & 0x0f0f...) + ((b >> 4) & 0x0f0f...);
14 c = (b + (b >> 4)) & ~0UL/0x11;
15 // d = (c & 0x00ff...) + ((c >> 8) & 0x00ff...);
16 d = (c + (c >> 8)) & ~0UL/0x101;
17 t = (d >> 32) + (d >> 48);
18 // Now do branchless select!
19 s = 64;
20 // if (r > t) {s -= 32; r -= t;}
21 s -= ((t - r) & 256) >> 3; r -= (t & ((t - r) >> 8));
22 t = (d >> (s - 16)) & 0xff;
23 // if (r > t) {s -= 16; r -= t;}
24 s -= ((t - r) & 256) >> 4; r -= (t & ((t - r) >> 8));
25 t = (c >> (s - 8)) & 0xf;
26 // if (r > t) {s -= 8; r -= t;}
27 s -= ((t - r) & 256) >> 5; r -= (t & ((t - r) >> 8));
28 t = (b >> (s - 4)) & 0x7;
29 // if (r > t) {s -= 4; r -= t;}
30 s -= ((t - r) & 256) >> 6; r -= (t & ((t - r) >> 8));
31 t = (a >> (s - 2)) & 0x3;
32 // if (r > t) {s -= 2; r -= t;}
33 s -= ((t - r) & 256) >> 7; r -= (t & ((t - r) >> 8));
34 t = (v >> (s - 1)) & 0x1;
35 // if (r > t) s--;
36 s -= ((t - r) & 256) >> 8;
37 s = 65 - s;
如果你的目标CPU上的分支速度很快,可以考虑取消对if-statements的注释,并对它们后面的行进行注释。
以天真的方式计算奇偶性
unsigned int v; // word value to compute the parity of
bool parity = false; // parity will be the parity of v
while (v)
{
parity = !parity;
v = v & (v - 1);
}
上面的代码使用了类似上面Brian Kernigan的位数计算的方法。它所需要的时间与设置的位数成正比。
通过查找表计算奇偶性
1 static const bool ParityTable256[256] =
2 {
3 # define P2(n) n, n^1, n^1, n
4 # define P4(n) P2(n), P2(n^1), P2(n^1), P2(n)
5 # define P6(n) P4(n), P4(n^1), P4(n^1), P4(n)
6 P6(0), P6(1), P6(1), P6(0)
7 };
8
9 unsigned char b; // byte value to compute the parity of
10 bool parity = ParityTable256[b];
11
12 // OR, for 32-bit words:
13 unsigned int v;
14 v ^= v >> 16;
15 v ^= v >> 8;
16 bool parity = ParityTable256[v & 0xff];
17
18 // Variation:
19 unsigned char * p = (unsigned char *) &v;
20 parity = ParityTable256[p[0] ^ p[1] ^ p[2] ^ p[3]];
使用64位乘法和模数除法计算字节的奇偶性
unsigned char b; // byte value to compute the parity of
bool parity =
(((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;
上面的方法大约需要4次操作,但只对字节有效。
用乘法计算字的奇偶性
下面的方法用乘法只用8次操作就计算出32位值的奇偶性。
1 unsigned int v; // 32-bit word
2 v ^= v >> 1;
3 v ^= v >> 2;
4 v = (v & 0x11111111U) * 0x11111111U;
5 return (v >> 28) & 1;
同样是64位,8次操作还是足够的。
unsigned long long v; // 64-bit word
v ^= v >> 1;
v ^= v >> 2;
v = (v & 0x1111111111111111UL) * 0x1111111111111111UL;
return (v >> 60) & 1;
并行计算奇偶性
unsigned int v; // word value to compute the parity of
v ^= v >> 16;
v ^= v >> 8;
v ^= v >> 4;
v &= 0xf;
return (0x6996 >> v) & 1;
上面的方法需要大约9次操作,并且适用于32位的字。通过删除 "unsigned int v; "后面的两行代码,可以优化为只对字节进行5次操作。该方法首先将32位值的8个字头一起移位和XOR,将结果留在v的最低字头中。接下来,二进制数0110 1001 1001 0110(十六进制中的0x6996)被v的最低字头所代表的值向右移位,这个数字就像一个微型的16位奇偶表,以v中的低4位为索引。
用减法和加法交换数值
#define SWAP(a, b) ((&(a) == &(b)) || \
(((a) -= (b)), ((b) += (a)), ((a) = (b) - (a))))
这将不使用临时变量来交换a和b的值。当你知道这不可能发生时,可以省略a和b在内存中的同一位置的初始检查。(编译器可能会作为一种优化而省略它)。(编译器可能会作为一种优化而省略它。)如果你启用了溢出异常,那么传递无符号值,这样就不会抛出异常。下面的XOR方法在某些机器上可能会稍微快一些。不要对浮点数使用这种方法(除非你对它们的原始整数表示进行操作)
用XOR交换值
#define SWAP(a, b) (((a) ^= (b)), ((b) ^= (a)), ((a) ^= (b)))
这是一个古老的技巧,可以交换变量a和b的值,而不需要为临时变量使用额外的空间。
用XOR交换单个位
unsigned int i, j; // positions of bit sequences to swap
unsigned int n; // number of consecutive bits in each sequence
unsigned int b; // bits to swap reside in b
unsigned int r; // bit-swapped result goes here
unsigned int x = ((b >> i) ^ (b >> j)) & ((1U << n) - 1); // XOR temporary
r = b ^ ((x << i) | (x << j));
作为一个交换比特范围的例子,假设我们有b=00101111(用二进制表示),我们想把从i=1(右起第二位)开始的n=3个连续比特与从j=5开始的3个连续比特交换;结果将是r=11100011(二进制)。
这种交换方法类似于通用的XOR交换技巧,但旨在对单个位进行操作。 变量x存储了我们要交换的比特值对的XOR结果,然后比特被设置为自己与x的XOR结果。当然,如果序列重叠,结果是未定义的。
用明显的方法反转位
1 unsigned int v; // input bits to be reversed
2 unsigned int r = v; // r will be reversed bits of v; first get LSB of v
3 int s = sizeof(v) * CHAR_BIT - 1; // extra shift needed at end
4
5 for (v >>= 1; v; v >>= 1)
6 {
7 r <<= 1;
8 r |= v & 1;
9 s--;
10 }
11 r <<= s; // shift when v's highest bits are zero
通过查找表反转字中的位
1 static const unsigned char BitReverseTable256[256] =
2 {
3 # define R2(n) n, n + 2*64, n + 1*64, n + 3*64
4 # define R4(n) R2(n), R2(n + 2*16), R2(n + 1*16), R2(n + 3*16)
5 # define R6(n) R4(n), R4(n + 2*4 ), R4(n + 1*4 ), R4(n + 3*4 )
6 R6(0), R6(2), R6(1), R6(3)
7 };
8
9 unsigned int v; // reverse 32-bit value, 8 bits at time
10 unsigned int c; // c will get v reversed
11
12 // Option 1:
13 c = (BitReverseTable256[v & 0xff] << 24) |
14 (BitReverseTable256[(v >> 8) & 0xff] << 16) |
15 (BitReverseTable256[(v >> 16) & 0xff] << 8) |
16 (BitReverseTable256[(v >> 24) & 0xff]);
17
18 // Option 2:
19 unsigned char * p = (unsigned char *) &v;
20 unsigned char * q = (unsigned char *) &c;
21 q[3] = BitReverseTable256[p[0]];
22 q[2] = BitReverseTable256[p[1]];
23 q[1] = BitReverseTable256[p[2]];
24 q[0] = BitReverseTable256[p[3]];
第一种方法大约需要17次操作,第二种方法大约需要12次操作,假设你的CPU可以轻松加载和存储字节。
用3个操作(64位乘法和模数除法)反转一个字节中的位
unsigned char b; // reverse this (8-bit) byte
b = (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;
乘法运算为8位字节模式创建5个独立的副本,以扇出成一个64位的值。AND操作选择相对于每个10位比特组的正确(反向)位置的比特。乘法和AND操作将原始字节中的位复制出来,使它们各自只出现在10位组中的一个位中。原始字节中的比特的反转位置与它们在任何10位组中的相对位置一致。最后一步,涉及模数除以2^10 - 1,效果是将64位值中的每一组10位(从0-9,10-19,20-29,...的位置)合并在一起。它们并不重叠,所以模数除法所依据的加法步骤就像或运算一样。
用4次操作将一个字节中的位数反转(64位乘,不除)
unsigned char b; // reverse this byte
b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;
下面显示了布尔变量a、b、c、d、e、f、g和h的位值流程,这些变量组成了一个8位字节。请注意第一个乘法是如何将比特模式扇出到多个副本,而最后一个乘法是如何将它们合并到右边第五个字节中。
1 abcd efgh (-> hgfe dcba)
2 * 1000 0000 0010 0000 0000 1000 0000 0010 (0x80200802)
3 -------------------------------------------------------------------------------------------------
4 0abc defg h00a bcde fgh0 0abc defg h00a bcde fgh0
5 & 0000 1000 1000 0100 0100 0010 0010 0001 0001 0000 (0x0884422110)
6 -------------------------------------------------------------------------------------------------
7 0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
8 * 0000 0001 0000 0001 0000 0001 0000 0001 0000 0001 (0x0101010101)
9 -------------------------------------------------------------------------------------------------
10 0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
11 0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
12 0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
13 0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
14 0000 d000 h000 0c00 0g00 00b0 00f0 000a 000e 0000
15 -------------------------------------------------------------------------------------------------
16 0000 d000 h000 dc00 hg00 dcb0 hgf0 dcba hgfe dcba hgfe 0cba 0gfe 00ba 00fe 000a 000e 0000
17 >> 32
18 -------------------------------------------------------------------------------------------------
19 0000 d000 h000 dc00 hg00 dcb0 hgf0 dcba hgfe dcba
20 & 1111 1111
21 -------------------------------------------------------------------------------------------------
22 hgfe dcba
请注意,在某些处理器上,最后两步可以结合起来,因为寄存器可以作为字节访问;只需乘法,使一个寄存器存储结果的上32位,取低字节。因此,可能只需要6次操作。
用7个操作反转一个字节中的位(没有64位)
1 unsigned int v; // 32-bit word to reverse bit order
2
3 // swap odd and even bits
4 v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
5 // swap consecutive pairs
6 v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
7 // swap nibbles ...
8 v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
9 // swap bytes
10 v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
11 // swap 2-byte long pairs
12 v = ( v >> 16 ) | ( v << 16);
下面的变化也是O(lg(N)),但是它需要更多的操作来反转v,它的优点是通过在飞行中计算常数,占用较少的内存。
unsigned int s = sizeof(v) * CHAR_BIT; // bit size; must be power of 2
unsigned int mask = ~0;
while ((s >>= 1) > 0)
{
mask ^= (mask << s);
v = ((v >> s) & mask) | ((v << s) & ~mask);
}
以上这些方法最适合N大的情况。如果你在64位ints(或更大)的情况下使用上面的方法,那么你需要增加更多的行数(按照模式),否则只有低32位会被反转,结果会在低32位。
计算模数除以1<<s,不需要除法运算符
1 const unsigned int n; // numerator
2 const unsigned int s;
3 const unsigned int d = 1U << s; // So d will be one of: 1, 2, 4, 8, 16, 32, ...
4 unsigned int m; // m will be n % d
5 m = n & (d - 1);
大多数程序员很早就学会了这一招,但为了完整起见,还是把它包括进去了。
计算模数除以(1<<s)-1,不含除法运算符
1 unsigned int n; // numerator
2 const unsigned int s; // s > 0
3 const unsigned int d = (1 << s) - 1; // so d is either 1, 3, 7, 15, 31, ...).
4 unsigned int m; // n % d goes here.
5
6 for (m = n; n > d; n = m)
7 {
8 for (m = 0; n; n >>= s)
9 {
10 m += n & d;
11 }
12 }
13 // Now m is a value from 0 to d, but since with modulus division
14 // we want m to be 0 when it is d.
15 m = m == d ? 0 : m;
这种模数除以小于2次方的整数的方法最多需要5+(4+5*ceil(N/s)) * ceil(lg(N / s))操作,其中N是分子中的位数。换句话说,它最多需要O(N * lg(N))时间。
在没有除法运算符的情况下,并行计算(1<<s)-1的模数除法
1 // The following is for a word size of 32 bits!
2
3 static const unsigned int M[] =
4 {
5 0x00000000, 0x55555555, 0x33333333, 0xc71c71c7,
6 0x0f0f0f0f, 0xc1f07c1f, 0x3f03f03f, 0xf01fc07f,
7 0x00ff00ff, 0x07fc01ff, 0x3ff003ff, 0xffc007ff,
8 0xff000fff, 0xfc001fff, 0xf0003fff, 0xc0007fff,
9 0x0000ffff, 0x0001ffff, 0x0003ffff, 0x0007ffff,
10 0x000fffff, 0x001fffff, 0x003fffff, 0x007fffff,
11 0x00ffffff, 0x01ffffff, 0x03ffffff, 0x07ffffff,
12 0x0fffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff
13 };
14
15 static const unsigned int Q[][6] =
16 {
17 { 0, 0, 0, 0, 0, 0}, {16, 8, 4, 2, 1, 1}, {16, 8, 4, 2, 2, 2},
18 {15, 6, 3, 3, 3, 3}, {16, 8, 4, 4, 4, 4}, {15, 5, 5, 5, 5, 5},
19 {12, 6, 6, 6 , 6, 6}, {14, 7, 7, 7, 7, 7}, {16, 8, 8, 8, 8, 8},
20 { 9, 9, 9, 9, 9, 9}, {10, 10, 10, 10, 10, 10}, {11, 11, 11, 11, 11, 11},
21 {12, 12, 12, 12, 12, 12}, {13, 13, 13, 13, 13, 13}, {14, 14, 14, 14, 14, 14},
22 {15, 15, 15, 15, 15, 15}, {16, 16, 16, 16, 16, 16}, {17, 17, 17, 17, 17, 17},
23 {18, 18, 18, 18, 18, 18}, {19, 19, 19, 19, 19, 19}, {20, 20, 20, 20, 20, 20},
24 {21, 21, 21, 21, 21, 21}, {22, 22, 22, 22, 22, 22}, {23, 23, 23, 23, 23, 23},
25 {24, 24, 24, 24, 24, 24}, {25, 25, 25, 25, 25, 25}, {26, 26, 26, 26, 26, 26},
26 {27, 27, 27, 27, 27, 27}, {28, 28, 28, 28, 28, 28}, {29, 29, 29, 29, 29, 29},
27 {30, 30, 30, 30, 30, 30}, {31, 31, 31, 31, 31, 31}
28 };
29
30 static const unsigned int R[][6] =
31 {
32 {0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000},
33 {0x0000ffff, 0x000000ff, 0x0000000f, 0x00000003, 0x00000001, 0x00000001},
34 {0x0000ffff, 0x000000ff, 0x0000000f, 0x00000003, 0x00000003, 0x00000003},
35 {0x00007fff, 0x0000003f, 0x00000007, 0x00000007, 0x00000007, 0x00000007},
36 {0x0000ffff, 0x000000ff, 0x0000000f, 0x0000000f, 0x0000000f, 0x0000000f},
37 {0x00007fff, 0x0000001f, 0x0000001f, 0x0000001f, 0x0000001f, 0x0000001f},
38 {0x00000fff, 0x0000003f, 0x0000003f, 0x0000003f, 0x0000003f, 0x0000003f},
39 {0x00003fff, 0x0000007f, 0x0000007f, 0x0000007f, 0x0000007f, 0x0000007f},
40 {0x0000ffff, 0x000000ff, 0x000000ff, 0x000000ff, 0x000000ff, 0x000000ff},
41 {0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff, 0x000001ff},
42 {0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff, 0x000003ff},
43 {0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff, 0x000007ff},
44 {0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff, 0x00000fff},
45 {0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff, 0x00001fff},
46 {0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff, 0x00003fff},
47 {0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff, 0x00007fff},
48 {0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff, 0x0000ffff},
49 {0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff, 0x0001ffff},
50 {0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff, 0x0003ffff},
51 {0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff, 0x0007ffff},
52 {0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff, 0x000fffff},
53 {0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff, 0x001fffff},
54 {0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff, 0x003fffff},
55 {0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff, 0x007fffff},
56 {0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff, 0x00ffffff},
57 {0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff, 0x01ffffff},
58 {0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff},
59 {0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff, 0x07ffffff},
60 {0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff, 0x0fffffff},
61 {0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff},
62 {0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff, 0x3fffffff},
63 {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff}
64 };
65
66 unsigned int n; // numerator
67 const unsigned int s; // s > 0
68 const unsigned int d = (1 << s) - 1; // so d is either 1, 3, 7, 15, 31, ...).
69 unsigned int m; // n % d goes here.
70
71 m = (n & M[s]) + ((n >> s) & M[s]);
72
73 for (const unsigned int * q = &Q[s][0], * r = &R[s][0]; m > d; q++, r++)
74 {
75 m = (m >> *q) + (m & *r);
76 }
77 m = m == d ? 0 : m; // OR, less portably: m = m & -((signed)(m - d) >> s);
这种寻找小于2次幂的整数模数除法的方法最多需要O(lg(N))时间,其中N是分子中的位数(32位,上面的代码)。操作次数最多为12 + 9 * ceil(lg(N))。如果你在编译时知道分母,就可以去掉表格;只需提取几个相关的条目并展开循环。它可以很容易地扩展到更多的位。
它通过将基数(1<<s)中的值平行相加来计算结果。首先,每隔一个基数(1<<s)的值都要加到前一个值上。想象一下,结果是写在一张纸上的。将纸切成两半,这样每张纸上就有一半的值。将这些值对齐,并将它们加到一张新的纸上。重复将这张纸切成两半(将是前一张纸的四分之一大小),然后求和,直到不能再切下去。在进行lg(N/s/2)次切割后,我们就不再切割了,只需继续将数值相加,并将结果放到一张新的纸上,就像之前一样,同时至少有两个s位数值。
用O(N)运算(显而易见的方法)求一个整数的对数基数2与MSB N的集合
unsigned int v; // 32-bit word to find the log base 2 of
unsigned int r = 0; // r will be lg(v)
while (v >>= 1) // unroll for more speed...
{
r++;
}
一个整数的对数基数2与最高位集(或最重要位集,MSB)的位置相同。以下的对数基数2方法比这个方法更快。
求一个64位IEEE浮点数的整数对数基数2
1 int v; // 32-bit integer to find the log base 2 of
2 int r; // result of log_2(v) goes here
3 union { unsigned int u[2]; double d; } t; // temp
4
5 t.u[__FLOAT_WORD_ORDER==LITTLE_ENDIAN] = 0x43300000;
6 t.u[__FLOAT_WORD_ORDER!=LITTLE_ENDIAN] = v;
7 t.d -= 4503599627370496.0;
8 r = (t.u[__FLOAT_WORD_ORDER==LITTLE_ENDIAN] >> 20) - 0x3FF;
上面的代码将一个64位(IEEE-754浮点)的双倍数与一个32位的整数(没有垫位)一起加载,在指数设置为252的同时,将整数存储在尾数中。从这个新造的双数中减去252(用双数表示),将所得指数设置为输入值v的对数基数2,剩下的就是将指数位移到位置上(向右移20位),并减去偏置位0x3FF(即1023十进制)。这种技术只需要5次操作,但是很多CPU在处理双数时速度很慢,必须要照顾到架构的endianess。
用查找表求一个整数的对数基数2
1 static const char LogTable256[256] =
2 {
3 #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
4 -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
5 LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
6 LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
7 };
8
9 unsigned int v; // 32-bit word to find the log of
10 unsigned r; // r will be lg(v)
11 register unsigned int t, tt; // temporaries
12
13 if (tt = v >> 16)
14 {
15 r = (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
16 }
17 else
18 {
19 r = (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
20 }
查表法只需要大约7次操作就可以找到一个32位值的对数。如果扩展到64位的量,大概需要9次操作。通过使用四个表,并将可能的加法纳入到每个表中,可以减去另一个操作。使用int表元素可能更快,这取决于你的架构。
上面的代码被调整为均匀分布的输出值。如果你的输入是均匀分布在所有32位的值上,那么可以考虑使用下面的代码:
1 if (tt = v >> 24)
2 {
3 r = 24 + LogTable256[tt];
4 }
5 else if (tt = v >> 16)
6 {
7 r = 16 + LogTable256[tt];
8 }
9 else if (tt = v >> 8)
10 {
11 r = 8 + LogTable256[tt];
12 }
13 else
14 {
15 r = LogTable256[v];
16 }
要通过算法初步生成日志表:
1 LogTable256[0] = LogTable256[1] = 0;
2 for (int i = 2; i < 256; i++)
3 {
4 LogTable256[i] = 1 + LogTable256[i / 2];
5 }
6 LogTable256[0] = -1; // if you want log(0) to return -1
在O(lg(N))运算中求一个N位整数的对数基数2
1 unsigned int v; // 32-bit value to find the log2 of
2 const unsigned int b[] = {0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000};
3 const unsigned int S[] = {1, 2, 4, 8, 16};
4 int i;
5
6 register unsigned int r = 0; // result of log2(v) will go here
7 for (i = 4; i >= 0; i--) // unroll for speed...
8 {
9 if (v & b[i])
10 {
11 v >>= S[i];
12 r |= S[i];
13 }
14 }
15
16
17 // OR (IF YOUR CPU BRANCHES SLOWLY):
18
19 unsigned int v; // 32-bit value to find the log2 of
20 register unsigned int r; // result of log2(v) will go here
21 register unsigned int shift;
22
23 r = (v > 0xFFFF) << 4; v >>= r;
24 shift = (v > 0xFF ) << 3; v >>= shift; r |= shift;
25 shift = (v > 0xF ) << 2; v >>= shift; r |= shift;
26 shift = (v > 0x3 ) << 1; v >>= shift; r |= shift;
27 r |= (v >> 1);
28
29
30 // OR (IF YOU KNOW v IS A POWER OF 2):
31
32 unsigned int v; // 32-bit value to find the log2 of
33 static const unsigned int b[] = {0xAAAAAAAA, 0xCCCCCCCC, 0xF0F0F0F0,
34 0xFF00FF00, 0xFFFF0000};
35 register unsigned int r = (v & b[0]) != 0;
36 for (i = 4; i > 0; i--) // unroll for speed...
37 {
38 r |= ((v & b[i]) != 0) << i;
39 }
当然,如果要将代码扩展到查找33位到64位数的对数,我们会将另一个元素0xFFFFFFFF00000000追加到b上,将32追加到S上,然后从5循环到0,这种方法比早期的查表版本要慢得多,但如果你不想要大表,或者你的架构对内存的访问速度很慢,它是一个不错的选择。第二种变化涉及的操作略多,但在分支成本高的机器上可能更快(如PowerPC)。
用O(lg(N))运算求一个N位整数的对数基数2,并进行乘法和查找
1 uint32_t v; // find the log base 2 of 32-bit v
2 int r; // result goes here
3
4 static const int MultiplyDeBruijnBitPosition[32] =
5 {
6 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
7 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
8 };
9
10 v |= v >> 1; // first round down to one less than a power of 2
11 v |= v >> 2;
12 v |= v >> 4;
13 v |= v >> 8;
14 v |= v >> 16;
15
16 r = MultiplyDeBruijnBitPosition[(uint32_t)(v * 0x07C4ACDDU) >> 27];
上面的代码用小表查找和乘法计算一个32位整数的对数基数2。它只需要13次操作,而前一种方法需要20次操作(最多)。纯粹基于表的方法需要最少的操作,但这在表的大小和速度之间提供了一个合理的折中。
如果你知道v是2的幂,那么你只需要以下几点:
1 static const int MultiplyDeBruijnBitPosition2[32] =
2 {
3 0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
4 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
5 };
6 r = MultiplyDeBruijnBitPosition2[(uint32_t)(v * 0x077CB531U) >> 27];
寻找整数对数基数10的整数
1 unsigned int v; // non-zero 32-bit integer value to compute the log base 10 of
2 int r; // result goes here
3 int t; // temporary
4
5 static unsigned int const PowersOf10[] =
6 {1, 10, 100, 1000, 10000, 100000,
7 1000000, 10000000, 100000000, 1000000000};
8
9 t = (IntegerLogBase2(v) + 1) * 1233 >> 12; // (use a lg2 method from above)
10 r = t - (v < PowersOf10[t]);
整数对数基数10的计算方法是先用上面求对数基数2的技术之一。通过log10(v)=log2(v)/log2(10)的关系,我们需要将其乘以1/log2(10),大约是1233/4096,也就是1233后右移12。因为IntegerLogBase2四舍五入,所以需要加一。最后,由于值t只是一个近似值,可能会偏离1,所以要减去v<PowersOf10[t]的结果,才能找到准确的值。
这个方法比IntegerLogBase2多了6次操作。可以通过修改上面的对数基数2表查找方法来加快速度(在内存访问速度快的机器上),使条目保持为t计算的内容(即预添加、-mulitply和-shift)。这样做总共只需要9次操作就可以找到对数基10,假设使用4个表(v的每个字节一个表)。
寻找整数对数基数10的明显方法
unsigned int v; // non-zero 32-bit integer value to compute the log base 10 of
int r; // result goes here
r = (v >= 1000000000) ? 9 : (v >= 100000000) ? 8 : (v >= 10000000) ? 7 :
(v >= 1000000) ? 6 : (v >= 100000) ? 5 : (v >= 10000) ? 4 :
(v >= 1000) ? 3 : (v >= 100) ? 2 : (v >= 10) ? 1 : 0;
当输入均匀分布在32位值上时,这种方法很好用,因为76%的输入被第一次比较捕获,21%的输入被第二次比较捕获,2%的输入被第三次比较捕获,以此类推(每次比较将剩余的输入砍掉90%)。因此,平均需要不到2.6次操作。
查找32位IEEE浮点数的整数对数基数2
1 const float v; // find int(log2(v)), where v > 0.0 && finite(v) && isnormal(v)
2 int c; // 32-bit int c gets the result;
3
4 c = *(const int *) &v; // OR, for portability: memcpy(&c, &v, sizeof c);
5 c = (c >> 23) - 127;
6 The above is fast, but IEEE 754-compliant architectures utilize subnormal (also called denormal) floating point numbers. These have the exponent bits set to zero (signifying pow(2,-127)), and the mantissa is not normalized, so it contains leading zeros and thus the log2 must be computed from the mantissa. To accomodate for subnormal numbers, use the following:
7 const float v; // find int(log2(v)), where v > 0.0 && finite(v)
8 int c; // 32-bit int c gets the result;
9 int x = *(const int *) &v; // OR, for portability: memcpy(&x, &v, sizeof x);
10
11 c = x >> 23;
12
13 if (c)
14 {
15 c -= 127;
16 }
17 else
18 { // subnormal, so recompute using mantissa: c = intlog2(x) - 149;
19 register unsigned int t; // temporary
20 // Note that LogTable256 was defined earlier
21 if (t = x >> 16)
22 {
23 c = LogTable256[t] - 133;
24 }
25 else
26 {
27 c = (t = x >> 8) ? LogTable256[t] - 141 : LogTable256[x] - 149;
28 }
29 }
求一个32位IEEE浮点数的pow(2, r)根的整数对数基2(对于无符号整数r)
1 const int r;
2 const float v; // find int(log2(pow((double) v, 1. / pow(2, r)))),
3 // where isnormal(v) and v > 0
4 int c; // 32-bit int c gets the result;
5
6 c = *(const int *) &v; // OR, for portability: memcpy(&c, &v, sizeof c);
7 c = ((((c - 0x3f800000) >> r) + 0x3f800000) >> 23) - 127;
所以,例如,如果r为0,我们有c=int(log2((double)v))。如果r是1,那么我们有c=int(log2(sqrt((double)v))。如果r是2,那么我们有c = int(log2(pow((double) v, 1./4)))
以线性方式计算右侧的连续零位(尾部)
1 unsigned int v; // input to count trailing zero bits
2 int c; // output: c will count v's trailing zero bits,
3 // so if v is 1101000 (base 2), then c will be 3
4 if (v)
5 {
6 v = (v ^ (v - 1)) >> 1; // Set v's trailing 0s to 1s and zero rest
7 for (c = 0; v; c++)
8 {
9 v >>= 1;
10 }
11 }
12 else
13 {
14 c = CHAR_BIT * sizeof(v);
15 }
(均匀分布的)随机二进制数的平均尾部零位数是1,所以这个O(尾部零位)的解法与下面更快的方法相比并不差。
并行计算右边的连续零位(尾部)
1 unsigned int v; // 32-bit word input to count zero bits on right
2 unsigned int c = 32; // c will be the number of zero bits on the right
3 v &= -signed(v);
4 if (v) c--;
5 if (v & 0x0000FFFF) c -= 16;
6 if (v & 0x00FF00FF) c -= 8;
7 if (v & 0x0F0F0F0F) c -= 4;
8 if (v & 0x33333333) c -= 2;
9 if (v & 0x55555555) c -= 1;
在这里,我们基本上是在并行地做着和找对数基数2一样的操作,但我们先隔离出最低的1位,然后从最大的c开始递减。大致上,对于N位字,操作次数最多为3*lg(N)+4。
通过二进制搜索计算右边的连续零位(尾部)
1 unsigned int v; // 32-bit word input to count zero bits on right
2 unsigned int c; // c will be the number of zero bits on the right,
3 // so if v is 1101000 (base 2), then c will be 3
4 // NOTE: if 0 == v, then c = 31.
5 if (v & 0x1)
6 {
7 // special case for odd v (assumed to happen half of the time)
8 c = 0;
9 }
10 else
11 {
12 c = 1;
13 if ((v & 0xffff) == 0)
14 {
15 v >>= 16;
16 c += 16;
17 }
18 if ((v & 0xff) == 0)
19 {
20 v >>= 8;
21 c += 8;
22 }
23 if ((v & 0xf) == 0)
24 {
25 v >>= 4;
26 c += 4;
27 }
28 if ((v & 0x3) == 0)
29 {
30 v >>= 2;
31 c += 2;
32 }
33 c -= v & 0x1;
34 }
上面的代码与前面的方法类似,但它通过累加c的方式计算尾部零的数量,其方式类似于二进制搜索。在第一步中,它检查v的底层16位是否为0,如果是,则将v右移16位,并将16加到c中,这样v中要考虑的位数减少一半。随后的每一个条件步骤同样将位数减半,直到只有1位。这个方法比上一个方法快(大约33%),因为if语句的主体被执行的次数较少。
通过投掷到浮点数来计算右边的连续零位(尾部)
1 unsigned int v; // find the number of trailing zeros in v 2 int r; // the result goes here 3 float f = (float)(v & -v); // cast the least significant bit in v to a float 4 r = (*(uint32_t *)&f >> 23) - 0x7f;
虽然这只需要大约6次操作,但在某些机器上,将整数转换为浮点数的时间可能很高。将32位IEEE浮点表示法的指数下移,再减去偏置,得到v中设置的最不重要的1位的位置,如果v为零,那么结果就是-127。
用模数除法和查找法计算右边的连续零位(尾数)
1 unsigned int v; // find the number of trailing zeros in v
2 int r; // put the result in r
3 static const int Mod37BitPosition[] = // map a bit value mod 37 to its position
4 {
5 32, 0, 1, 26, 2, 23, 27, 0, 3, 16, 24, 30, 28, 11, 0, 13, 4,
6 7, 17, 0, 25, 22, 31, 15, 29, 10, 12, 6, 0, 21, 14, 9, 5,
7 20, 8, 19, 18
8 };
9 r = Mod37BitPosition[(-v & v) % 37];
上面的代码找到了右边尾部的零的数量,所以二进制0100会产生2。它利用了前32位位置值与37相对质数的事实,所以用37执行模数除法,每个位置值都有一个从0到36的唯一数字。然后,这些数字可以使用一个小的查找表映射到0的数量。它只使用了4个操作,然而索引到表中并执行模数除法可能使它不适合某些情况。我独立想出了这个办法,然后搜索表值的子序列,发现它是Reiser较早发明的,根据Hacker's Delight。
乘法和查找来计算右边的连续零位(尾部)
1 unsigned int v; // find the number of trailing zeros in 32-bit v
2 int r; // result goes here
3 static const int MultiplyDeBruijnBitPosition[32] =
4 {
5 0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
6 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
7 };
8 r = MultiplyDeBruijnBitPosition[((uint32_t)((v & -v) * 0x077CB531U)) >> 27];
用浮投法将其四舍五入到2的下一个最高倍数
1 unsigned int const v; // Round this 32-bit value to the next highest power of 2
2 unsigned int r; // Put the result here. (So v=3 -> r=4; v=8 -> r=8)
3
4 if (v > 1)
5 {
6 float f = (float)v;
7 unsigned int const t = 1U << ((*(unsigned int *)&f >> 23) - 0x7f);
8 r = t << (t < v);
9 }
10 else
11 {
12 r = 1;
13 }
14 The code above uses 8 operations, but works on all v <= (1<<31).
15 Quick and dirty version, for domain of 1 < v < (1<<25):
16
17 float f = (float)(v - 1);
18 r = 1U << ((*(unsigned int*)(&f) >> 23) - 126);
虽然快速和肮脏的版本只使用了大约6个操作,但在Athlon™ XP 2100+ CPU上进行基准测试时,它比下面的技术(涉及12个操作)慢了大约三倍。不过,有些CPU的表现会更好。
四舍五入到2的下一个最高倍数
1 unsigned int v; // compute the next highest power of 2 of 32-bit v
2
3 v--;
4 v |= v >> 1;
5 v |= v >> 2;
6 v |= v >> 4;
7 v |= v >> 8;
8 v |= v >> 16;
9 v++;
在12次操作中,这个代码计算一个32位整数的2的次高次幂。结果可以用公式1U << (lg(v - 1) + 1)来表示。请注意,在v为0的边缘情况下,它返回的是0,这不是2的幂;如果重要的话,你可以附加表达式v += (v == 0)来补救。使用公式和使用查找表的log base 2方法会快2个操作,但在某些情况下,查找表并不适合,所以上述代码可能是最好的。(在Athlon™ XP 2100+上,我发现上述shift-left然后OR代码的速度与使用单条BSR汇编语言指令的速度一样快,该指令反向扫描以找到最高集位)。它的工作原理是将最高集位复制到所有的低位,然后加1,结果是将所有的低位都设置为0,最高集位以外的1位设置为1,如果原来的数字是2的幂,那么减法就会将其减少1,这样我们就可以四舍五入到相同的原值。
以明显的方式交错位
1 unsigned short x; // Interleave bits of x and y, so that all of the
2 unsigned short y; // bits of x are in the even positions and y in the odd;
3 unsigned int z = 0; // z gets the resulting Morton Number.
4
5 for (int i = 0; i < sizeof(x) * CHAR_BIT; i++) // unroll for more speed...
6 {
7 z |= (x & 1U << i) << i | (y & 1U << i) << (i + 1);
8 }
交错位(又名莫顿数)对于二维整数坐标的线性化非常有用,因此x和y被组合成一个单一的数字,可以很容易地进行比较,并且具有这样的特性:如果一个数字的x和y值很接近,那么它通常与另一个数字很接近。
通过表格查找交错位
1 static const unsigned short MortonTable256[256] =
2 {
3 0x0000, 0x0001, 0x0004, 0x0005, 0x0010, 0x0011, 0x0014, 0x0015,
4 0x0040, 0x0041, 0x0044, 0x0045, 0x0050, 0x0051, 0x0054, 0x0055,
5 0x0100, 0x0101, 0x0104, 0x0105, 0x0110, 0x0111, 0x0114, 0x0115,
6 0x0140, 0x0141, 0x0144, 0x0145, 0x0150, 0x0151, 0x0154, 0x0155,
7 0x0400, 0x0401, 0x0404, 0x0405, 0x0410, 0x0411, 0x0414, 0x0415,
8 0x0440, 0x0441, 0x0444, 0x0445, 0x0450, 0x0451, 0x0454, 0x0455,
9 0x0500, 0x0501, 0x0504, 0x0505, 0x0510, 0x0511, 0x0514, 0x0515,
10 0x0540, 0x0541, 0x0544, 0x0545, 0x0550, 0x0551, 0x0554, 0x0555,
11 0x1000, 0x1001, 0x1004, 0x1005, 0x1010, 0x1011, 0x1014, 0x1015,
12 0x1040, 0x1041, 0x1044, 0x1045, 0x1050, 0x1051, 0x1054, 0x1055,
13 0x1100, 0x1101, 0x1104, 0x1105, 0x1110, 0x1111, 0x1114, 0x1115,
14 0x1140, 0x1141, 0x1144, 0x1145, 0x1150, 0x1151, 0x1154, 0x1155,
15 0x1400, 0x1401, 0x1404, 0x1405, 0x1410, 0x1411, 0x1414, 0x1415,
16 0x1440, 0x1441, 0x1444, 0x1445, 0x1450, 0x1451, 0x1454, 0x1455,
17 0x1500, 0x1501, 0x1504, 0x1505, 0x1510, 0x1511, 0x1514, 0x1515,
18 0x1540, 0x1541, 0x1544, 0x1545, 0x1550, 0x1551, 0x1554, 0x1555,
19 0x4000, 0x4001, 0x4004, 0x4005, 0x4010, 0x4011, 0x4014, 0x4015,
20 0x4040, 0x4041, 0x4044, 0x4045, 0x4050, 0x4051, 0x4054, 0x4055,
21 0x4100, 0x4101, 0x4104, 0x4105, 0x4110, 0x4111, 0x4114, 0x4115,
22 0x4140, 0x4141, 0x4144, 0x4145, 0x4150, 0x4151, 0x4154, 0x4155,
23 0x4400, 0x4401, 0x4404, 0x4405, 0x4410, 0x4411, 0x4414, 0x4415,
24 0x4440, 0x4441, 0x4444, 0x4445, 0x4450, 0x4451, 0x4454, 0x4455,
25 0x4500, 0x4501, 0x4504, 0x4505, 0x4510, 0x4511, 0x4514, 0x4515,
26 0x4540, 0x4541, 0x4544, 0x4545, 0x4550, 0x4551, 0x4554, 0x4555,
27 0x5000, 0x5001, 0x5004, 0x5005, 0x5010, 0x5011, 0x5014, 0x5015,
28 0x5040, 0x5041, 0x5044, 0x5045, 0x5050, 0x5051, 0x5054, 0x5055,
29 0x5100, 0x5101, 0x5104, 0x5105, 0x5110, 0x5111, 0x5114, 0x5115,
30 0x5140, 0x5141, 0x5144, 0x5145, 0x5150, 0x5151, 0x5154, 0x5155,
31 0x5400, 0x5401, 0x5404, 0x5405, 0x5410, 0x5411, 0x5414, 0x5415,
32 0x5440, 0x5441, 0x5444, 0x5445, 0x5450, 0x5451, 0x5454, 0x5455,
33 0x5500, 0x5501, 0x5504, 0x5505, 0x5510, 0x5511, 0x5514, 0x5515,
34 0x5540, 0x5541, 0x5544, 0x5545, 0x5550, 0x5551, 0x5554, 0x5555
35 };
36
37 unsigned short x; // Interleave bits of x and y, so that all of the
38 unsigned short y; // bits of x are in the even positions and y in the odd;
39 unsigned int z; // z gets the resulting 32-bit Morton Number.
40
41 z = MortonTable256[y >> 8] << 17 |
42 MortonTable256[x >> 8] << 16 |
43 MortonTable256[y & 0xFF] << 1 |
44 MortonTable256[x & 0xFF];
为了提高速度,可以使用一个额外的表,其值是MortonTable256预先向左移动一位。这第二个表可以用于y查找,从而减少了两个操作,但所需内存几乎增加了一倍。扩展这个同样的想法,可以使用四个表,其中两个表向左预移16位,这样我们总共只需要11次操作。
64位乘法的交错位
在11个操作中,这个版本将两个字节的位进行交错运算(而不是像其他版本那样进行短线运算),但很多操作都是64位的乘法,所以并不适合所有机器。输入参数x和y应该小于256。
unsigned char x; // Interleave bits of (8-bit) x and y, so that all of the
unsigned char y; // bits of x are in the even positions and y in the odd;
unsigned short z; // z gets the resulting 16-bit Morton Number.
z = ((x * 0x0101010101010101ULL & 0x8040201008040201ULL) *
0x0102040810204081ULL >> 49) & 0x5555 |
((y * 0x0101010101010101ULL & 0x8040201008040201ULL) *
0x0102040810204081ULL >> 48) & 0xAAAA;
Holger Bettag在2004年10月10日阅读了这里的基于乘法的位反转后,受到启发,提出了这个技术。
通过二进制魔数交错位
1 static const unsigned int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF};
2 static const unsigned int S[] = {1, 2, 4, 8};
3
4 unsigned int x; // Interleave lower 16 bits of x and y, so the bits of x
5 unsigned int y; // are in the even positions and bits from y in the odd;
6 unsigned int z; // z gets the resulting 32-bit Morton Number.
7 // x and y must initially be less than 65536.
8
9 x = (x | (x << S[3])) & B[3];
10 x = (x | (x << S[2])) & B[2];
11 x = (x | (x << S[1])) & B[1];
12 x = (x | (x << S[0])) & B[0];
13
14 y = (y | (y << S[3])) & B[3];
15 y = (y | (y << S[2])) & B[2];
16 y = (y | (y << S[1])) & B[1];
17 y = (y | (y << S[0])) & B[0];
18
19 z = x | (y << 1);
判断一个字是否有零字节
1 // Fewer operations: 2 unsigned int v; // 32-bit word to check if any 8-bit byte in it is 0 3 bool hasZeroByte = ~((((v & 0x7F7F7F7F) + 0x7F7F7F7F) | v) | 0x7F7F7F7F);
上面的代码在进行快速字符串复制时可能很有用,其中一次复制一个字;它使用了5次操作。另一方面,用明显的方法(下面的方法)测试空字节至少有7次操作(当用最节省的方法计算时),最多有12次。
1 // More operations:
2 bool hasNoZeroByte = ((v & 0xff) && (v & 0xff00) && (v & 0xff0000) && (v & 0xff000000))
3 // OR:
4 unsigned char * p = (unsigned char *) &v;
5 bool hasNoZeroByte = *p && *(p + 1) && *(p + 2) && *(p + 3);
本节开头的代码(标注为 "较少操作")的工作原理是首先将字中4个字节的高位清零。随后,它添加一个数字,如果任何一个低位被初始化,将导致一个字节的高位溢出。接下来,原始字的高位与这些值进行OR,因此,如果一个字节中的任何位被设置,则该字节的高位被设置。最后,我们通过在除高位以外的所有地方用1进行OR,并对结果进行反转来确定这些高位是否为零。扩展到64位是很简单的,只需将常量增加到0x7F7F7F7F7F7F7F7F7F7F。
为了进行额外的改进,可以执行一个只需要4次操作的快速预测试,以确定该字是否可能有一个零字节。如果高字节是0x80,测试也会返回真值,因此偶尔会出现假阳性,但上述较慢和更可靠的版本可以用于候选者,以便在正确输出的情况下全面提高速度。
bool hasZeroByte = ((v + 0x7efefeff) ^ ~v) & 0x81010100;
if (hasZeroByte) // or may just have 0x80 in the high byte
{
hasZeroByte = ~((((v & 0x7F7F7F7F) + 0x7F7F7F7F) | v) | 0x7F7F7F7F);
}
还有一种更快的方法--使用下面定义的hasless(v, 1);它在4个操作中工作,并且不需要子序验证。它简化为
#define haszero(v) (((v) - 0x01010101UL) & ~(v) & 0x80808080UL)
子表达式(v - 0x010101UL),每当v中对应的字节为0或大于0x80时,就会评估到任何字节中设置的高位。子表达式~v & 0x80808080UL,在v的字节没有设置高位的情况下(所以该字节小于0x80),就会评估到字节中设置的高位。最后,通过对这两个子表达式进行AND运算,结果是v中的字节为零的地方设置的高位,因为第一个子表达式中由于大于0x80的值而设置的高位被第二个子表达式掩盖掉了。
判断一个字的字节是否等于n
我们可能想知道一个词中的任何字节是否有特定的值。要做到这一点,我们可以用一个已经被我们感兴趣的字节值填满的字来XOR测试这个值。因为将一个值与自身进行 XOR 的结果是一个零字节,否则为非零,我们可以将结果传递给 haszero。
#define hasvalue(x,n) \
(haszero((x) ^ (~0UL/255 * (n))))
判断一个字的字节数是否小于n
测试一个字x是否包含一个值<n的无符号字节,具体来说对于n=1,可以通过一次检查一个长的字来找到0字节,或者通过先用掩码对x进行XOR,找到任何一个字节。当n为常数时,使用4次算术/逻辑运算。
要求:x>=0;0<=n<=128。
#define hasless(x,n) (((x)-~0UL/255*(n))&~(x)&~0UL/255*128)
要计算x中在7次操作中小于n的字节数,使用
#define countless(x,n) \
(((~0UL/255*(127+(n))-((x)&~0UL/255*127))&~(x)&~0UL/255*128)/128%255)
判断一个字是否有大于n的字节
测试一个字x是否包含一个值>n的无符号字节,当n为常数时,使用3个算术/逻辑运算。
要求:x>=0;0<=n<=127。
#define hasmore(x,n) (((x)+~0UL/255*(127-(n))|(x))&~0UL/255*128)
要计算x中在6次操作中超过n的字节数,使用:
#define countmore(x,n) \
(((((x)&~0UL/255*127)+~0UL/255*(127-(n))|(x))&~0UL/255*128)/128%255)
确定一个字是否有一个介于m和n之间的字节
当m<n时,该技术测试一个字x是否包含一个无符号字节值,这样m<值<n,当n和m为常数时,它使用7次算术/逻辑运算。
注意:等于n的字节会被 likelyhasbetween报告为假阳性,所以如果需要某种结果,应该按字符进行检查。
要求:x>=0;0<=m<=127;0<=n<=128。
#define likelyhasbetween(x,m,n) \
((((x)-~0UL/255*(n))&~(x)&((x)&~0UL/255*127)+~0UL/255*(127-(m)))&~0UL/255*128)
这种技术适合于快速预试。一种需要多做一次操作(对于常数m和n,共8次),但能提供准确答案的变式是:
#define hasbetween(x,m,n) \
((~0UL/255*(127+(n))-((x)&~0UL/255*127)&~(x)&((x)&~0UL/255*127)+~0UL/255*(127-(m)))&~0UL/255*128)
要在10次操作中计算x中介于m和n(独占)之间的字节数,使用:
#define countbetween(x,m,n) (hasbetween(x,m,n)/128%255)
计算词法上的下一个位排列组合
假设我们有一个整数中N位设为1的模式,我们想得到N个1位在词法意义上的下一个排列方式。例如,如果N为3,比特模式为00010011,那么接下来的模式将是00010101,00010110,00011001,00011010,00011100,00100011,以此类推。下面是一种快速计算下一个排列方式的方法。
unsigned int v; // current permutation of bits
unsigned int w; // next permutation of bits
unsigned int t = v | (v - 1); // t gets v's least significant 0 bits set to 1
// Next set to 1 the most significant bit to change,
// set to 0 the least significant ones, and add the necessary 1 bits.
w = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));
适用于 x86 CPU 的 __builtin_ctz(v) GNU C 编译器固有函数返回尾部零的数目。如果你使用的是微软的 x86 编译器,内在函数是 _BitScanForward。这两个编译器都会发出bsf指令,但其他架构可能会有等价的指令。如果没有,那么可以考虑使用前面提到的计算连续零位的方法之一
这里是另一个版本,由于它的除法运算符,速度往往比较慢,但它不需要计算尾部的零
unsigned int t = (v | (v - 1)) + 1;
w = t | ((((t & -t) / (v & -v)) >> 1) - 1);
浙公网安备 33010602011771号