---
created_at: '2017-11-19T14:55:13.000Z'
title: Bit Twiddling Hacks (2005)
url: https://graphics.stanford.edu/~seander/bithacks.html
author: kercker
points: 66
story_text: 
comment_text: 
num_comments: 3
story_id: 
story_title: 
story_url: 
parent_id: 
created_at_i: 1511103313
_tags:
- story
- author_kercker
- story_15734119
objectID: '15734119'

---
## 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](http://www-2.cs.cmu.edu/~bryant/), the Dean of Computer Science
at Carnegie Mellon University, has personally tested almost everything
with his [Uclid code verification
system](http://www-2.cs.cmu.edu/~uclid/). 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

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.

    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](http://rpg-314.blogspot.com/)
suggested the version for non-negative integers on September 12, 2009.

    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.

    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](http://patft.uspto.gov/netacgi/nph-Parser?Sect1=PTO2&Sect2=HITOFF&p=1&u=/netahtml/search-adv.htm&r=1&f=G&l=50&d=ptxt&S1=6073150&OS=6073150&RS=6073150)
in the USA on June 6, 2000 by Vladimir Yu Volkonsky and assigned to [Sun
Microsystems](http://www.sun.com/). 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 in [How to
Optimize for the Pentium
Processor](http://www.goof.com/pcg/doc/pentopt.txt) 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 old
[link](http://web.archive.org/web/19961201174141/www.x86.org/ftp/articles/pentopt/PENTOPT.TXT)
to it. On January 30, 2007, Peter Kankowski shared with me an [abs
version](http://smallcode.weblogs.us/2007/01/31/microsoft-probably-uses-the-abs-function-patented-by-sun/)
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).

    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  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.

    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 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.

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 provided `m = (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.

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.

    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.

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.

    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.

    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.

    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.

    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.)"

    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 of [Beeler, M., Gosper, R. W., and
Schroeppel, R. HAKMEM. MIT AI Memo 239, Feb. 29,
1972.](http://www.inwap.com/pdp10/hbaker/hakmem/hakmem.html) 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.

    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](http://groups.google.com/groups?q=reverse+bits&num=100&hl=en&group=comp.graphics.algorithms&imgsafe=off&safe=off&rnum=2&ic=1&selm=4fulhm%248dn%40atlas.uniserve.com)
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](http://onezero.org); he
found it in pages 187-188 of [Software Optimization Guide for AMD
Athlon™ 64 and Opteron™
Processors](http://www.amd.com/us-en/assets/content_type/white_papers_and_tech_docs/25112.PDF).
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.

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.

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.

    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.

    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.

    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.

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.

    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.

    #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

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

This is an old trick to exchange the values of the variables a and b
without 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.

    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 =
**001**0**111**1 (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 =
**111**0**001**1 (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.

    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.

    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.

    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 of [Beeler, M., Gosper, R. W., and Schroeppel, R. HAKMEM. MIT AI
Memo 239, Feb. 29,
1972.](http://www.inwap.com/pdp10/hbaker/hakmem/hakmem.html)

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

The following shows the flow of the bit values with the boolean
variables `a, 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.

    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.

    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.

    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.

    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 assign `m = ((m + 1) & d) - 1;` at the end. Michael Miller
spotted a typo in the code April 25, 2005.

``` 

// 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 assign `m = ((m + 1) &
d) - 1;` at the end, and Don Knuth corrected me on April 19, 2006 and
suggested `m = 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.

    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.

    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.

    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.

    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 version, but if you don't want big table or your
architecture is slow to access memory, it's a good choice. The second
variation involves slightly more operations, but it may be faster on
machines with high branch costs (e.g. PowerPC).

The second version was sent to me by [Eric
Cole](http://www.balance-software.com/ec/) on January 7, 2006. Andrew
Shapira subsequently trimmed a few operations off of it and sent me his
variation (above) on Sept. 1, 2007. The third variation was suggested to
me by [John Owens](http://www.ece.ucdavis.edu/~jowens/) on April 24,
2002; it's faster, but it is only suitable when the input is known to be
a power of 2. On May 25, 2003, Ken Raeburn suggested improving the
general case by using smaller numbers for b\[\], which load faster on
some architectures (for instance if the word size is 16 bits, then only
one load instruction may be needed). These values work for the general
version, but not for the special-case version below it, where v is a
power of 2; Glenn Slayden brought this oversight to my attention on
December 12, 2003.

    uint32_t v; // find the log base 2 of 32-bit v
    int r;      // result goes here
    
    static const int MultiplyDeBruijnBitPosition[32] = 
    {
      0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
      8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
    };
    
    v |= v >> 1; // first round down to one less than a power of 2 
    v |= v >> 2;
    v |= v >> 4;
    v |= v >> 8;
    v |= v >> 16;
    
    r = MultiplyDeBruijnBitPosition[(uint32_t)(v * 0x07C4ACDDU) >> 27];

The code above computes the log base 2 of a 32-bit integer with a small
table lookup and multiply. It requires only 13 operations, compared to
(up to) 20 for the previous method. The purely table-based method
requires the fewest operations, but this offers a reasonable compromise
between table size and speed.

If you know that v is a power of 2, then you only need the following:

    static const int MultiplyDeBruijnBitPosition2[32] = 
    {
      0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 
      31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
    };
    r = MultiplyDeBruijnBitPosition2[(uint32_t)(v * 0x077CB531U) >> 27];

Eric Cole devised this January 8, 2006 after reading about the entry
below to [round up to a power of 2](#RoundUpPowerOf2) and the method
below for [computing the number of trailing bits with a multiply and
lookup](#ZerosOnRightMultLookup) using a DeBruijn sequence. On December
10, 2009, Mark Dickinson shaved off a couple operations by requiring v
be rounded up to one less than the next power of 2 rather than the power
of
    2.

    unsigned int v; // non-zero 32-bit integer value to compute the log base 10 of 
    int r;          // result goes here
    int t;          // temporary
    
    static unsigned int const PowersOf10[] = 
        {1, 10, 100, 1000, 10000, 100000,
         1000000, 10000000, 100000000, 1000000000};
    
    t = (IntegerLogBase2(v) + 1) * 1233 >> 12; // (use a lg2 method from above)
    r = t - (v < PowersOf10[t]);

The integer log base 10 is computed by first using one of the techniques
above for finding the log base 2. By the relationship log10(v) = log2(v)
/ log2(10), we need to multiply it by 1/log2(10), which is approximately
1233/4096, or 1233 followed by a right shift of 12. Adding one is needed
because the IntegerLogBase2 rounds down. Finally, since the value t is
only an approximation that may be off by one, the exact value is found
by subtracting the result of v \< PowersOf10\[t\].

This method takes 6 more operations than IntegerLogBase2. It may be sped
up (on machines with fast memory access) by modifying the log base 2
table-lookup method above so that the entries hold what is computed for
t (that is, pre-add, -mulitply, and -shift). Doing so would require a
total of only 9 operations to find the log base 10, assuming 4 tables
were used (one for each byte of v).

Eric Cole suggested I add a version of this on January 7,
    2006.

    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;

This method works well when the input is uniformly distributed over
32-bit values because 76% of the inputs are caught by the first compare,
21% are caught by the second compare, 2% are caught by the third, and so
on (chopping the remaining down by 90% with each comparision). As a
result, less than 2.6 operations are needed on average.

On April 18, 2007, Emanuel Hoogeveen suggested a variation on this where
the conditions used divisions, which were not as fast as simple
comparisons.

    const float v; // find int(log2(v)), where v > 0.0 && finite(v) && isnormal(v)
    int c;         // 32-bit int c gets the result;
    
    c = *(const int *) &v;  // OR, for portability:  memcpy(&c, &v, sizeof c);
    c = (c >> 23) - 127;

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:

    const float v;              // find int(log2(v)), where v > 0.0 && finite(v)
    int c;                      // 32-bit int c gets the result;
    int x = *(const int *) &v;  // OR, for portability:  memcpy(&x, &v, sizeof x);
    
    c = x >> 23;          
    
    if (c)
    {
      c -= 127;
    }
    else
    { // subnormal, so recompute using mantissa: c = intlog2(x) - 149;
      register unsigned int t; // temporary
      // Note that LogTable256 was defined earlier
      if (t = x >> 16)
      {
        c = LogTable256[t] - 133;
      }
      else
      {
        c = (t = x >> 8) ? LogTable256[t] - 141 : LogTable256[x] - 149;
      }
    }

On June 20, 2004, Sean A. Irvine suggested that I include code to handle
subnormal numbers. On June 11, 2005, Falk Hüffner pointed out that ISO
C99 6.5/7 specified undefined behavior for the common type punning idiom
\*(int \*)&, though it has worked on 99.9% of C compilers. He proposed
using memcpy for maximum portability or a union with a float and an int
for better code generation than memcpy on some compilers.

    const int r;
    const float v; // find int(log2(pow((double) v, 1. / pow(2, r)))), 
                   // where isnormal(v) and v > 0
    int c;         // 32-bit int c gets the result;
    
    c = *(const int *) &v;  // OR, for portability:  memcpy(&c, &v, sizeof c);
    c = ((((c - 0x3f800000) >> r) + 0x3f800000) >> 23) - 127;

So, if r is 0, for example, we have c = int(log2((double) v)). If r is
1, then we have c = int(log2(sqrt((double) v))). If r is 2, then we have
c = int(log2(pow((double) v, 1./4))).

On June 11, 2005, Falk Hüffner pointed out that ISO C99 6.5/7 left the
type punning idiom \*(int \*)& undefined, and he suggested using memcpy.

    unsigned int v;  // input to count trailing zero bits
    int c;  // output: c will count v's trailing zero bits,
            // so if v is 1101000 (base 2), then c will be 3
    if (v)
    {
      v = (v ^ (v - 1)) >> 1;  // Set v's trailing 0s to 1s and zero rest
      for (c = 0; v; c++)
      {
        v >>= 1;
      }
    }
    else
    {
      c = CHAR_BIT * sizeof(v);
    }

The average number of trailing zero bits in a (uniformly distributed)
random binary number is one, so this O(trailing zeros) solution isn't
that bad compared to the faster methods below.

Jim Cole suggested I add a linear-time method for counting the trailing
zeros on August 15, 2007. On October 22, 2007, Jason Cunningham pointed
out that I had neglected to paste the unsigned modifier for
    v.

    unsigned int v;      // 32-bit word input to count zero bits on right
    unsigned int c = 32; // c will be the number of zero bits on the right
    v &= -signed(v);
    if (v) c--;
    if (v & 0x0000FFFF) c -= 16;
    if (v & 0x00FF00FF) c -= 8;
    if (v & 0x0F0F0F0F) c -= 4;
    if (v & 0x33333333) c -= 2;
    if (v & 0x55555555) c -= 1;

Here, we are basically doing the same operations as finding the log base
2 in parallel, but we first isolate the lowest 1 bit, and then proceed
with c starting at the maximum and decreasing. The number of operations
is at most 3 \* lg(N) + 4, roughly, for N bit words.

Bill Burdick suggested an optimization, reducing the time from 4 \*
lg(N) on February 4, 2011.

``` 
unsigned int v;     // 32-bit word input to count zero bits on right
unsigned int c;     // c will be the number of zero bits on the right,
                    // so if v is 1101000 (base 2), then c will be 3
// NOTE: if 0 == v, then c = 31.
if (v & 0x1) 
{
  // special case for odd v (assumed to happen half of the time)
  c = 0;
}
else
{
  c = 1;
  if ((v & 0xffff) == 0) 
  {  
    v >>= 16;  
    c += 16;
  }
  if ((v & 0xff) == 0) 
  {  
    v >>= 8;  
    c += 8;
  }
  if ((v & 0xf) == 0) 
  {  
    v >>= 4;
    c += 4;
  }
  if ((v & 0x3) == 0) 
  {  
    v >>= 2;
    c += 2;
  }
  c -= v & 0x1;
}   
```

The code above is similar to the previous method, but it computes the
number of trailing zeros by accumulating c in a manner akin to binary
search. In the first step, it checks if the bottom 16 bits of v are
zeros, and if so, shifts v right 16 bits and adds 16 to c, which reduces
the number of bits in v to consider by half. Each of the subsequent
conditional steps likewise halves the number of bits until there is only
1. This method is faster than the last one (by about 33%) because the
bodies of the if statements are executed less often.

Matt Whitlock suggested this on January 25, 2006. Andrew Shapira shaved
a couple operations off on Sept. 5, 2007 (by setting c=1 and
unconditionally subtracting at the end).

    unsigned int v;            // find the number of trailing zeros in v
    int r;                     // the result goes here
    float f = (float)(v & -v); // cast the least significant bit in v to a float
    r = (*(uint32_t *)&f >> 23) - 0x7f;

Although this only takes about 6 operations, the time to convert an
integer to a float can be high on some machines. The exponent of the
32-bit IEEE floating point representation is shifted down, and the bias
is subtracted to give the position of the least significant 1 bit set in
v. If v is zero, then the result is -127.

    unsigned int v;  // find the number of trailing zeros in v
    int r;           // put the result in r
    static const int Mod37BitPosition[] = // map a bit value mod 37 to its position
    {
      32, 0, 1, 26, 2, 23, 27, 0, 3, 16, 24, 30, 28, 11, 0, 13, 4,
      7, 17, 0, 25, 22, 31, 15, 29, 10, 12, 6, 0, 21, 14, 9, 5,
      20, 8, 19, 18
    };
    r = Mod37BitPosition[(-v & v) % 37];

The code above finds the number of zeros that are trailing on the right,
so binary 0100 would produce 2. It makes use of the fact that the first
32 bit position values are relatively prime with 37, so performing a
modulus division with 37 gives a unique number from 0 to 36 for each.
These numbers may then be mapped to the number of zeros using a small
lookup table. It uses only 4 operations, however indexing into a table
and performing modulus division may make it unsuitable for some
situations. I came up with this independently and then searched for a
subsequence of the table values, and found it was invented earlier by
Reiser, according to [Hacker's
Delight](http://www.hackersdelight.org/HDcode/ntz.c.txt).

    unsigned int v;  // find the number of trailing zeros in 32-bit v 
    int r;           // result goes here
    static const int MultiplyDeBruijnBitPosition[32] = 
    {
      0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, 
      31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
    };
    r = MultiplyDeBruijnBitPosition[((uint32_t)((v & -v) * 0x077CB531U)) >> 27];

Converting bit vectors to indices of set bits is an example use for
this. It requires one more operation than the earlier one involving
modulus division, but the multiply may be faster. The expression (v &
-v) extracts the least significant 1 bit from v. The constant
0x077CB531UL is a de Bruijn sequence, which produces a unique pattern of
bits into the high 5 bits for each possible bit position that it is
multiplied against. When there are no bits set, it returns 0. More
information can be found by reading the paper [Using de Bruijn Sequences
to Index 1 in a Computer
Word](http://citeseer.ist.psu.edu/leiserson98using.html) by Charles E.
Leiserson, Harald Prokof, and Keith H. Randall.

On October 8, 2005 [Andrew Shapira](http://onezero.org) suggested I add
this. Dustin Spicuzza asked me on April 14, 2009 to cast the result of
the multiply to a 32-bit type so it would work when compiled with 64-bit
ints.

    unsigned int const v; // Round this 32-bit value to the next highest power of 2
    unsigned int r;       // Put the result here. (So v=3 -> r=4; v=8 -> r=8)
    
    if (v > 1) 
    {
      float f = (float)v;
      unsigned int const t = 1U << ((*(unsigned int *)&f >> 23) - 0x7f);
      r = t << (t < v);
    }
    else 
    {
      r = 1;
    }

The code above uses 8 operations, but works on all v \<= (1\<\<31).

Quick and dirty version, for domain of 1 \< v \< (1\<\<25):

    float f = (float)(v - 1);  
    r = 1U << ((*(unsigned int*)(&f) >> 23) - 126);

Although the quick and dirty version only uses around 6 operations, it
is roughly three times slower than the [technique
below](#RoundUpPowerOf2) (which involves 12 operations) when benchmarked
on an Athlon™ XP 2100+ CPU. Some CPUs will fare better with it, though.

On September 27, 2005 Andi Smithers suggested I include a technique for
casting to floats to find the lg of a number for rounding up to a power
of 2. Similar to the quick and dirty version here, his version worked
with values less than (1\<\<25), due to mantissa rounding, but it used
one more operation.

    unsigned int v; // compute the next highest power of 2 of 32-bit v
    
    v--;
    v |= v >> 1;
    v |= v >> 2;
    v |= v >> 4;
    v |= v >> 8;
    v |= v >> 16;
    v++;

In 12 operations, this code computes the next highest power of 2 for a
32-bit integer. The result may be expressed by the formula 1U \<\< (lg(v
- 1) + 1). Note that in the edge case where v is 0, it returns 0, which
isn't a power of 2; you might append the expression v += (v == 0) to
remedy this if it matters. It would be faster by 2 operations to use the
formula and the log base 2 method that uses a lookup table, but in some
situations, lookup tables are not suitable, so the above code may be
best. (On a Athlon™ XP 2100+ I've found the above shift-left and then OR
code is as fast as using a single BSR assembly language instruction,
which scans in reverse to find the highest set bit.) It works by copying
the highest set bit to all of the lower bits, and then adding one, which
results in carries that set all of the lower bits to 0 and one bit
beyond the highest set bit to 1. If the original number was a power of
2, then the decrement will reduce it to one less, so that we round up to
the same original value.

You might alternatively compute the next higher power of 2 in only 8 or
9 operations using a lookup table for floor(lg(v)) and then evaluating
1\<\<(1+floor(lg(v))); Atul Divekar suggested I mention this on
September 5, 2010.

Devised by Sean Anderson, Sepember 14, 2001. Pete Hart pointed me to [a
couple newsgroup
posts](http://groups.google.com/group/comp.lang.python/browse_thread/thread/c4d3aae0df917df5/6fdae3872f9de79d?lnk=st&q=comp.lang.python+zeddy&rnum=6#6fdae3872f9de79d)
by him and William Lewis in February of 1997, where they arrive at the
same
    algorithm.

    unsigned short x;   // Interleave bits of x and y, so that all of the
    unsigned short y;   // bits of x are in the even positions and y in the odd;
    unsigned int z = 0; // z gets the resulting Morton Number.
    
    for (int i = 0; i < sizeof(x) * CHAR_BIT; i++) // unroll for more speed...
    {
      z |= (x & 1U << i) << i | (y & 1U << i) << (i + 1);
    }

Interleaved bits (aka Morton numbers) are useful for linearizing 2D
integer coordinates, so x and y are combined into a single number that
can be compared easily and has the property that a number is usually
close to another if their x and y values are close.

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

For more speed, use an additional table with values that are
MortonTable256 pre-shifted one bit to the left. This second table could
then be used for the y lookups, thus reducing the operations by two, but
almost doubling the memory required. Extending this same idea, four
tables could be used, with two of them pre-shifted by 16 to the left of
the previous two, so that we would only need 11 operations total. In 11
operations, this version interleaves bits of two bytes (rather than
shorts, as in the other versions), but many of the operations are 64-bit
multiplies so it isn't appropriate for all machines. The input
parameters, x and y, should be less than
    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 was inspired to suggest this technique on October 10, 2004
after reading the multiply-based bit reversals
    here.

    static const unsigned int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF};
    static const unsigned int S[] = {1, 2, 4, 8};
    
    unsigned int x; // Interleave lower 16 bits of x and y, so the bits of x
    unsigned int y; // are in the even positions and bits from y in the odd;
    unsigned int z; // z gets the resulting 32-bit Morton Number.  
                    // x and y must initially be less than 65536.
    
    x = (x | (x << S[3])) & B[3];
    x = (x | (x << S[2])) & B[2];
    x = (x | (x << S[1])) & B[1];
    x = (x | (x << S[0])) & B[0];
    
    y = (y | (y << S[3])) & B[3];
    y = (y | (y << S[2])) & B[2];
    y = (y | (y << S[1])) & B[1];
    y = (y | (y << S[0])) & B[0];
    
    z = x | (y << 1);

    // Fewer operations:
    unsigned int v; // 32-bit word to check if any 8-bit byte in it is 0
    bool hasZeroByte = ~((((v & 0x7F7F7F7F) + 0x7F7F7F7F) | v) | 0x7F7F7F7F);

The code above may be useful when doing a fast string copy in which a
word is copied at a time; it uses 5 operations. On the other hand,
testing for a null byte in the obvious ways (which follow) have at least
7 operations (when counted in the most sparing way), and at most 12.

    // More operations:
    bool hasNoZeroByte = ((v & 0xff) && (v & 0xff00) && (v & 0xff0000) && (v & 0xff000000))
    // OR:
    unsigned char * p = (unsigned char *) &v;  
    bool hasNoZeroByte = *p && *(p + 1) && *(p + 2) && *(p + 3);

The code at the beginning of this section (labeled "Fewer operations")
works by first zeroing the high bits of the 4 bytes in the word.
Subsequently, it adds a number that will result in an overflow to the
high bit of a byte if any of the low bits were initialy set. Next the
high bits of the original word are ORed with these values; thus, the
high bit of a byte is set iff any bit in the byte was set. Finally, we
determine if any of these high bits are zero by ORing with ones
everywhere except the high bits and inverting the result. Extending to
64 bits is trivial; simply increase the constants to be
0x7F7F7F7F7F7F7F7F.

For an additional improvement, a fast pretest that requires only 4
operations may be performed to determine if the word may have a zero
byte. The test also returns true if the high byte is 0x80, so there are
occasional false positives, but the slower and more reliable version
above may then be used on candidates for an overall increase in speed
with correct output.

    bool hasZeroByte = ((v + 0x7efefeff) ^ ~v) & 0x81010100;
    if (hasZeroByte) // or may just have 0x80 in the high byte
    {
      hasZeroByte = ~((((v & 0x7F7F7F7F) + 0x7F7F7F7F) | v) | 0x7F7F7F7F);
    }

There is yet a faster method — use [`hasless`](#HasLessInWord)(v, 1),
which is defined below; it works in 4 operations and requires no
subsquent verification. It simplifies to

    #define haszero(v) (((v) - 0x01010101UL) & ~(v) & 0x80808080UL)

The subexpression (v - 0x01010101UL), evaluates to a high bit set in any
byte whenever the corresponding byte in v is zero or greater than 0x80.
The sub-expression ~v & 0x80808080UL evaluates to high bits set in bytes
where the byte of v doesn't have its high bit set (so the byte was less
than 0x80). Finally, by ANDing these two sub-expressions the result is
the high bits set where the bytes in v were zero, since the high bits
set due to a value greater than 0x80 in the first sub-expression are
masked off by the second.

Paul Messmer suggested the fast pretest improvement on October 2, 2004.
Juha Järvi later suggested `hasless(v, 1)` on April 6, 2005, which he
found on [Paul Hsieh's Assembly
Lab](http://www.azillionmonkeys.com/qed/asmexample.html); previously it
was written in a newsgroup post on April 27, 1987 by Alan Mycroft.

We may want to know if any byte in a word has a specific value. To do
so, we can XOR the value to test with a word that has been filled with
the byte values in which we're interested. Because XORing a value with
itself results in a zero byte and nonzero otherwise, we can pass the
result to `haszero`.

    #define hasvalue(x,n) \
    (haszero((x) ^ (~0UL/255 * (n))))

Stephen M Bennet suggested this on December 13, 2009 after reading the
entry for `haszero`.

Test if a word x contains an unsigned byte with value \&lt n.
Specifically for n=1, it can be used to find a 0-byte by examining one
long at a time, or any byte by XORing x with a mask first. Uses 4
arithmetic/logical operations when n is constant.

Requirements: x\>=0; 0\<=n\<=128

    #define hasless(x,n) (((x)-~0UL/255*(n))&~(x)&~0UL/255*128)

To count the number of bytes in x that are less than n in 7 operations,
use

    #define countless(x,n) \
    (((~0UL/255*(127+(n))-((x)&~0UL/255*127))&~(x)&~0UL/255*128)/128%255)

Juha Järvi sent this clever technique to me on April 6, 2005. The
`countless` macro was added by Sean Anderson on April 10, 2005, inspired
by Juha's `countmore`, below.

Test if a word x contains an unsigned byte with value \> n. Uses 3
arithmetic/logical operations when n is constant.

Requirements: x\>=0; 0\<=n\<=127

    #define hasmore(x,n) (((x)+~0UL/255*(127-(n))|(x))&~0UL/255*128)

To count the number of bytes in x that are more than n in 6 operations,
use:

    #define countmore(x,n) \
    (((((x)&~0UL/255*127)+~0UL/255*(127-(n))|(x))&~0UL/255*128)/128%255)

The macro `hasmore` was suggested by Juha Järvi on April 6, 2005, and he
added `countmore` on April 8, 2005.

When m \< n, this technique tests if a word x contains an unsigned byte
value, such that m \< value \< n. It uses 7 arithmetic/logical
operations when n and m are constant.

Note: Bytes that equal n can be reported by `likelyhasbetween` as false
positives, so this should be checked by character if a certain result is
needed.

Requirements: x\>=0; 0\<=m\<=127; 0\&lt=n\<=128

    #define likelyhasbetween(x,m,n) \
    ((((x)-~0UL/255*(n))&~(x)&((x)&~0UL/255*127)+~0UL/255*(127-(m)))&~0UL/255*128)

This technique would be suitable for a fast pretest. A variation that
takes one more operation (8 total for constant m and n) but provides the
exact answer is:

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

To count the number of bytes in x that are between m and n (exclusive)
in 10 operations, use:

    #define countbetween(x,m,n) (hasbetween(x,m,n)/128%255)

Juha Järvi suggested `likelyhasbetween` on April 6, 2005. From there,
Sean Anderson created `hasbetween` and `countbetween` on April 10, 2005.

Suppose we have a pattern of N bits set to 1 in an integer and we want
the next permutation of N 1 bits in a lexicographical sense. For
example, if N is 3 and the bit pattern is 00010011, the next patterns
would be 00010101, 00010110, 00011001,00011010, 00011100, 00100011, and
so forth. The following is a fast way to compute the next permutation.

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

The \_\_builtin\_ctz(v) GNU C compiler intrinsic for x86 CPUs returns
the number of trailing zeros. If you are using Microsoft compilers for
x86, the intrinsic is \_BitScanForward. These both emit a bsf
instruction, but equivalents may be available for other architectures.
If not, then consider using one of the methods for counting the
consecutive zero bits mentioned earlier.

Here is another version that tends to be slower because of its division
operator, but it does not require counting the trailing zeros.

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

Thanks to Dario Sneidermanis of Argentina, who provided this on November
28, 2009.

[A Belorussian
translation](http://webhostingrating.com/libs/bithacks-be) (provided by
[Webhostingrating](http://webhostingrating.com/)) is available.
