Bit Twiddling Hacks

本文为转载,计划翻译中


《二进制位运算教程》

Bit Twiddling Hacks

By Sean Eron Anderson
seander@cs.stanford.edu

Individually, the code snippets here are in the
public domain
(unless otherwise noted) — feel free to use them however you please.  
The aggregate collection and descriptions are
© 1997-2005 Sean Eron Anderson.  
The code and descriptions are distributed
in the hope that they will be useful,
but WITHOUT ANY WARRANTY and
without even the implied warranty of
merchantability or fitness for a particular purpose.  
As of May 5, 2005, all the code has been tested thoroughly.  Thousands of
people have read it.  Moreover, Professor Randal Bryant, the Dean of Computer Science
at Carnegie Mellon University, has personally tested almost everything
with his Uclid code
verification system
.  What he hasn't tested, I have checked against all
possible inputs on a 32-bit machine.To the first person to inform me of a legitimate bug
in the code, I'll pay a bounty of US$10 (by check or Paypal)
.
If directed to a charity, I'll pay US$20.  

Contents

  • Computing parity (1 if an odd number of bits set, 0 otherwise)

  • Swapping Values

  • Reversing bit sequences

  • Modulus division (aka computing remainders)

  • Finding integer log base 2 of an integer (aka the position of the highest bit set)

  • Find integer log base 10 of an integer

  • Find integer log base 10 of an integer the obvious way

  • Find integer log base 2 of a 32-bit IEEE float

  • Find integer log base 2 of the pow(2, r)-root of a 32-bit IEEE float (for unsigned integer r)

  • Counting consecutive trailing zero bits (or finding bit indices)

  • Round up to the next highest power of 2 by float casting

  • Round up to the next highest power of 2

  • Interleaving bits (aka computing Morton Numbers)

  • Testing for ranges of bytes in a word (and counting occurances found)

  • Compute the lexicographically next bit permutation


  • About the operation counting methodology

    When totaling the number of operations for algorithms here,
    any C operator is counted as one operation.  
    Intermediate assignments, which need not be written to RAM, are
    not counted.
    Of course, this operation counting approach only serves as an
    approximation of the actual number of machine instructions and CPU time.
    All operations are assumed to take the same amount of time, which
    is not true in reality, but CPUs have been heading increasingly in this
    direction over time.  There are many nuances that determine
    how fast a system will run a given sample of code, such as cache sizes,
    memory bandwidths, instruction sets, etc.  In the end, benchmarking is
    the best way to determine whether one method is really faster than another,
    so consider the techniques below as possibilities to test on your target
    architecture.


    Compute the sign of an integer

    int v;      // we want to find the sign of v
    int sign;   // the result goes here 
    
    // CHAR_BIT is the number of bits per byte (normally 8).
    sign = -(v < 0);  // if v < 0 then -1, else 0. 
    // or, to avoid branching on CPUs with flag registers (IA32):
    sign = -(int)((unsigned int)((int)v) >> (sizeof(int) * CHAR_BIT - 1));
    // or, for one less instruction (but not portable):
    sign = v >> (sizeof(int) * CHAR_BIT - 1);

    The last expression above evaluates to sign = v >> 31
    for 32-bit integers.  
    This is one operation faster than the obvious way, sign = -(v < 0).
    This trick works because when signed integers are shifted right, the value
    of the far left bit is copied to the other bits.  The far left bit is 1
    when the value is negative and 0 otherwise; all 1 bits gives -1.  
    Unfortunately, this behavior is architecture-specific.

    Alternatively, if you prefer the result be either -1 or +1, then use:

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

    On the other hand, if you prefer the result be either -1, 0, or +1, then use:

    sign = (v != 0) | -(int)((unsigned int)((int)v) >> (sizeof(int) * CHAR_BIT - 1));
    // Or, for more speed but less portability:
    sign = (v != 0) | (v >> (sizeof(int) * CHAR_BIT - 1));  // -1, 0, or +1
    // Or, for portability, brevity, and (perhaps) speed:
    sign = (v > 0) - (v < 0); // -1, 0, or +1

    If instead you want to know if something is non-negative, resulting in +1 or
    else 0, then use:

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

    Caveat:  On March 7, 2003, Angus Duggan pointed out that the 1989 ANSI C
    specification leaves the result of signed right-shift implementation-defined,
    so on some systems this hack might not work.  For greater portability,
    Toby Speight suggested on September 28, 2005 that CHAR_BIT be used here
    and throughout rather than assuming bytes were 8 bits long.  Angus recommended
    the more portable versions above, involving casting on March 4, 2006.Rohit Garg suggested the version
    for non-negative integers on September 12, 2009.


    Detect if two integers have opposite signs

    int x, y;               // input values to compare signs
    
    bool f = ((x ^ y) < 0); // true iff x and y have opposite signs

    Manfred Weis suggested I add this entry on November 26, 2009.


    Compute the integer absolute value (abs) without branching

    int v;           // we want to find the absolute value of v
    unsigned int r;  // the result goes here 
    int const mask = v >> sizeof(int) * CHAR_BIT - 1;
    
    r = (v + mask) ^ mask;

    Patented variation:

    r = (v ^ mask) - mask;

    Some CPUs don't have an integer absolute value instruction (or the
    compiler fails to use them).  On machines where branching is expensive,
    the above expression can be faster than the obvious approach,
    r = (v < 0) ? -(unsigned)v : v, even though the number of operations
    is the same.

    On March 7, 2003, Angus Duggan pointed out that the 1989 ANSI C
    specification leaves the result of signed right-shift implementation-defined,
    so on some systems this hack might not work.  I've read that ANSI C does not
    require values to be represented as two's complement, so it may not work
    for that reason as well (on a diminishingly small number of old machines
    that still use one's complement).

    On March 14, 2004, Keith H. Duggar sent me the patented variation above; it is
    superior to the one I initially came up with,r=(+1|(v>>(sizeof(int)*CHAR_BIT-1)))*v,
    because a multiply is not used.  
    Unfortunately, this method has been patented in the USA on June 6, 2000 by Vladimir Yu Volkonsky and
    assigned to Sun Microsystems.

    On August 13, 2006, Yuriy
    Kaminskiy told me that the patent is likely invalid because the method
    was published well before the patent was even filed, such as inHow to Optimize for
    the Pentium Processor
    by Agner Fog, dated November, 9, 1996.  Yuriy also
    mentioned that this document was translated to Russian in 1997, which
    Vladimir could have read.  Moreover, the Internet Archive also has an oldlink to it.

    On January 30, 2007, Peter Kankowski shared with me anabs version he discovered that was inspired by Microsoft's Visual C++ compiler output.  
    It is featured here as the primary solution.

    On December 6, 2007, Hai Jin complained that the result was signed, so when
    computing the abs of the most negative value, it was still negative.

    On April 15, 2008 Andrew Shapira pointed out that the obvious approach
    could overflow, as it lacked an (unsigned) cast then;
    for maximum portability he suggested(v < 0) ? (1 + ((unsigned)(-1-v))) : (unsigned)v.  

    But citing the ISO C99 spec on July 9, 2008,
    Vincent Lefèvre convinced me
    to remove it becasue even on non-2s-complement machines -(unsigned)v
    will do the right thing.  The evaluation of -(unsigned)v first converts
    the negative value of v to an unsigned by adding 2**N,
    yielding a 2s complement representation of v's value that I'll call U.  
    Then, U is negated, giving the desired result,
    -U = 0 – U = 2**N – U = 2**N – (v+2**N) = -v = abs(v).


    Compute the minimum (min) or maximum (max) of two integers without branching

    int x;  // we want to find the minimum of x and y
    int y;   
    int r;  // the result goes here 
    
    r = y ^ ((x ^ y) & -(x < y)); // min(x, y)

    On some rare machines where branching is very expensive and no condition
    move instructions exist, the above expression
    might be faster than the obvious approach, r = (x < y) ? x : y, even
    though it involves two more instructions.  
    (Typically, the obvious approach is best, though.)
    It works because if x < y,
    then -(x < y) will be all ones,
    so r = y ^ (x ^ y) & ~0 = y ^ x ^ y = x.  
    Otherwise, if x >= y,
    then -(x < y) will be all zeros,
    so r = y ^ ((x ^ y) & 0) = y.  
    On some machines, evaluating (x < y) as 0
    or 1 requires a branch instruction, so there may be no
    advantage.

    To find the maximum, use:

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

    Quick and dirty versions:

    If you know that INT_MIN <= x – y <= INT_MAX,
    then you can use the following, which
    are faster because (x – y) only needs to be evaluated once.

    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)

    Note that the 1989 ANSI C specification doesn't specify the result of
    signed right-shift, so these aren't portable.
    If exceptions are thrown on overflows, then the values of x
    and y should be unsigned or cast to unsigned for the subtractions
    to avoid unnecessarily throwing an exception, however the right-shift
    needs a signed operand to produce all one bits when negative, so cast
    to signed there.

    On March 7, 2003, Angus Duggan pointed out the right-shift portability issue.
    On May 3, 2005, Randal E. Bryant alerted me to the need for the
    precondition, INT_MIN <= x – y <= INT_MAX,
    and suggested the non-quick and dirty version as a fix.  
    Both of these issues concern only the quick and dirty version.  
    Nigel Horspoon observed on July 6, 2005 that gcc produced the
    same code on a Pentium as the obvious solution because of how it
    evaluates (x < y).  On July 9, 2008 Vincent Lefèvre pointed out
    the potential for overflow exceptions with subtractions in
    r = y + ((x – y) & -(x < y)), which was the previous version.
    Timothy B. Terriberry suggested using xor rather than add and subract
    to avoid casting and the risk of overflows on June 2, 2009.


    Determining if an integer is a power of 2

    unsigned int v; // we want to see if v is a power of 2
    bool f;         // the result goes here 
    
    f = (v & (v - 1)) == 0;

    Note that 0 is incorrectly considered a power of 2 here.  To remedy
    this, use:

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

    Sign extending from a constant bit-width

    Sign extension is automatic for built-in types, such as chars and ints.
    But suppose you have a signed two's complement number, x, that is stored
    using only b bits.  Moreover, suppose you want to convert x to an int,
    which has more than b bits.  A simple copy will work if x
    is positive, but if negative, the sign must be extended.  For example,
    if we have only 4 bits to store a number, then -3 is represented as 1101
    in binary.  If we have 8 bits, then -3 is 11111101.  The most-significant
    bit of the 4-bit representation is replicated sinistrally to fill
    in the destination when we convert to a representation with more bits;
    this is sign extending.  
    In C, sign extension from a constant bit-width is trivial, since bit
    fields may be specified in structs or unions.  
    For example, to convert from 5 bits to an full integer:

    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;

    The following is a C++ template function that uses the same
    language feature to convert from B bits in one operation (though
    the compiler is generating more, of course).

    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

    John Byrd caught a typo in the code (attributed to html formatting)
    on May 2, 2005.  On March 4, 2006, Pat Wood pointed out that the ANSI C
    standard requires that the bitfield have the keyword "signed" to be signed;
    otherwise, the sign is undefined.


    Sign extending from a variable bit-width

    Sometimes we need to extend the sign of a number but we don't know a priori
    the number of bits, b, in which it is represented.  (Or we could be
    programming in a language like Java, which lacks bitfields.)

    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;

    The code above requires four operations, but when the bitwidth is a
    constant rather than variable, it requires only two fast operations,
    assuming the upper bits are already zeroes.

    A slightly faster but less portable method that doesn't depend on
    the bits in x above position b being zero is:

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

    Sean A. Irvine suggested that I
    add sign extension methods to this page on June 13, 2004, and he providedm = (1 << (b - 1)) - 1; r = -(x & ~m) | x;as a starting point from which I optimized to get
    m = 1U << (b – 1); r = -(x & m) | x.  
    But then on May 11, 2007, Shay Green suggested the version above,
    which requires one less operation than mine.  Vipin Sharma suggested
    I add a step to deal with situations where x had possible ones in bits
    other than the b bits we wanted to sign-extend on Oct. 15, 2008.
    On December 31, 2009 Chris Pirazzi suggested I add the faster version,
    which requires two operations for constant bit-widths and three
    for variable widths.


    Sign extending from a variable bit-width in 3 operations

    The following may be slow on some machines, due to the effort required for
    multiplication and division.  This version is 4 operations.  If you
    know that your initial bit-width, b, is greater than 1, you might do this
    type of sign extension in 3 operations by using
    r = (x * multipliers[b]) / multipliers[b],
    which requires only one array lookup.

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

    The following variation is not portable,
    but on architectures that employ an arithmetic right-shift,
    maintaining the sign, it should be fast.

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

    Randal E. Bryant pointed out a bug on May 3, 2005 in an earlier version
    (that used multipliers[] for divisors[]), where it failed on the case of
    x=1 and b=1.


    Conditionally set or clear bits without branching

    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;
    
    // OR, for superscalar CPUs:
    w = (w & ~m) | (-f & m);

    On some architectures, the lack of branching can more than make up for
    what appears to be twice as many operations.  For instance, informal
    speed tests on an AMD Athlon™ XP 2100+ indicated it was 5-10%
    faster.  An Intel Core 2 Duo ran the superscalar version about 16%
    faster than the first.
    Glenn Slayden informed me of the first expression on
    December 11, 2003.  Marco Yu shared the superscalar version with me on
    April 3, 2007 and alerted me to a typo 2 days later.


    Conditionally negate a value without branching

    If you need to negate only when a flag is false, then use the following
    to avoid branching:

    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;

    If you need to negate only when a flag is true, then use this:

    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;

    Avraham Plotnitzky suggested I add the first version on June 2, 2009.  
    Motivated to avoid the multiply, I came up with the second version on
    June 8, 2009.  Alfonso De Gregorio pointed out that some parens were
    missing on November 26, 2009, and received a bug bounty.


    Merge bits from two values according to a mask

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

    This shaves one operation from the obvious way of combining two sets of bits
    according to a bit mask.  If the mask is a constant, then there may be no
    advantage.

    Ron Jeffery sent this to me on February 9, 2006.


    Counting bits set (naive way)

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

    The naive approach requires one iteration per bit, until no more
    bits are set.  So on a 32-bit word with only the high set,
    it will go through 32 iterations.


    Counting bits set by lookup table

    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];
    }

    On July 14, 2009 Hallvard Furuseth suggested the macro compacted table.


    Counting bits set, Brian Kernighan's way

    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's method goes through as many iterations
    as there are set bits.  So if we have a 32-bit word with only
    the high bit set, then it will only go once through the loop.

    Published in 1988, the C Programming Language 2nd Ed.
    (by Brian W. Kernighan and Dennis M. Ritchie)
    mentions this in exercise 2-9.
    On April 19, 2006 Don Knuth pointed out to me that this method
    "was first published by Peter Wegner in CACM 3 (1960), 322.
    (Also discovered independently by Derrick Lehmer and published
    in 1964 in a book edited by Beckenbach.)"


    Counting bits set in 14, 24, or 32-bit words using 64-bit instructions

    unsigned int v; // count the number of bits set in v
    unsigned int c; // c accumulates the total bits set in v
    
    // option 1, for at most 14-bit values in v:
    c = (v * 0x200040008001ULL & 0x111111111111111ULL) % 0xf;
    
    // option 2, for at most 24-bit values in v:
    c =  ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
    c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) 
         % 0x1f;
    
    // option 3, for at most 32-bit values in v:
    c =  ((v & 0xfff) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;
    c += (((v & 0xfff000) >> 12) * 0x1001001001001ULL & 0x84210842108421ULL) % 
         0x1f;
    c += ((v >> 24) * 0x1001001001001ULL & 0x84210842108421ULL) % 0x1f;

    This method requires a 64-bit CPU with fast modulus division to be efficient.
    The first option takes only 3 operations; the second option takes 10; and
    the third option takes 15.

    Rich Schroeppel originally created a 9-bit version, similiar to option 1;
    see the Programming Hacks section ofBeeler, M., Gosper, R. W., and Schroeppel, R.
    HAKMEM. MIT AI Memo 239, Feb. 29, 1972.
    His method was the inspiration for the variants above,
    devised by Sean Anderson.  Randal E. Bryant offered a couple bug fixes on
    May 3, 2005.  Bruce Dawson tweaked what had been a 12-bit version and made
    it suitable for 14 bits using the same number of operations on
    Feburary 1, 2007.


    Counting bits set, in parallel

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

    The B array, expressed as binary, is:

    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

    We can adjust the method for larger integer sizes by continuing with
    the patterns for the Binary Magic Numbers, B and S.  

    If there are k bits, then we need the arrays S and B to be ceil(lg(k))
    elements long, and we must compute the same number of expressions for c as
    S or B are long.  For a 32-bit v, 16 operations are used.

    The best method for counting bits in a 32-bit integer v is the following:

    v = v - ((v >> 1) & 0x55555555);                    // reuse input as temporary
    v = (v & 0x33333333) + ((v >> 2) & 0x33333333);     // temp
    c = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24; // count

    The best bit counting method takes only 12 operations,
    which is the same as the lookup-table method,
    but avoids the memory and potential cache misses of a table.  
    It is a hybrid between the purely parallel method above
    and the earlier methods using multiplies (in the section on counting bits
    with 64-bit instructions), though it doesn't use 64-bit instructions.
    The counts of bits set in the bytes is done in parallel, and the sum
    total of the bits set in the bytes is computed by multiplying by 0x1010101
    and shifting right 24 bits.

    A generalization of the best bit counting method to integers
    of bit-widths upto 128 (parameterized by type T) is this:

    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

    See Ian Ashdown's nice newsgroup post for more information
    on counting the number of bits set (also known as sideways addition).
    The best bit counting method was brought to my attention on October 5, 2005
    by Andrew Shapira;
    he found it in pages 187-188 ofSoftware
    Optimization Guide for AMD Athlon™ 64 and Opteron™ Processors
    .
    Charlie Gordon suggested a way to shave off one operation from the
    purely parallel version on December 14, 2005, and Don Clugston trimmed three
    more from it on December 30, 2005.  I made a typo with Don's suggestion that
    Eric Cole spotted on January 8, 2006.  Eric later suggested the
    arbitrary bit-width generalization to the best method on November 17, 2006.
    On April 5, 2007, Al Williams observed that I had a line of dead code at the
    top of the first method.


    Count bits set (rank) from the most-significant bit upto a given position

    The following finds the the rank of a bit, meaning it returns the sum of
    bits that are set to 1 from the most-signficant bit downto the bit at
    the given position.

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

    Juha Järvi sent this to me on November 21, 2009 as an inverse operation
    to the computing the bit position with the given rank, which follows.


    Select the bit position (from the most-significant bit) with the given count (rank)

    The following 64-bit code selects the position of the rth 1 bit
    when counting from the left.  In other words if we start
    at the most significant bit and proceed to the right,
    counting the number of bits set to 1 until we reach the desired rank, r,
    then the position where we stop is returned.  If the rank requested exceeds
    the count of bits set, then 64 is returned.
    The code may be modified for 32-bit or counting from the right.

      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;

    If branching is fast on your target CPU, consider uncommenting the
    if-statements and commenting the lines that follow them.

    Juha Järvi sent this to me on November 21, 2009.


    Computing parity the naive way

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

    The above code uses an approach like Brian Kernigan's bit counting, above.  
    The time it takes is proportional to the number of bits set.


    Compute parity by lookup table

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

    Randal E. Bryant encouraged the addition of the (admittedly) obvious
    last variation with variable p on May 3, 2005.  Bruce Rawles found a typo
    in an instance of the table variable's name on September 27, 2005, and
    he received a $10 bug bounty.  On October 9, 2006, Fabrice Bellard
    suggested the 32-bit variations above, which require only one table lookup;
    the previous version had four lookups (one per byte) and were slower.
    On July 14, 2009 Hallvard Furuseth suggested the macro compacted table.


    Compute parity of a byte using 64-bit multiply and modulus division

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

    The method above takes around 4 operations, but only works on bytes.


    Compute parity of word with a multiply

    The following method computes the parity of the 32-bit value in only 8
    operations using a multiply.

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

    Also for 64-bits, 8 operations are still enough.

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

    Andrew Shapira came up with this and sent it to me on Sept. 2, 2007.


    Compute parity in parallel

    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;

    The method above takes around 9 operations, and works for 32-bit words.  
    It may be optimized to work just on bytes in 5 operations by removing the
    two lines immediately following "unsigned int v;".  The method first shifts
    and XORs the eight nibbles of the 32-bit value together, leaving the result
    in the lowest nibble of v.  Next, the binary number 0110 1001 1001 0110
    (0x6996 in hex) is shifted to the right by the value represented in the lowest
    nibble of v.  This number is like a miniature 16-bit parity-table indexed by
    the low four bits in v.  The result has the parity of v in bit 1, which is
    masked and returned.

    Thanks to Mathew Hendry for pointing out the shift-lookup idea at the end on
    Dec. 15, 2002.  
    That optimization shaves two operations off using only shifting and XORing
    to find the parity.


    Swapping values with subtraction and addition

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

    This swaps the values of a and b without using a temporary variable.The initial check for a and b being the same location in memory may be
    omitted when you know this can't happen.  (The compiler may omit it anyway
    as an optimization.)  If you enable overflows exceptions, then pass unsigned
    values so an exception isn't thrown.

    The XOR method that follows may be slightly faster on some machines.
    Don't use this with floating-point numbers (unless you operate on
    their raw integer representations).

    Sanjeev Sivasankaran suggested I add this on June 12, 2007.  
    Vincent Lefèvre pointed out the potential for overflow exceptions on
    July 9, 2008


    Swapping values with XOR

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

    This is an old trick to exchange the values of the variables a and bwithout using extra space for a temporary variable.

    On January 20, 2005, Iain A. Fleming pointed out that the macro above
    doesn't work when you swap with the same memory location, such as
    SWAP(a[i], a[j]) with i == j.  So if that may occur, consider
    defining the macro as
    (((a) == (b)) || (((a) ^= (b)), ((b) ^= (a)), ((a) ^= (b)))).
    On July 14, 2009, Hallvard Furuseth suggested that on some machines,
    (((a) ^ (b)) && ((b) ^= (a) ^= (b), (a) ^= (b))) might be faster,
    since the (a) ^ (b) expression is reused.


    Swapping individual bits with 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));

    As an example of swapping ranges of bits
    suppose we have have b = 00101111
    (expressed in binary) and we want to swap the n = 3 consecutive bits
    starting at i = 1 (the second bit from the right) with the 3 consecutive
    bits starting at j = 5; the result would be r =11100011 (binary).

    This method of swapping is similar to the general purpose XOR swap
    trick, but intended for operating on individual bits.   The variable x
    stores the result of XORing the pairs of bit values we want to swap,
    and then the bits are set to the result of themselves XORed with x.   
    Of course, the result is undefined if the sequences overlap.

    On July 14, 2009 Hallvard Furuseth suggested that I change the
    1 << n to 1U << n
    because the value was being assigned to an unsigned and to avoid shifting into
    a sign bit.


    Reverse bits the obvious way

    unsigned int v;     // input bits to be reversed
    unsigned int r = v; // r will be reversed bits of v; first get LSB of v
    int s = sizeof(v) * CHAR_BIT - 1; // extra shift needed at end
    
    for (v >>= 1; v; v >>= 1)
    {   
      r <<= 1;
      r |= v & 1;
      s--;
    }
    r <<= s; // shift when v's highest bits are zero

    On October 15, 2004, Michael Hoisie pointed out a bug in the original version.
    Randal E. Bryant suggested removing an extra operation on May 3, 2005.  
    Behdad Esfabod suggested a slight change that eliminated one iteration of the
    loop on May 18, 2005.  Then, on February 6, 2007, Liyong Zhou suggested a
    better version that loops while v is not 0, so rather than iterating over
    all bits it stops early.


    Reverse bits in word by lookup table

    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; // reverse 32-bit value, 8 bits at time
    unsigned int c; // c will get v reversed
    
    // 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]];

    The first method takes about 17 operations, and the second takes about 12,
    assuming your CPU can load and store bytes easily.

    On July 14, 2009 Hallvard Furuseth suggested the macro compacted table.


    Reverse the bits in a byte with 3 operations (64-bit multiply and modulus division):

    unsigned char b; // reverse this (8-bit) byte
     
    b = (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;

    The multiply operation creates five separate copies of the 8-bit
    byte pattern to fan-out into a 64-bit value.
    The AND operation selects the bits that are in the correct (reversed)
    positions, relative to each 10-bit groups of bits.
    The multiply and the AND operations copy the bits from the original
    byte so they each appear in only one of the 10-bit sets.  
    The reversed positions of the bits from the original byte coincide with
    their relative positions within any 10-bit set.

    The last step, which involves modulus division by 2^10 – 1, has the
    effect of merging together each set of 10 bits
    (from positions 0-9, 10-19, 20-29, …) in the 64-bit value.  
    They do not overlap, so the addition steps underlying the
    modulus division behave like or operations.

    This method was attributed to Rich Schroeppel in the
    Programming Hacks section ofBeeler, M., Gosper, R. W., and Schroeppel, R.
    HAKMEM. MIT AI Memo 239, Feb. 29, 1972.


    Reverse the bits in a byte with 4 operations (64-bit multiply, no division):

    unsigned char b; // reverse this byte
     
    b = ((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32;

    The following shows the flow of the bit values with the boolean variablesa, b, c, d, e, f, g, and h, which
    comprise an 8-bit byte.  Notice how the first multiply fans out the
    bit pattern to multiple copies, while the last multiply combines them
    in the fifth byte from the right.

                                                                                            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

    Note that the last two steps can be combined on some processors because
    the registers can be accessed as bytes;
    just multiply so that a register stores the upper 32 bits of the result
    and the take the low byte.  Thus, it may take only 6 operations.

    Devised by Sean Anderson, July 13, 2001.


    Reverse the bits in a byte with 7 operations (no 64-bit):

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

    Make sure you assign or cast the result to an unsigned char to remove
    garbage in the higher bits.  
    Devised by Sean Anderson, July 13, 2001.  Typo spotted and correction supplied
    by Mike Keith, January 3, 2002.


    Reverse an N-bit quantity in parallel in 5 * lg(N) operations:

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

    The following variation is also O(lg(N)), however it requires more
    operations to reverse v.  Its virtue is in taking less slightly memory
    by computing the constants on the fly.

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

    These methods above are best suited to situations where N is large.
    If you use the above with 64-bit ints (or larger), then you need
    to add more lines (following the pattern); otherwise only the lower
    32 bits will be reversed and the result will be in the lower 32 bits.

    See Dr. Dobb's Journal 1983, Edwin Freed's article on Binary Magic
    Numbers for more information.  The second variation was suggested
    by Ken Raeburn on September 13, 2005.  Veldmeijer mentioned that
    the first version could do without ANDS in the last line
    on March 19, 2006.


    Compute modulus division by 1 << s without a division operator

    const unsigned int n;          // numerator
    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);

    Most programmers learn this trick early, but it was included for the
    sake of completeness.


    Compute modulus division by (1 << s) – 1 without a division operator

    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.
    
    for (m = n; n > d; n = m)
    {
      for (m = 0; n; n >>= s)
      {
        m += n & d;
      }
    }
    // Now m is a value from 0 to d, but since with modulus division
    // we want m to be 0 when it is d.
    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 assignm = ((m + 1) & d) - 1; at the end.  Michael Miller spotted a
    typo in the code April 25, 2005.


    Compute modulus division by (1 << s) – 1 in parallel without a division
    operator

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

    This method of finding modulus division by an integer that is one less
    than a power of 2 takes at most O(lg(N)) time, where N is the number of
    bits in the numerator (32 bits, for the code above).  
    The number of operations is at most 12 + 9 * ceil(lg(N)).

    The tables may be removed if you know
    the denominator at compile time; just extract the few relevent entries
    and unroll the loop.  It may be easily extended to more bits.

    It finds the result by summing the values in base (1 << s) in parallel.
    First every other base (1 << s) value is added to the previous one.
    Imagine that the result is written on a piece of paper.  Cut the paper
    in half, so that half the values are on each cut piece.  Align the values
    and sum them onto a new piece of paper.  Repeat by cutting this paper
    in half (which will be a quarter of the size of the previous one)
    and summing, until you cannot cut further.  After performing lg(N/s/2)
    cuts, we cut no more; just continue to add the values and put the result
    onto a new piece of paper as before, while there are at least two s-bit
    values.

    Devised by Sean Anderson, August 20, 2001.  A typo was spotted by
    Randy E. Bryant on May 3, 2005 (after pasting the code, I had later
    added "unsinged" to a variable declaration).  As in the previous hack,
    I mistakenly commented that we could alternatively assignm = ((m + 1) & d) - 1; at the end, and Don Knuth corrected
    me on April 19, 2006 and suggestedm = m & -((signed)(m - d) >> s).  
    On June 18, 2009 Sean Irvine proposed a change that used((n >> s) & M[s]) instead of((n & ~M[s]) >> s),
    which typically requires fewer operations because the M[s] constant is already
    loaded.


    Find the log base 2 of an integer with the MSB N set in O(N) operations
    (the obvious way)

    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++;
    }

    The log base 2 of an integer is the same as the position of the highest
    bit set (or most significant bit set, MSB).  The following log base 2
    methods are faster than this one.


    Find the integer log base 2 of an integer with an 64-bit IEEE float

    int v; // 32-bit integer to find the log base 2 of
    int r; // result of log_2(v) goes here
    union { unsigned int u[2]; double d; } t; // temp
    
    t.u[__FLOAT_WORD_ORDER==LITTLE_ENDIAN] = 0x43300000;
    t.u[__FLOAT_WORD_ORDER!=LITTLE_ENDIAN] = v;
    t.d -= 4503599627370496.0;
    r = (t.u[__FLOAT_WORD_ORDER==LITTLE_ENDIAN] >> 20) - 0x3FF;

    The code above loads a 64-bit (IEEE-754 floating-point) double with
    a 32-bit integer (with no paddding bits) by storing the integer in the
    mantissa while the exponent is set to 252.  
    From this newly minted double,
    252 (expressed as a double) is subtracted, which sets the
    resulting exponent to the log base 2 of the input value, v.  All that is
    left is shifting the exponent bits into position (20 bits right) and
    subtracting the bias, 0x3FF (which is 1023 decimal).  This technique
    only takes 5 operations, but many CPUs are slow at manipulating doubles,
    and the endianess of the architecture must be accommodated.

    Eric Cole sent me this on January 15, 2006.  Evan Felix pointed out a typo
    on April 4, 2006.  Vincent Lefèvre told me on July 9, 2008 to
    change the endian check to use the float's endian, which could differ
    from the integer's endian.


    Find the log base 2 of an integer with a lookup table

    static const char LogTable256[256] = 
    {
    #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
        -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
        LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
        LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
    };
    
    unsigned int v; // 32-bit word to find the log of
    unsigned r;     // r will be lg(v)
    register unsigned int t, tt; // temporaries
    
    if (tt = v >> 16)
    {
      r = (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
    }
    else 
    {
      r = (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
    }

    The lookup table method takes only about 7 operations to find the log of
    a 32-bit value.  If extended for 64-bit quantities, it would take roughly
    9 operations.  Another operation can be trimmed off by using four tables,
    with the possible additions incorporated into each.  Using int table
    elements may be faster, depending on your architecture.

    The code above is tuned to uniformly distributed output values.  
    If your inputs are evenly distributed across all 32-bit values,
    then consider using the following:

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

    To initially generate the log table algorithmically:

    LogTable256[0] = LogTable256[1] = 0;
    for (int i = 2; i < 256; i++) 
    {
      LogTable256[i] = 1 + LogTable256[i / 2];
    }
    LogTable256[0] = -1; // if you want log(0) to return -1

    Behdad Esfahbod and I shaved off a fraction of an operation (on average) on
    May 18, 2005.  Yet another fraction of an operation was removed on November
    14, 2006 by Emanuel Hoogeveen.  The variation that is tuned to evenly
    distributed input values was suggested by David A. Butterfield on September
    19, 2008.  Venkat Reddy told me on January 5, 2009 that log(0) should return
    -1 to indicate an error, so I changed the first entry in the table to that.


    Find the log base 2 of an N-bit integer in O(lg(N)) operations

    unsigned int v;  // 32-bit value to find the log2 of 
    const unsigned int b[] = {0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000};
    const unsigned int S[] = {1, 2, 4, 8, 16};
    int i;
    
    register unsigned int r = 0; // result of log2(v) will go here
    for (i = 4; i >= 0; i--) // unroll for speed...
    {
      if (v & b[i])
      {
        v >>= S[i];
        r |= S[i];
      } 
    }
    
    
    // OR (IF YOUR CPU BRANCHES SLOWLY):
    
    unsigned int v;	         // 32-bit value to find the log2 of 
    register unsigned int r; // result of log2(v) will go here
    register unsigned int shift;
    
    r =     (v > 0xFFFF) << 4; v >>= r;
    shift = (v > 0xFF  ) << 3; v >>= shift; r |= shift;
    shift = (v > 0xF   ) << 2; v >>= shift; r |= shift;
    shift = (v > 0x3   ) << 1; v >>= shift; r |= shift;
                                            r |= (v >> 1);
    
    
    // OR (IF YOU KNOW v IS A POWER OF 2):
    
    unsigned int v;  // 32-bit value to find the log2 of 
    static const unsigned int b[] = {0xAAAAAAAA, 0xCCCCCCCC, 0xF0F0F0F0, 
                                     0xFF00FF00, 0xFFFF0000};
    register unsigned int r = (v & b[0]) != 0;
    for (i = 4; i > 0; i--) // unroll for speed...
    {
      r |= ((v & b[i]) != 0) << i;
    }

    Of course, to extend the code to find the log of a 33- to 64-bit
    number, we would append another element, 0xFFFFFFFF00000000, to b,
    append 32 to S, and loop from 5 to 0.  This method is much slower
    than the earlier table-lookup