What is the efficient way to count set bits at a position or lower?

by

This C++ gets g++ to emit very good x86 ASM (godbolt compiler explorer). I expect it will compile efficiently on other 64bit architectures, too (if there’s a HW popcount for std::bitset::count to use, otherwise that’ll always be the slow part; e.g. sure to use g++ -march=nehalem or higher, or -mpopcnt if you don’t want to enable anything else, if you can limit your code to only running on CPUs that support that x86 instruction):

#include <bitset>

int popcount_subset(std::bitset<64> A, int pos) {
  int high_bits_to_eliminate = 63 - pos;
  A <<= (high_bits_to_eliminate & 63);  // puts A[pos] at A[63].

  return (A[63]? ~0ULL : 0) & A.count();  // most efficient way: great code with gcc and clang
  // see the godbolt link for some #ifdefs with other ways to do the check, like
    // return A[BSET_SIZE-1] ? A.count() : 0;
}

This probably isn’t optimal on 32bit architectures, so compare other alternatives if you need to make a 32bit build.

This will work for other sizes of bitset, as long as you do something about the hard-coded 63s, and change the & 63 mask for the shift count into a more general range-check. For optimal performance with strange size bitsets, make a template function with a specialization for size <= register width of the target machine. In that case, extract the bitset to an unsigned type of the appropriate width, and shift to the top of the register instead of the top of the bitset.

You’d expect this to also generate ideal code for bitset<32>, but it doesn’t quite. gcc/clang still use 64bit registers on x86-64.

For large bitsets, shifting the whole thing will be slower than just popcounting the words below the one containing pos, and using this on that word. (This is where a vectorized popcount really shines on x86 if you can assume SSSE3 but not the popcnt insn hardware support, or for 32bit targets. AVX2 256bit pshufb is the fastest way to do bulk popcounts, but without AVX2 I think 64bit popcnt is pretty close to a 128-bit pshufb implementation. See the comments for more discussion.)

If you have an array of 64-bit elements, and want to count bits below a certain position in each one separately, then you definitely should use SIMD. The shift parts of this algorithm vectorize, not just the popcnt part. Use psadbw against an all-zero register to horizontal-sum bytes in 64-bit chunks after a pshufb-based popcnt that produces counts for the bits in each byte separately. SSE/AVX doesn’t have 64-bit arithmetic right shift, but you can use a different technique to blend on the high bit of each element.

How I came up with this:

The asm instructions you want to get the compiler to output will:

  1. remove the unwanted bits from the 64bit value
  2. test the highest of the wanted bits.
  3. popcount it.
  4. return 0 or popcount, depending on the result of the test. (Branchless or branching implementations both have advantages. If the branch is predictable, a branchless implementation tends to be slower.)

The obvious way to do 1 is to generate a mask ((1<<(pos+1)) -1) and & it. A more efficient way is to left-shift by 63-pos, leaving the bits you want packed at the top of a register.

This also has the interesting side-effect of putting the bit you want to test as the top bit in the register. Testing the sign bit, rather than any other arbitrary bit, takes slightly fewer instructions. An arithmetic right shift can broadcast the sign bit to the rest of the register, allowing for more-efficient-than-usual branchless code.

Doing the popcount is a much-discussed problem, but is actually the trickier part of the puzzle. On x86, there is extremely efficient hardware support for it, but only on recent-enough hardware. On Intel CPUs, the popcnt instruction is only available on Nehalem and newer. I forget when AMD added support.

So to use it safely, you need to either do CPU dispatching with a fallback that doesn’t use popcnt. Or, make separate binaries that do/don’t depend on some CPU features.

popcount without the popcnt instruction can be done a few ways. One uses SSSE3 pshufb to implement a 4-bit LUT. This is most effective when used on a whole array, rather than a single 64b at a time, though. Scalar bithacks might be best here, and wouldn’t require SSSE3 (and so would be compatible with ancient AMD CPUs that have 64bit but not pshufb.)

The Bitbroadcast:

