Reversing a 64-bit word

Here’s a cute way of reversing a 64-bit word (there is a similar, slightly faster, but slightly more obscure method in TAOCP 4a, so if you really need to save the cycles, use that one).

We use ternary swaps: for example, to bit-reverse a 3-bit sequence, we can just swap the high and low bits. To bit-reverse a 9-bit segment, break into [3,3,3] and swap the top and bottom 3-bit sections, and then bit-reverse each 3-bit section separately.

More generally, we can reverse any bit sequence by breaking it up into [M,N,M], swapping top and bottom sections (by shifting up and down by M+N), and recursing.

So, to reverse a 63 bit number, since 63 = 3*3*7, we can nest two 3-swaps and a 7-swap. We can do the actual swapping of sections with another classy Knuth function, nicely expressible as a C++ template:

template <typename T, T m, int k>
static inline T swapbits(T p) {
  T q = ((p>>k)^p)&m;
  return p^q^(q<<k);
}

m is a mask, k a shift amount, the function swaps the bits at positions given by the set bits of m and m<<k (clearly these two should be disjoint, ie. m & m<<k == 0)

All we need to do now is calculate the masks and the shift amounts.

First 3-swap; the binary mask is ...001001001 and we need to shift the mask up by 2. We can see that m1 + m1<<1 + m1<<2 = 2^63-1, so we can calculate the mask with a (compile-time) division.

Second 3-swap: the binary mask is ...000000111000000111 and we need to shift the mask up by 6. Once again, we can easily compute the correct mask.

The 7-swap (ie. reverse 7 sections of 9 bits), we do in two stages, do a 3-swap for the top and bottom sections, so the mask is 111111111 + 111111111<<36, finally we swap the top and bottom 27 bits, so the mask is just 2^27-1, and the shift is 36.

This reverses the bottom 63 bits, a simple rotate by one then puts everything into the right place.

uint64_t bitreverse (uint64_t n)
{
  static const uint64_t m1 = ((uint64_t(1)<<63)-1)/(1+(1<<1)+(1<<2));
  static const uint64_t m2 = ((uint64_t(1)<<63)-1)/(1+(1<<3)+(1<<6));
  static const uint64_t m3 = ((uint64_t(1)<<9)-1)+(((uint64_t(1)<<9)-1)<<36);
  static const uint64_t m4 = (uint64_t(1)<<27)-1;
  n = swapbits<uint64_t, m1, 2>(n);
  n = swapbits<uint64_t, m2, 6>(n);
  n = swapbits<uint64_t, m3, 18>(n);
  n = swapbits<uint64_t, m4, 36>(n);
  n = (n >> 63) | (n << 1);
  return n;
}

Here is what gcc makes of that:

_Z10bitreversey:
	movq	%rdi, %rdx
	movabsq	$1317624576693539401, %rax
	movabsq	$126347562148695559, %rcx
	shrq	$2, %rdx
	xorq	%rdi, %rdx
	andq	%rax, %rdx
	movq	%rdx, %rax
	salq	$2, %rdx
	xorq	%rdi, %rax
	xorq	%rdx, %rax
	movq	%rax, %rdx
	shrq	$6, %rdx
	xorq	%rax, %rdx
	andq	%rcx, %rdx
	movabsq	$35115652612607, %rcx
	xorq	%rdx, %rax
	salq	$6, %rdx
	xorq	%rdx, %rax
	movq	%rax, %rdx
	shrq	$18, %rdx
	xorq	%rax, %rdx
	andq	%rcx, %rdx
	xorq	%rdx, %rax
	salq	$18, %rdx
	xorq	%rdx, %rax
	movq	%rax, %rdx
	shrq	$36, %rdx
	xorq	%rax, %rdx
	andl	$134217727, %edx
	xorq	%rdx, %rax
	salq	$36, %rdx
	xorq	%rdx, %rax
	rorq	$63, %rax
	ret

For comparison, here is Knuth’s 64-bit reverse (I’ve just hard-coded the constants this time). It’s based on a 32-bit reverse that breaks one 17-bit segment into [7,3,7] and the remaining 15-bit segment into [3,7,3] – we can do both swaps with the same shift of 10. First step is to swap adjacent bits which can be done slightly faster than a general swap. Very cunning:

uint64_t kbitreverse (uint64_t n)
{
  static const uint64_t m0 = 0x5555555555555555LLU;
  static const uint64_t m1 = 0x0300c0303030c303LLU;
  static const uint64_t m2 = 0x00c0300c03f0003fLLU;
  static const uint64_t m3 = 0x00000ffc00003fffLLU;
  n = ((n>>1)&m0) | (n&m0)<<1;
  n = swapbits<uint64_t, m1, 4>(n);
  n = swapbits<uint64_t, m2, 8>(n);
  n = swapbits<uint64_t, m3, 20>(n);
  n = (n >> 34) | (n << 30);
  return n;
}

and the corresponding compiler output:

