zoukankan      html  css  js  c++  java
  • [翻译]位运算暗黑魔导书

    原文地址:http://graphics.stanford.edu/~seander/bithacks.html

    有本书叫hack‘s delight也是主讲位运算的。

    下面正文开始翻译吧,不定时更新。

    声明:

    如无特殊说明,所有的代码段都是不受版权限制的,如果喜欢,大家可以随便使用。文章内容由Sean Eron Anderson收集整理于1997-2005年。位运算的相关代码、说明以及分类希望能够给大家带来帮助,但是不对于该内容是否适合特定用途提供任何担保(or授权?)。截止2005年5月5日,所有的代码都已经彻底的测试过了。成千上万人阅读过这些内容。此外,卡内基梅隆大学计算机科学学院院长Randal Bryant教授亲自使用他的Uclid代码检证系统测试了几乎所有的代码。对于他没有测试的部分,我使用所有可能的输入在32位机上进行了检证。第一个发现代码中确实存在BUG的人,作者原意通过支票或者Paypal支付10美刀作为赏金,或者是向慈善机构捐赠20美刀。

    关于运算次数的计算方法

    当计算某一算法的总操作次数时,任何一个C语言的运算符被当做一条操作。写内存的中间操作不被计算在内。当然,这种操作数计算方法仅仅是为了得到机器指令和CPU占用时间的近似值。我们只是假设每一条不同操作耗费的时间是相同的(事实上并不是这样的),但是CPU的确已经在逐渐的朝着这个方向发展了。有很多细节会影响到在某个系统上一段代码的运行时间,比如如缓存大小,内存带宽,指令集等。最终,基准测试(benchmarking)是判断一个算法是比其他算法高效的最好办法。所以请在你的目标机器上测试使用下面技巧的可能性。(后半段一大堆,总结一下,实践是检验真理的唯一标准)。

    判断一个整数的符号

    int v; // 判断v的符号
    int sign; // 判断的结果存在这里 

    // CHAR_BIT 一个字节(byte)包含的位数(bit) (通常是8).
    sign = -(v < 0); // if v < 0 then -1, else 0. 
    // 方法2 避免通过使用CPU寄存器的标志位进行分支(IA32):
    sign = -(int)((unsigned int)((int)v) >> (sizeof(int) * CHAR_BIT - 1));
    // 方法3 再减少一条指令 (但不具有可移植性):
    sign = v >> (sizeof(int) * CHAR_BIT - 1); 

    上面最后一个表达式等效于:sign = v >> 31(用于计算32位整数的符号),这种方法只进行一次运算,比通常的sign = -(v < 0)要快的多。这种方法可以奏效是因为当有符号的整数v进行右移时,最左边的位被拷贝到了其他位.当v的值为负数时,最左边的位为1,否则的话最左边的位为0(v为整数或者0时);所有的位为1表示的数为-1(-1在内存中为0xffffffff)。不幸的是,这种方法适用于特定的架构(architecture-specific)。

    除此之外,如果你想要的结果是-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

    注意:2003年3月7日,Angus Duggan指出,1989年的ANSI C规范允许有符号数右移的结果由实际情况决定(implementation-defined,详见下文扩展说明),即取决于编译器,所以在一些系统这个方法(hack)可能不奏效。Toby Speight于2005年9月28日建议为了提高代码的可移植性,此处使用CHAR_BIT而不是认为一个字节的长度就是8位。Angus于2006年3月4日推荐了上面更具有移植性的计算方法。Rohit Garg于2009年9月12日给出了判断非负数版本的算法。

    【扩展说明】(源自互联网资料,进行了简单的总结归纳,未做深入探究):

    C标准中没有做明确规定的地方会用Implementation-defined、Unspecified或Undefined来表述,中文资料中很多地方把这三种情况统称为"未明确定义"的。但实际上是由区别的:

    [undefined]
    未定义行为是在某些不正确的情况下,标准并未规定应该怎样做,实现可以采取任何行动,也可以什么都不做。
    如当有符号整数溢出时该采取什么行动。
    Undefined的情况是完全不确定的,C标准没规定怎么处理,编译器很可能也没规定,甚至也没做出错处理,有很多Undefined的情况编译器是检查不出来的,最终会导致运行时错误,比如数组访问越界就是Undefined的。

    [unspecified]
    未确定行为是在某些正确的情况下,标准并未规定应该怎样做。
    对于Unspecified(未详细说明,未指定)的情况,往往有几种可选的处理方式,C标准没有明确规定按哪种方式处理,编译器可以自己决定,并且也不必写在编译器的文档中,这样即便用同一个编译器的不同版本来编译也可能得到不同的结果,因为编译器没有在文档中明确写它会怎么处理,那么不同版本的编译器就可以选择不同的处理方式,比如一个函数调用的各个实参表达式按什么顺序求值是Unspecified的(函数参数的求值顺序)。

    [implementation-defined]
    实现定义的行为是由编译器设计者决定采取何种行动,并写入实用手册。
    如当一个整型数向右移位时,要不要扩展符号位。
    如C标准没有明确规定char是有符号的还是无符号的,但是要求编译器必须对此做出明确规定,并写在编译器的文档中。

    比较两个整数的符号是否相反

    int x, y; // 待比较的两变量
    bool f = ((x ^ y) < 0); // 如果x和y的符号相反,结果为真

    Manfred Weis 于2009年11月26日建议作者添加了这个条目。

    不使用分支,计算一个整数的绝对值

    int v; // 求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 要快得多,虽然他们的操作数是相同的。
    2003年3月7日,Angus Duggan指出在1989年的ANSI C规范中,有符号的整数右移运算的结果是未明确定义的(implementation-defined,取决于编译器),所以在有些系统上这个方法是不奏效的。作者曾经在ANSI C规范中发现,ANSI C并没有要求数值以二进制补码(two's complement)的形式表示,所以出于这个原因上面的方法可能也不适用(再非常少数的老机器上仍然使用二进制反码(one's complement)的形式)。2004年3月14日,Keith H. Duggar发给了作者上面已被申请专利的衍生版的表达式(patented variation);这是我最早提出的r=(+1|(v>>(sizeof(int)*CHAR_BIT-1)))*v,的升级版,因为不需要使用乘法运算了。但是不幸的是Vladimir Yu Volkonsky已经于2000年6月6日在美国申请了注册了专利并分配给了Sun Microsystems。2006年8月13日,Yuriy Kaminskiy 告诉我专利可能是无效的因为在申请专利之前这个方法就已经被公开发表,比如1996年11月9日Agner Fog的《How to Optimize for the Pentium Processor》(如何为奔腾处理器做优化)。Yuriy还提到这个文档1997年的时候被翻译成了俄文,Vladimir可能读过了。此外,网站时光倒流机器(The Internet Archive)也存在到这里的旧链接。2007年1月30日,Peter Kankowski 给作者分享了一个版本的求绝对值,灵感来源于Microsoft's Visual C++编译器的输出。在这里作为最终解决方案。2007年12月6日Hai Jin抱怨运算的结果是有符号的,所以当计算most negative value时,得到的结果还是负数。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 的值首先将负数v加上2**N(2的N次方)转换为无符号的数,得到了一个v的二进制补码我们称之为U,然后U与我们想要的结果正好相反。即 -U = 0 - U = 2**N - U = 2**N - (v+2**N) = -v = abs(v). 

    不用分支结构,求两个整数中的最小值或最大值

    int x; // 求x y中的最小值
    int y; 
    int r; // 存放结果

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

    在很少一部分机器中,程序分支是非常耗费资源的,而且不支持条件跳转指令(condition move instructions),上面的方法可能会比显而易见的方法r = (x < y) ? x : y更高效,虽然上面的方法多两条指令。(通常情况下,还是显而易见的方法最好用)这个方法奏效是因为
    如果 x < y ,那么 - ( x < y ) 的每一位(bit)都是1,那么 r = y ^ (x ^ y) & ~0 = y ^ x ^ y = x 。在一些机器上,判断 (x < y) 是 0 还是 1 需要分支指令,所以这样可能就没有优势了。

    求最大值的方法为: 

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

    更快更贱(dirty,或者说是不稳定的?)的版本:

    如果可以确定 INT_MIN <= x - y <= INT_MAX, 那么你使用下面的方法会更快因为 (x - y) 只需要一次运算。

    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规范中对有符号数的右移结果没有明确规定吗,所以这个方法不具有可移植性。如果溢出时抛出异常,那么 x 和 y 的值需要是无符号的或者被转换成无符号的以便在进行减法时避免不必要的异常抛出,然后,右移需要需要有符号的数作为操作对象,这样才能使负数的所有位(bit)为1,所有再转换为有符号的。

    2003年3月7日,Angus Duggan 指出右移运算的可移植性问题。2005年5月3日,Randal E. Bryant 提醒了我需要 INT_MIN <= x - y <= INT_MAX 作为前提,同时建议了非快速、不稳定的版本方案(non-quick and dirty version)。 这些问题都需要在更快更贱(quick and dirty version)的方法中考虑。2005年7月6日 Nigel Horspoon 发现在奔腾处理器上使用GCC编译器编译显而易见的算法,因为处理器和编译器对(x < y)语句的处理方法而生成了同样的程序。2008年7月9日,Vincent Lefèvre 指出之前的 r = y + ((x - y) & -(x < y)) 版本在进行减法运算时有因为溢出抛出异常的潜在可能。2008年7月9日 Timothy B. Terriberry 建议使用异或运算(xor)代替加法和减法运算以避免类型转换和溢出的风险。

    确定一个整数是不是2的幂

    unsigned int v; // 判断 v 是不是 2 的幂
    bool f; // 结果放在这里 

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

    注意 0 在这里被错误的认为是 2 的幂了。想要避免这个问题,方法如下: 
    f = v && !(v & (v - 1));

    固定位宽(有符号数)的符号扩展

    对于内建数据类型的符号扩展是自动完成的,比如char型和int型,但是假设你有一个以二进制补码形式存储的数 x ,在存储时只存储 b bit。然后,假设你想要把 x 转换成一个超过 b bit的整数,对于整数使用简单的拷贝即可,但是如果是负数,符号位必须要扩展。举例说明,如果我们只使用4个bit存储一个数字,那么 -3 表示为二进制数 1101 。如果我们使用8 bit,那么 -3 为 1111 1101 。将 4-bit 表示法数的最高位,重复的填充在原数左边的各个位置上,以此实现使用更多的bit来表示一个数,这就是符号扩展。在C语言中,固定位宽数的扩展符号是非常简单的,因为我们可以在结构体(struct)或者联合体(union)中使用位字段(bit fields)。比如将一个 5 bit宽的数扩充为完全的整型:

    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++ 模板(template)功能通过一条操作将B bit位宽数据实现符号扩展(当然,编译器生成了更多)。

    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

    2005年5月2日John Byrd 发现了代码中的一个排版错误(由于HTML格式的原因)。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;

    上面的代码需要四步操作,但是当位宽是固定的而不是可变的时候,假设高位都已经是 0 了,那么只需要两步快速的操作。一个更简洁的但是不具有可移植性,并且不需要 x中 b 位的位高位全部为 0 的方法如下:

    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 除去 b bit 以外的部分为1的情况。2009年12月31日 Chris Pirazzi 建议作者增加了较快算法的版本,这样对于固定位宽数的符号扩展只需要两步操作,可变位宽需要三步操作。

    3步完成可变位宽(有符号数)的符号扩展

    由于这种实现方法需要使用乘法和除法运算,所以在一些机器上可能会很慢。这个版本需要4步操作。如果你已经知道初始的位宽 b是大于1的,你可以通过使用语句 r = (x * multipliers[b]) / multipliers[b] 只需要查一次表并在3步内完成符号扩展。

    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
    #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)
    }; // (add more if using more than 64 bits)
    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)
    }; // (add more for 64 bits)
    #undef M
    r = (x * multipliers[b]) / divisors[b];

    下面衍生出来的方法不具有可移植性,但是在一些构架中,进行算术右移运算,符号位会被保持,这种方法就会很快。

    const int s = -b; // OR: sizeof(x) * CHAR_BIT - b;
    r = (x << s) >> s;

    2005年5月3日 Randal E. Bryant 在早期的版本上(使用 multipliers[]代替 divisors[] )指出了当 x=1 并且 b=1 时存在BUG。

    PS:根据作者的意思,multipliers[b] 在进行查表时,为一次操作(operation)。乘法为一次操作,除法为一次操作。所以现在的完整版实际是4部完成可变位宽有符号数的符号扩展。只有初始位宽大于1的时候, multipliers[] 和 divisors[] 成了一样的表,减少了一次查表的操作,成为了3次操作完成的版本。

    不使用分支,判断进行设置位或清除位的运算

    bool f; // conditional flag
    unsigned int m; // the bit mask
    unsigned int w; // the word to modify: if (f) w |= m; else w &= ~m; 

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

    // 或者在超标量(superscalar)CPU上:
    w = (w & ~m) | (-f & m);

    在某些架构上,避免使用分支结构似乎可以提高两倍左右的效率。例如,在AMD 速龙™ XP 2100+ CPU上的通常速度测试可以提速5-10%。使用Intel 酷睿 2 双核处理器运行超标量版的程序会比第一个程序快16%。2003年12月11日 Glenn Slayden 向作者提供了第一个表达式。2007年4月3日 Marco Yu 分享了超标量版本的程序并且在2天后提醒了作者一个输入错误。

    PS:
    其实这个算法就是想完成如下语句:
    if (f) 
            w |= m; 
    else 
            w &= ~m; 
    如果f是真,目标数w通过掩码m完成相应位设置为1的操作;如果f为假,目标数w通过掩码m完成相应位的清零操作。
    使用作者提供的这种方法,可以避免分支结构的程序,在一些机器上分支结构是非常耗费资源的。。。

    超标量体系结构 (superscalar architectures)
    超标量体系结构描述一种微处理器设计,它能够在一个时钟周期执行多个指令。在超标量体系结构设计中,处理器或指令编译器能够判断指令能独立于其它顺序指令而执行,还是依赖于另一指令,必须跟其按顺序执行。处理器然后使用多个执行单元同时执行两个或更多独立指令。超标量体系结构设计有时称“第二代RISC”。

    不使用分支,根据判断条件求一个数的相反数

    如果你想要当Flag为false时求一个数的相反数,那么下面的方法可以避免分支结构:

    bool fDontNegate; // Flag indicating we should not negate v.
    int v; // Input value to negate if fDontNegate is false.
    int r; // result = fDontNegate ? v : -v;

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

    如果你需要在flag为True的时候求相反数,那么:

    bool fNegate; // Flag indicating if we should negate v.
    int v; // Input value to negate if fNegate is true.
    int r; // result = fNegate ? -v : v;

    r = (v ^ -fNegate) + fNegate;

    2009年6月2日 Avraham Plotnitzky 建议我添加了第一版的程序。为了避免乘法运算,2009年6月8日作者想出了第二版的程序。2009年11月26日Alfonso De Gregorio 指出代码中缺少括号,得到了一笔bug赏金。

    根据一个掩码完成两个数值的位合并

    unsigned int a; // value to merge in non-masked bits
    unsigned int b; // value to merge in masked bits
    unsigned int mask; // 1 where bits from b should be selected; 0 where from a.
    unsigned int r; // result of (a & ~mask) | (b & mask) goes here

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

    这种方法使根据一个掩码将两个变量位合并的通常算法减少了一步操作。((a & ~mask) | (b & mask)共四步操作,a ^ ((a ^ b) & mask)共三步操作。)如果掩码是一个常量,那么这种方法就没有优势了。本条内容为2006年2月9日Ron Jeffery向作者提供。

    有效位统计 (土方法)

    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次迭代。

    PS:有效位即二进制表示时该位(bit)为1

    查表法进行有效位统计

    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; // count the number of bits set in 32-bit value v
    unsigned int c; // c is the total bits set in v

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

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


    // To initially generate the table algorithmically:
    BitsSetTable256[0] = 0;
    for (int i = 0; i < 256; i++)
    {
      BitsSetTable256[i] = (i & 1) + BitsSetTable256[i / 2];
    }

    2009年7月14日 Hallvard Furuseth 建议使用宏对表进行压缩。


    PS:

    查表法不是很高端,标准的消耗便宜的ROM来节省RAM和CPU的手段,但是 Hallvard Furuseth 对表进行压缩的宏看起来很高端的样子,不知道是出于习惯还是什么原因作者把宏放进了table表里面定义,看起来怪怪的,使用GCC对宏的作用域实测,与括号无关,只要#define在有效代码段,在该文件#define行之后任意位置均可使用该宏。

    另外,表里面真实的样子的。。。:
    static const unsigned char BitsSetTable256[256] = 
    {
    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8 };

    Brian Kernighan 法进行有效位统计

    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位字长的变量最高位有效,那么也只需要进行一次循环即可得出结果。

    1988年出版的C《The C Programming Language 2nd Ed》(by Brian W. Kernighan and Dennis M. Ritchie《C 程序设计语言第二版》)在练习2-9(exercise 2-9)中提到了这种方法。2006年4月19日Don Knuth指出这种方法“最早由Peter Wegner发表于CACM 3 (1960),322(Communications of the ACM (CACM) 创刊于1957年的月刊)。(同样发现Derrick Lehmer 1964年独立发表在由Beckenbach编辑的书中)”。

    对14,24,,32位变量使用64位指令进行有效位统计

    unsigned int v; // count the number of bits set in v
    unsigned int c; // c accumulates the total bits set in v

    // 方法1, 最大 14位的变量 v:
    c = (v * 0x200040008001ULL & 0x111111111111111ULL) % 0xf;

    // 方法2, 最大 24位的变量 v:
    c = ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
    c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) 
         % 0x1f;

    // 方法3, 最大 32位的变量 v:
    c = ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
    c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) % 
         0x1f;
    c += ((v >> 24) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;

    这种方法需要使用能够进行快速取模和除法的64位CPU才能具有较高的效率。方法1只需要3步操作;方法2需要10部操作;方法3需要15步操作。
    Rich Schroeppel最早提出来9位数据的版本,和方法1类似;详见Beeler, M., Gosper, R. W., and Schroeppel, R. HAKMEM.1972年2月29日撰写的“ MIT AI Memo No.239”的Programming Hacks部分。由这些带来的灵感,Anderson. Randal E完成了上面的各种衍生方法。2005年5月3日Bryant提出了一对bug的修复方法。2007年2月1日Bruce Dawson对12位版本稍作修改,使14位变量可以使用相同的数据进行计算(适用于14位数)。

    并行计算有效位

    unsigned int v; // count bits set in this (32-bit value)
    unsigned int c; // store the total here
    static const int S[] = {1, 2, 4, 8, 16}; // Magic Binary Numbers 传说中的魔数,一个叫S,一个叫B。作者大学可能选修mandarin的。。。
    static const int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF, 0x0000FFFF};

    c = v - ((v >> 1) & B[0]); //我觉得是不是应该写成c = v - ((v >> S[0]) & 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以二进制表示如下: 
    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)),(ceil为进一法取整), 然后我们进行计算需要的表达式次数和S或者B一样大。 一个32位的数 v, 需要进行16次运算. 
    T计算32位整数 v 的有效位的最佳方法如下: 

    v = v - ((v >> 1) & 0x55555555); // 把输入值作为临时变量再次使用
    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

    请移步至Ian Ashdown's nice newsgroup post查看关于有效位计算的更多内容(也称为横向加法 sideways addition)。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 又减少了3步。2006年1月8日 Eric Cole 指出作者在 Don 提出的方法中出现了一个输入错误。之后 Eric 在2006年11月17日将有效位统计的最佳方法方法进行了任意位宽泛化。2007年4月5日 Al Williams 发现了作者在上面第一个方法中存在一行无效代码。

    计算从最高位到指定位之间的有效位

    下面的程序将返回二进制数从最高位到指定位之间1的个数。

     
     uint64_t v; // Compute the rank (bits set) in v from the MSB(最高位) to pos.
      unsigned int pos; // Bit position to count bits upto.
      uint64_t r; // Resulting rank of bit at pos goes here.

      // Shift out bits after given position.
      r = v >> (sizeof(v) * CHAR_BIT - pos);
      // Count set bits in parallel.
      // 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 向作者提供了这个方法,是下一节《从最高位开始,找到第 r 个有效位的位置(computing the bit position with the given rank)》的逆运算。

    从最高位开始,找到第 r 个有效位的位置。

    下面的基于64位的代码用于找到从左边开始第 r 个被设置为 1 的位所在的位置。换句话说,如果我们从最高位开始,向右进行检索,累计被设置为 1 的位的个数,直到我们想要的第 r 个,那么我们停下来的位置就是我们想要的结果(返回值)。如果我们想要找的 r 超过了有效位的总数,那么返回 64 。代码同样可以修改为 32 位数或者从右边开始计数。

      uint64_t v; // Input value to find position with rank r.
      unsigned int r; // Input: bit's desired rank [1-64].
      unsigned int s; // Output: Resulting position of bit with rank r [1-64]
      uint64_t a, b, c, d; // Intermediate temporaries for bit count.
      unsigned int t; // Bit count temporary.

      // Do a normal parallel bit count for a 64-bit integer, 
      // but store all intermediate steps. 
      // 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);
      // Now do branchless select! 
      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)上分支程序运行速度很快的话,可以考虑打开注释掉的 if 语句,然后把之前的语句注释掉。
    2009年11月21日 Juha Järvi 向作者提供了这个算法。

    判断有效位个数奇偶性的原始方法

    unsigned int v; // word value to compute the parity of
    bool parity = false; // parity will be the parity of v (奇数返回1,偶数返回0)

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

    上面的代码实现,使用了类似 Brian Kernigan 有效位统计法的算法。消耗的时间与有效位的个数成正比。

    查表法计算有效位个数的奇偶性

    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; // byte value to compute the parity of
    bool parity = ParityTable256[b];

    // OR, for 32-bit words:
    unsigned int v;
    v ^= v >> 16;
    v ^= v >> 8;
    bool parity = ParityTable256[v & 0xff];

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

    2005年5月3日 Randal E. Bryant 提出了使用指针 p的衍生版。2005年9月27日 Bruce Rawles 发现了表中变量名称的一个拼写错误,他因此得到了10美刀的bug赏金。2006年10月9日 Fabrice Bellard 建议了上面的 32位字宽变量, 只需要查一个表; 之前版本的程序需要4次查表(每个字节一次) 并且更慢. 2009年7月14日 Hallvard Furuseth 建议使用宏来压缩表。

    使用64位乘法和除法取模计算有效位个数的奇偶性

    unsigned char b; // byte value to compute the parity of
    bool parity = 
      (((b * 0x0101010101010101ULL) & 0x8040201008040201ULL) % 0x1FF) & 1;

    上面的方法大约需要 4 步操作, 但是只能计算一个字节位宽的数据.

    使用乘法运算计算一个字(word 此处为32位宽)中有效位的奇偶性

    下面的方法使用乘法运算,仅仅需要8步操作可以计算一个32位宽的变量中有效位的奇偶性。

        unsigned int v; // 32-bit word
        v ^= v >> 1;
        v ^= v >> 2;
        v = (v & 0x11111111U) * 0x11111111U;
        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;

    Andrew Shapira 2007年9月2日 实现了这种方法并提供给作者。

    并行计算有效位个数的奇偶性

    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 次操作完成对字节(byte)的计算。这种方法通过通过位移把 32 位变量的 8个 nibble(半字节,4比特)异或在一起。 然后, 二进制数 0110 1001 1001 0110 (16进制的 0x6996 ) 被右移的位数为变量 v 的低四位表示的数值 . 右移异或后的v的低四位,就像是索引一样,作用于小型的 16位奇偶性表。结果的所具有的有效位奇偶性与 v 是一样的,然后通过与掩码运算后返回。

    感谢 Mathew Hendry 2002年12月15日最后指出了移位查表的想法。查表的优化使得使用位移后异或的方法又消减了两步操作。

    使用减法和加法运算交换两个变量的值

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

    这种交换 a 和 b 值的方法不需要使用临时变量。初始时 a 和 b 的变量在内存中的地址是否相同的检查是可以省略的,前提是你确定这种情况不会发生。(编译器在进行优化时可能会不管什么情况都省略掉这部分)。如果你启用了溢出异常,那么请使用无符号的变量以保证不会抛出异常。下面使用异或运算的方法在某些机器上可能稍微快一点。 不要对浮点数使用 (除非你对他们的原始整数表示方式进行操作). 
    Sanjeev Sivasankaran 2007年6月12日建议作者添加该条. Vincent Lefèvre 2008年7月9日指出了可能溢出的问题。

    使用异或运算交换两个变量的值 

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

    这是一个不使用额外的存储空间作为临时变量来交换a 和b 的值的老把戏了。
     
    2005年1月20日 Iain A. Fleming 指出如果交换的变量存在相同的内存中,这个宏就不起作用了,比如 SWAP(a[i], a[j]) with i == j 。所以如果这种情况可能发生,请考虑将宏定义为:(((a) == (b)) || (((a) ^= (b)), ((b) ^= (a)), ((a) ^= (b)))). 2009年7月14日, Hallvard Furuseth 建议在一些机器上 (((a) ^ (b)) && ((b) ^= (a) ^= (b), (a) ^= (b))) 可能会更高效,由于 (a) ^ (b) 表达式会被复用。

    使用异或运算交换数据中的某一段连续位序列

    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 (二进制表示)连续的长度为 n = 3 的序列进行交换,从 i = 1 的起始位置(右数第二位)和 j = 5 的起始位置做长度为3的连续二进制序列位交换。结果为 r = 11100011 (二进制).这个交换方法和单纯的异或运算进行交换的把戏很相似,但是被用来实现交换特定的位段。变量 x 存储了待进行交换的两段位序列的异或结果,然后在通过使用这个异或的结果 x 将各段的数据交换完成。当然,如果待交换的序列有重叠的部分,那么结果是不定值。 

    2009年7月14日 Hallvard Furuseth 建议作者将 1 << n to 1U << n ,因为这样变量会被指定放在无符号的区域以避免位移运算时受到符号位的影响。

    二进制逆序的通常方法

    unsigned int v; // 待进行二进制逆序的数
    unsigned int r = v; // r 将会用来保存 v 的二进制逆序结果; 首先得到 v 的最低有效位(LSB = Least Significant Bit)
    int s = sizeof(v) * CHAR_BIT - 1; // 最后要额外进行的位移数

    for (v >>= 1; v; v >>= 1)

      r <<= 1;
      r |= v & 1;
      s--;
    }
    r <<= s; // 移动的位数s 即 v 的最高位为0的位的个数

    2004年10月15日 Michael Hoisie 指出了原始版本上得一个bug。Randal E. Bryant 2005年5月3日建议删除了一个额外的运算。Behdad Esfabod 2005年5月18日提出稍微进行修改可以使循环中减少一次迭代。接着在2007年2月6日 Liyong Zhou 建议提出了较大改进的版本,当 v 不是 0时进行循环,这样可以不需要对所有位进行迭代可以较早的结束循环。

    查表法对一个字(32bit)的二进制序列进行逆序 

    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; // 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 位乘法和除法取模运算): 

    unsigned char b; // 逆序这个 (8位) 字节
     
    b = (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;

    乘法运算产生了五个独立的8位数据的副本,用于扇出到一个64位的数据中。与运算删选出了需要进行翻转的位的位置,以每10个位为单位进行。乘法和与运算将原有的一个字节数据中的每一位会在长度为10的位组中存在一位。原字节中倒序之后的位置与某一个10位宽的域中的值是一致的(一一对应的)。最后一步的取模 2^10 - 1 运算的功能是将64位变量中各个10位宽的数据合并在一起(从0-9, 10-19, 20-29, ...)。他们没有重叠的部分,所以对多出的一步取模运算的效果和与运算类似。这种方法出自Rich Schroeppel的《Beeler, M., Gosper, R. W., and Schroeppel, R. HAKMEM. MIT AI Memo 239, Feb. 29, 1972》中的Programming Hacks 部分。

    PS:
    0x0884422010 = 0000000001 0000100010 0001000100 0010001000 0000010000

     4 步运算实现一个字节的二进制位序列逆序 (64 位乘法,不需要除法运算): 

    unsigned char b; // 逆序这个字节
     
    b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;

    下面看起来貌似火星文的部分展示了一个由布尔量 a, b, c, d, e, f, g, 和 h 组成的8位二进制数据的变化过程。注意第一个乘法运算将二进制数据扇出了多个拷贝的副本,然后最后的乘法将从右开始的连续五个字节合并到了一起。
     

     


                                                                                            

    注意最后两步在一些处理器上可以合并,因为在这些处理器上寄存器可以进行字节访问;仅进行乘法运算使寄存器存储了结果的高32位然后取出低8位。这样,可能只需要6部运算了。
    Devised by Sean Anderson, July 13, 2001. 

    7步运算实现一个字节的二进制位序列逆序(非64位):

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

    确保你的结果被存储或者转换为一个无符号的char型变量以确保将高位中的无用数据删除。
    Devised by Sean Anderson, July 13, 2001. 
    Mike Keith 2002年1月3日指出并修正了一个排版错误。 

    5 * lg(N) 步计算并行完成 N位二进制数据的逆序

    unsigned int v; // 32-bit word to reverse bit order

    // 交换奇数位和偶数位
    v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
    // 交换连续的一对数据
    v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
    // 交换半个字节
    v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
    // 交换字节
    v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
    // 以2字节为单位进行交换
    v = ( v >> 16 ) | ( v << 16);

    接下来的变化依然是 O(lg(N)), 但是这种方法需要较多的运算次数来逆序 v. 它的优点在于计算常量时运行速度较快并且可以使用更少的存储空间(memory)。

    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位整数(或者更大的数),那么你需要增加很多行代码(按照例子中的模式);否则只会低32位被逆序了并且结果存在低32位中。
    请参阅 Dobb's Journal 1983年撰写的《Edwin Freed's article on Binary Magic Numbers》得到更多的信息。衍生出来的第二个方法是由 Ken Raeburn 2005年9月13日提出的。Veldmeijer 2006年3月19日提出在第一个版本的程序上最后一行表达式的与运算可以省略。

    不使用除法运算计算一个变量取模 1 << s 的结果

    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); 

    大部分程序员很早就学会了这个技巧,此处包含这个是为了收录算法的全面性。

    不使用除法运算计算一个变量取模(1 << s)-1的结果

    unsigned int n;                      // numerator
    const unsigned int s;                // s > 0
    const unsigned int d = (1 << s) - 1; // d可能是后面其中一个数值: 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中其中一个值, but since with modulus division
    // 当m等于d的时候,m为0.
    m = m == d ? 0 : m;

    This method of modulus division by an integer that is one less than a power of 2 takes at most 5 + (4 + 5 * ceil(N / s)) * ceil(lg(N / s)) operations, where N is the number of bits in the numerator. In other words, it takes at most O(N * lg(N)) time.

    Devised by Sean Anderson, August 15, 2001. Before Sean A. Irvine corrected me on June 17, 2004, I mistakenly commented that we could alternatively assign m = ((m + 1) & d) - 1; at the end. Michael Miller spotted a typo in the code April 25, 2005.

  • 相关阅读:
    基于日志实现ssh服务防护脚本
    给多主机传输文件脚本
    Python map() 函数
    常用希腊字母读法
    Python中Numpy.nonzero()函数
    numpy.eye() 生成对角矩阵
    (sklearn)岭回归 sklearn.linear_model.Ridge用法
    pip配置永久国内源
    python numpy的var std cov研究
    numpy中的mean()函数
  • 原文地址:https://www.cnblogs.com/foohack/p/3772910.html
Copyright © 2011-2022 走看看