(A[63]? ~0ULL : 0) asks the compiler to broadcast the high bit to all other bit positions, allowing it to be used as an AND-mask to zero (or not) the popcount result. Note that even for large bitset sizes, it’s still only masking the output of popcnt, not the bitset itself, so ~0ULL is fine I used ULL to make sure wasn’t ever asking the compiler to broadcast the bit only to the low 32b of a register (with UL on Windows, for example).

This broadcast can be done with an arithmetic right shift by 63, which shifts in copies of the high bit.

clang generated this code from the original version. After some prodding from Glenn about different implementations for 4, I realized that I could lead gcc towards clang’s optimal solution by writing the source more like the ASM I want. The obvious ((int64_t)something) >> 63 to more directly request an arithmetic right shift would not be strictly portable, because signed right-shifts are implementation-defined as either arithmetic or logical. The standard doesn’t provide any portable arithmetic right-shift operator. (It’s not undefined behaviour, though.) Anyway, fortunately compilers are smart enough: gcc sees the best way once you give it enough of a hint.

This source makes great code on x86-64 and ARM64 with gcc and clang. Both simply use an arithmetic right shift on the input to popcnt (so the shift can run in parallel with the popcnt). It also compiles great on 32bit x86 with gcc, because the masking only happens to a 32bit variable (after multiple popcnt results are added). It’s the rest of the function that’s nasty on 32bit (when the bitset is larger than a register).

Original ternary-operator version with gcc

Compiled with gcc 5.3.0 -O3 -march=nehalem -mtune=haswell (older gcc, like 4.9.2, also still emit this):

; the original ternary-operator version.  See below for the optimal version we can coax gcc into emitting.
popcount_subset(std::bitset<64ul>, int):
    ; input bitset in rdi, input count in esi (SysV ABI)
    mov     ecx, esi    ; x86 variable-count shift requires the count in cl
    xor     edx, edx    ; edx=0 
    xor     eax, eax    ; gcc's workaround for popcnt's false dependency on the old value of dest, on Intel
    not     ecx         ; two's complement bithack for 63-pos (in the low bits of the register)
    sal     rdi, cl     ; rdi << ((63-pos) & 63);  same insn as shl (arithmetic == logical left shift)
    popcnt  rdx, rdi
    test    rdi, rdi    ; sets SF if the high bit is set.
    cmovs   rax, rdx    ; conditional-move on the sign flag
    ret

See How to prove that the C statement -x, ~x+1, and ~(x-1) yield the same results? for background on gcc’s use of the -x == ~x + 1 two’s complement identity. (And Which 2’s complement integer operations can be used without zeroing high bits in the inputs, if only the low part of the result is wanted? which tangentially mentions that shl masks the shift count, so we only need the low 6 bits of ecx to hold 63 - pos. Mostly linking that because I wrote it recently and anyone still reading this paragraph might find it interesting.)

Some of those instructions will go away when inlining. (e.g. gcc would generate the count in ecx in the first place.)

With Glenn’s multiply instead of ternary operator idea (enabled by USE_mul), gcc does

    shr     rdi, 63
    imul    eax, edi

at the end instead of xor / test / cmovs.

Haswell perf analysis, using microarch data from Agner Fog (Multiply version):

  • mov r,r: 1 fused-domain uop, 0 latency, no execution unit
  • xor-zeroing: 1 fused-domain uop, no execution unit
  • not: 1 uop for p0/p1/p5/p6, 1c latency, 1 per 0.25c throughput
  • shl (aka sal) with count in cl: 3 uops for p0/p6: 2c latency, 1 per 2c throughput. (Agner Fog’s data indicates that IvyBridge only takes 2 uops for this, strangely.)
  • popcnt: 1 uop for p1, 3c latency, 1 per 1c throughput
  • shr r,imm: 1 uop for p0/p6, 1c latency. 1 per 0.5c throughput.
  • imul r,r: 1uop for p1, 3c latency.
  • not counting the ret

Totals:

  • 9 fused-domain uops, can issue in 2.25 cycles (in theory; uop cache-line effects usually bottleneck the frontend slightly).
  • 4 uops (shifts) for p0/p6. 2 uops for p1. 1 any-ALU-port uop. Can execute at one per 2c (saturating the shift ports), so the frontend is the worst bottleneck.

