title: What is the fastest way to return the positions of all set bits in a 64-bit integer?
body: <p>I need a fast way to get the position of all one bits in a 64-bit integer.  For example, given <code>x = 123703</code>, I'd like to fill an array <code>idx[] = {0, 1, 2, 4, 5, 8, 9, 13, 14, 15, 16}</code>. We can assume we know the number of bits a priori.  This will be called 10<sup>12</sup> - 10<sup>15</sup> times, so speed is of the essence.  The fastest answer I've come up with so far is the following monstrosity, which uses each byte of the 64-bit integer as an index into tables that give the number of bits set in that byte and the positions of the ones:</p>
  <pre><code>int64_t x;            // this is the input
  unsigned char idx[K]; // this is the array of K bits that are set
  unsigned char *dst=idx, *src;
  unsigned char zero, one, two, three, four, five;  // these hold the 0th-5th bytes
  zero  =  x &amp; 0x0000000000FFUL;
  one   = (x &amp; 0x00000000FF00UL) &gt;&gt; 8;
  two   = (x &amp; 0x000000FF0000UL) &gt;&gt; 16;
  three = (x &amp; 0x0000FF000000UL) &gt;&gt; 24;
  four  = (x &amp; 0x00FF00000000UL) &gt;&gt; 32;
  five  = (x &amp; 0xFF0000000000UL) &gt;&gt; 40;
  src=tab0+tabofs[zero ]; COPY(dst, src, n[zero ]);
  src=tab1+tabofs[one  ]; COPY(dst, src, n[one  ]);
  src=tab2+tabofs[two  ]; COPY(dst, src, n[two  ]);
  src=tab3+tabofs[three]; COPY(dst, src, n[three]);
  src=tab4+tabofs[four ]; COPY(dst, src, n[four ]);
  src=tab5+tabofs[five ]; COPY(dst, src, n[five ]);
  </code></pre>
  <p>where <code>COPY</code> is a switch statement to copy up to 8 bytes, <code>n</code> is array of the number of bits set in a byte and <code>tabofs</code> gives the offset into <code>tabX</code>, which holds the positions of the set bits in the X-th byte.  <del>This is about 3x faster than unrolled loop-based methods with <code>__builtin_ctz()</code> on my Xeon E5-2609.</del> (See below.)  I am currently iterating <code>x</code> in lexicographical order for a given number of bits set.</p>
  <p>Is there a better way?</p>
  <p><strong>EDIT</strong>: Added an example (that I have subsequently fixed).  Full code is available here: <a href="http://pastebin.com/79X8XL2P" rel="nofollow noreferrer">http://pastebin.com/79X8XL2P</a> .  Note: GCC with -O2 seems to optimize it away, but Intel's compiler (which I used to compose it) doesn't...</p>
  <p>Also, let me give some additional background to address some of the comments below.  The  goal is to perform a statistical test on every possible subset of K variables out of a universe of N possible explanatory variables; the specific target right now is N=41, but I can see some projects needing N up to 45-50.  The test basically involves factorizing the corresponding data submatrix.  In pseudocode, something like this:</p>
  <pre><code>double doTest(double *data, int64_t model) {
    int nidx, idx[];
    double submatrix[][];
    nidx = getIndices(model, idx);  // get the locations of ones in model
    // copy data into submatrix
    for(int i=0; i&lt;nidx; i++) {
      for(int j=0; j&lt;nidx; j++) {
        submatrix[i][j] = data[idx[i]][idx[j]];
      }
    }
    factorize(submatrix, nidx);
    return the_answer;
  }
  </code></pre>
  <p>I coded up a version of this for an Intel Phi board that should complete the N=41 case in about 15 days, of which ~5-10% of the time is spent in a naive <code>getIndices()</code> so right off the bat a faster version could save a day or more.  I'm working on an implementation for NVidia Kepler too, but unfortunately the problem I have (ludicrous numbers of small matrix operations) is not ideally suited to the hardware (ludicrously large matrix operations).  That said, <a href="http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6267819" rel="nofollow noreferrer">this paper</a> presents a solution that seems to achieve hundreds of GFLOPS/s on matrices of my size by aggressively unrolling loops and performing the entire factorization in registers, with the caveat that the dimensions of the matrix be defined at compile-time.  (This loop unrolling should help reduce overhead and improve vectorization in the Phi version too, so <code>getIndices()</code> will become more important!)  So now I'm thinking my kernel should look more like:</p>
  <pre><code>double *data;  // move data to GPU/Phi once into shared memory
  template&lt;unsigned int K&gt; double doTestUnrolled(int *idx) {
    double submatrix[K][K];
    // copy data into submatrix
    #pragma unroll
    for(int i=0; i&lt;K; i++) {
      #pragma unroll
      for(int j=0; j&lt;K; j++) {
        submatrix[i][j] = data[idx[i]][idx[j]];
      }
    }
    factorizeUnrolled&lt;K&gt;(submatrix);
    return the_answer;
  }
  </code></pre>
  <p>The Phi version solves each model in a `cilk_for' loop from model=0 to 2<sup>N</sup> (or, rather, a subset for testing), but now in order to batch work for the GPU and amortize the kernel launch overhead I have to iterate model numbers in lexicographical order for each of K=1 to 41 bits set (as doynax noted).</p>
  <p><strong>EDIT 2:</strong>  Now that vacation is over, here are some results on my Xeon E5-2602 using icc version 15.  The code that I used to benchmark is here: <a href="http://pastebin.com/XvrGQUat" rel="nofollow noreferrer">http://pastebin.com/XvrGQUat</a>.  I perform the bit extraction on integers that have exactly K bits set, so there is some overhead for the lexicographic iteration measured in the &quot;Base&quot; column in the table below.  These are performed 2<sup>30</sup> times with N=48 (repeating as necessary).</p>
  <p>&quot;CTZ&quot; is a loop that uses the the gcc intrinsic <code>__builtin_ctzll</code> to get the lowest order bit set:</p>
  <pre><code>for(int i=0; i&lt;K; i++) {
      idx[i] = __builtin_ctzll(tmp);
      lb = tmp &amp; -tmp;    // get lowest bit
      tmp ^= lb;      // remove lowest bit from tmp
  } 
  </code></pre>
  <p>Mark is Mark's branchless for loop:</p>
  <pre><code>for(int i=0; i&lt;K; i++) {
      *dst = i;
      dst += x &amp; 1;
      x &gt;&gt;= 1;
  } 
  </code></pre>
  <p>Tab1 is my original table-based code with the following copy macro:</p>
  <pre><code>#define COPY(d, s, n) \
  switch(n) { \
  case 8: *(d++) = *(s++); \
  case 7: *(d++) = *(s++); \
  case 6: *(d++) = *(s++); \
  case 5: *(d++) = *(s++); \
  case 4: *(d++) = *(s++); \
  case 3: *(d++) = *(s++); \
  case 2: *(d++) = *(s++); \
  case 1: *(d++) = *(s++); \
  case 0: break;        \
  }
  </code></pre>
  <p>Tab2 is the same code as Tab1, but the copy macro just moves 8 bytes  as a single copy (taking ideas from doynax and Lưu Vĩnh Phúc... but note this does <em>not</em> ensure alignment):</p>
  <pre><code>#define COPY2(d, s, n) { *((uint64_t *)d) = *((uint64_t *)s); d+=n; }
  </code></pre>
  <p>Here are the results.  I guess my initial claim that Tab1 is 3x faster than CTZ only holds for large K (where I was testing).  Mark's loop is faster than my original code, but getting rid of the branch in the <code>COPY2</code> macro takes the cake for K &gt; 8.</p>
  <pre><code> K    Base    CTZ   Mark   Tab1   Tab2
  001  4.97s  6.42s  6.66s 18.23s 12.77s
  002  4.95s  8.49s  7.28s 19.50s 12.33s
  004  4.95s  9.83s  8.68s 19.74s 11.92s
  006  4.95s 16.86s  9.53s 20.48s 11.66s
  008  4.95s 19.21s 13.87s 20.77s 11.92s
  010  4.95s 21.53s 13.09s 21.02s 11.28s
  015  4.95s 32.64s 17.75s 23.30s 10.98s
  020  4.99s 42.00s 21.75s 27.15s 10.96s
  030  5.00s 100.64s 35.48s 35.84s 11.07s
  040  5.01s 131.96s 44.55s 44.51s 11.58s
  </code></pre>