_Z11kbitreversey:
	movq	%rdi, %rdx
	movabsq	$6148914691236517205, %rax
	movabsq	$216384095313249027, %rcx
	shrq	%rdx
	andq	%rax, %rdx
	andq	%rdi, %rax
	addq	%rax, %rax
	orq	%rdx, %rax
	movq	%rax, %rdx
	shrq	$4, %rdx
	xorq	%rax, %rdx
	andq	%rcx, %rdx
	movabsq	$54096023692247103, %rcx
	xorq	%rdx, %rax
	salq	$4, %rdx
	xorq	%rdx, %rax
	movq	%rax, %rdx
	shrq	$8, %rdx
	xorq	%rax, %rdx
	andq	%rcx, %rdx
	movabsq	$17575006191615, %rcx
	xorq	%rdx, %rax
	salq	$8, %rdx
	xorq	%rdx, %rax
	movq	%rax, %rdx
	shrq	$20, %rdx
	xorq	%rax, %rdx
	andq	%rcx, %rdx
	xorq	%rdx, %rax
	salq	$20, %rdx
	xorq	%rdx, %rax
	rorq	$34, %rax
	ret

Knuth wins here by 1 instruction! Oh well, maybe I’ll have better luck next time…

Templated bit twiddling

Sometimes it’s nice to combine C++ template metaprogramming with some bit twiddling, here’s a classic parallelized bit counting algorithm with all the magic numbers constructed as part of template instantiation:

#include <stdio.h>
#include <stdlib.h>

template <typename T, int k, T m>
struct A
{
  static T f(T n)
  {
    static const int k1 = k/2;
    n = A<T, k1, m^(m<<k1)>::f(n);
    return ((n>>k)&m) + (n&m);
  }
};

template <typename T, T m>
struct A<T, 0, m>
{
  static T f(T n) { return n; }
};

template<typename T>
static inline int bitcount (T n)
{
  static const int k = 4 * sizeof(n);
  return A<T, k, (T(1)<<k)-1>::f(n);
}

int main(int argc, char *argv[])
{
  unsigned long n = strtoul(argv[1], NULL, 10);
  printf("bitcount(%lu) = %d\n", n, bitcount(n));
}

We could extract the magic number construction out as a separate template, and there are some extra optimizations that this code doesn’t do, for example, the level 1 code can be simplified to:

i = i - ((i >> 1) & 0x55555555);

rather than what we get from the above:

i = ((i >> 1)&0x55555555) + (i&0x55555555);

Adding an appropriate partial template instantiation for A<T,1,m> is left as an exercise.

In fact, there are lots of shortcuts that can be made here, we aren’t really being competitive with Bit Twiddling Hacks:

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

which would need 3 template specializations and we are starting to get a bit complicated…

In the same style then, here’s a bit reverse function, with a separate magic number calculating template:

#include <stdio.h>
#include <stdlib.h>

template <typename T, int k>
struct Magic {
  static const T m = ~T(0)/(T(1)+(T(1)<<k));
};

template <typename T, int k>
struct RevAux {
  static inline T f(T n) {
    T n1 = RevAux<T,k/2>::f(n);
    return ((n1>>k) & Magic<T,k>::m) | ((n1 & Magic<T,k>::m)<<k);
  }
};

template <typename T>
struct RevAux<T,0> {
  static inline T f(T n) { return n; }
};

template <typename T>
T rev(T n) { return RevAux<T,4*sizeof(T)>::f(n); }

int main(int argc, char *argv[])
{
  unsigned long n = strtoul(argv[1],NULL,16);
  printf("%lx %lx\n", n, rev(n));
}

Now, how do we check that the template parameter T is unsigned (otherwise the right shifts are probably going to go wrong)…

A simple optimization is to replace the top level function with:

template <typename T>
static inline T rev(T n) { 
  static const int k = 4*sizeof(T);
  n = RevAux<T,k/2>::f(n); 
  return (n>>k)|(n<<k);
}

but GCC seems to have worked that out for itself (compiling gcc 4.4.3 , -O3 -S):

.globl _Z7reversej
        pushl   %ebp
        movl    %esp, %ebp
        movl    8(%ebp), %eax
        popl    %ebp
        movl    %eax, %edx
        andl    $1431655765, %edx
        shrl    %eax
        addl    %edx, %edx
        andl    $1431655765, %eax
        orl     %eax, %edx
        movl    %edx, %eax
        andl    $858993459, %eax
        shrl    $2, %edx
        andl    $858993459, %edx
        sall    $2, %eax
        orl     %edx, %eax
        movl    %eax, %edx
        andl    $252645135, %edx
        shrl    $4, %eax
        andl    $252645135, %eax
        sall    $4, %edx
        orl     %eax, %edx
        movl    %edx, %eax
        andl    $16711935, %eax
        shrl    $8, %edx
        sall    $8, %eax
        andl    $16711935, %edx
        orl     %edx, %eax
        rorl    $16, %eax
        ret

This sort of thing isn’t going to work very well if the compiler isn’t up to scratch, but GCC seems to have done it right here.

There are slightly faster ways to reverse bits using ternary swap (swap top and bottom thirds, recurse), but they need to be more ad hoc unless you have a ternary computer with word length a power of 3, so less amenable to this sort of generic treatment (or at least, it becomes a lot harder).