Latency: Critical path from when the bitset is ready to when the result is: shl(2) -> popcnt(3) -> imul(3). Total 8 cycles. Or 9c from when pos is ready, because the not is an extra 1c latency for it.

The optimal bitbroadcast version replaces shr with sar (same perf), and imul with and (1c latency instead of 3c, runs on any port). So the only perf change is reducing the critical path latency to 6 cycles. Throughput is still bottlenecked on the frontend. and being able to run on any port doesn’t make a difference, unless you’re mixing this with code that bottlenecks on port1 (instead of looking at the throughput for running just this code in a tight loop).

cmov (ternary operator) version: 11 fused-domain uops (frontend: one per 2.75c). execution units: still bottlenecked on the shift ports (p0/p6) at one per 2c. Latency: 7c from bitset to result, 8c from pos to result. (cmov is 2c latency, 2 uops for any of p0/p1/p5/p6.)

Clang has some different tricks up its sleeve: Instead of test/cmovs, it generates a mask of either all-ones or all-zeros by using an arithmetic right-shift to broadcast the sign bit to all positions of a register. I love it: Using and instead of cmov is more efficient on Intel. It still has the data-dependency and does the work for both sides of the branch (which is the main downside to cmov in general), though. Update: with the right source code, gcc will use this method, too.

clang 3.7 -O3 -Wall -march=nehalem -mtune=haswell 

popcount_subset(std::bitset<64ul>, int):
    mov     ecx, 63
    sub     ecx, esi      ; larger code size, but faster on CPUs without mov-elimination
    shl     rdi, cl       ; rdi << ((63-pos) & 63)
    popcnt  rax, rdi      ; doesn't start a fresh dep chain before this, like gcc does
    sar     rdi, 63       ; broadcast the sign bit
    and     eax, edi      ; eax = 0 or its previous value
    ret

sar / and replaces xor / test / cmov, and cmov is a 2-uop instruction on Intel CPUs, so that’s really nice. (For the ternary-operator version).

Clang still does the sar / and trick instead of an actual imul when using the multiply source version, or the “bitbroadcast” source version. So those help gcc without hurting clang. (sar/and is definitely better than shr/imul: 2c less latency on the critical path.) The pow_of_two_sub version does hurt clang (see the first godbolt link: omitted from this answer to avoid clutter with ideas that didn’t pan out).

The mov ecx, 63 / sub ecx, esi is actually faster on CPUs without mov-elimination for reg,reg moves (zero latency and no execution port, handled by register renaming). This includes Intel pre-IvyBridge, but not more recent Intel and AMD CPUs.

Clang’s mov imm / sub method puts only one cycle of latency for pos onto the critical path (beyond the bitset->result latency), instead of two for a mov ecx, esi / not ecx on CPUs where mov r,r has 1c latency.

With BMI2 (Haswell and later), an optimal ASM version can save a mov to ecx. Everything else works the same, because shlx masks its shift-count input register down to the operand-size, just like shl.

x86 shift instructions have crazy CISC semantics where if the shift count is zero, the flags aren’t affected. So variable-count shift instructions have a (potential) dependency on the old value of the flags. “Normal” x86 shl r, cl decodes to 3 uops on Haswell, but BMI2 shlx r, r, r is only 1. So it’s too bad that gcc still emits sal with -march=haswell, instead of using shlx (which it does use in some other cases).

// hand-tuned BMI2 version using the NOT trick and the bitbroadcast
popcount_subset(std::bitset<64ul>, int):
    not     esi           ; The low 6 bits hold 63-pos.  gcc's two-s complement trick
    xor     eax, eax      ; break false dependency on Intel.  maybe not needed when inlined.
    shlx    rdi, rdi, rsi ; rdi << ((63-pos) & 63)
    popcnt  rax, rdi
    sar     rdi, 63       ; broadcast the sign bit: rdi=0 or -1
    and     eax, edi      ; eax = 0 or its previous value
    ret

Perf analysis for Intel Haswell: 6 fused-domain uops (frontend: one per 1.5c). Execution units: 2 p0/p6 shift uops. 1 p1 uop. 2 any-port uops: (one per 1.25c from total execution port limits). Critical path latency: shlx(1) -> popcnt(3) -> and(1) = 5c bitset->result. (or 6c from pos->result).

Note that when inlining, a human (or smart compiler) could avoid the need for the xor eax, eax. It’s only there because of popcnt‘s false dependency on the output register (on Intel), and we need the output in eax (which the caller may have used recently for a long dep chain). With -mtune=bdver2 or something, gcc won’t zero the register it’s going to use for popcnt output.

When inlining, we could use an output register that already has to be ready at least as early as popcnt‘s source reg to avoid the problem. Compilers will do an in-place popcnt rdi,rdi when the source isn’t needed later, but that’s not the case here. Instead, we can pick another register that already has to be ready before the source. popcnt‘s input depends on 63-pos, and we can clobber it, so popcnt rsi,rdi‘s dependency on rsi can’t delay it. Or if we had 63 in a register, we could popcnt rsi,rdi / sarx rax, rsi, reg_63 / and eax, esi. Or BMI2 3-operand shift instructions would also let us not clobber inputs in case they’re needed afterwards.

This is so light-weight that loop overhead and setting up the input operands / storing the results are going to be major factors. (And the 63-pos can optimize away with a compile-time constant, or into wherever a variable count comes from.)

The Intel compiler amusingly shoots itself in the foot and doesn’t take advantage of the fact that A[63] is the sign bit. shl / bt rdi, 63 / jc. It even sets up the branches in a really dumb way. It could zero eax, and then jump over popcnt or not based on the sign flag set by shl.

An optimal branching implementation, starting from ICC13 output from -O3 -march=corei7 on godbolt:

   // hand-tuned, not compiler output
        mov       ecx, esi    ; ICC uses neg/add/mov :/
        not       ecx
        xor       eax, eax    ; breaks the false dep, or is the return value in the taken-branch case
        shl       rdi, cl
        jns    .bit_not_set
        popcnt    rax, rdi
.bit_not_set:
        ret

That’s pretty much optimal: The A[pos] == true case has one not-taken branch. It doesn’t save very much over the branchless method, though.

If the A[pos] == false case is more common: jump over a ret instruction, to a popcnt / ret. (Or after inlining: jump to a block at the end that does the popcnt and jumps back).

More Related Contents:

  1. How to add two numbers without using ++ or + or another arithmetic operator
  2. Generating m distinct random numbers in the range [0..n-1]
  3. What’s the fastest way to pack 32 0/1 values into the bits of a single 32-bit variable?
  4. Is it possible to map string to int faster than using hashmap?
  5. Fastest way to get the integer part of sqrt(n)?
  6. How can I test whether a number is a power of 2?
  7. Invalid solution for code challenge with operator restrictions
  8. Why does C++ code for testing the Collatz conjecture run faster than hand-written assembly?
  9. Position of least significant bit that is set
  10. Is there any advantage of using map over unordered_map in case of trivial keys?
  11. Static linking vs dynamic linking
  12. Finding all the subsets of a set
  13. Efficient implementation of log2(__m256d) in AVX2
  14. Is inline assembly language slower than native C++ code?
  15. read line by line in the most efficient way *platform specific*
  16. Efficient Cartesian Product algorithm
  17. 5 years later, is there something better than the “Fastest Possible C++ Delegates”?
  18. Printing prime numbers from 1 through 100
  19. Circular lock-free buffer
  20. Efficient bitwise operations for counting bits or find the right|left most ones
  21. c++11 regex slower than python
  22. Efficiently getting all divisors of a given number
  23. What’s the quickest way to compute log2 of an integer in C#?
  24. Using bts assembly instruction with gcc compiler
  25. What is a good random number generator for a game?
  26. Branchless code that maps zero, negative, and positive to 0, 1, 2
  27. Is it more efficient to copy a vector by reserving and copying, or by creating and swapping? [duplicate]
  28. Why do I see 400x outlier timings when calling clock_gettime repeatedly?
  29. What is the modern, correct way to do type punning in C++?
  30. Why is this C++ wrapper class not being inlined away?

Leave a Comment