Portable efficient alternative to PDEP without using BMI2?

Question

The documentation for the parallel deposit instruction (PDEP) in Intel's Bit Manipulation Instruction Set 2 (BMI2) describes the following serial implementation for the instruction (C-like pseudocode):

U64 _pdep_u64(U64 val, U64 mask) {
  U64 res = 0;
  for (U64 bb = 1; mask; bb += bb) {
    if (val & bb)
      res |= mask & -mask;
    mask &= mask - 1;
  }
  return res;
}

See also Intel's pdep insn ref manual entry.

This algorithm is O(n), where n is the number of set bits in mask, which obviously has a worst case of O(k) where k is the total number of bits in mask.

Is a more efficient worst case algorithm possible?

Is it possible to make a faster version that assumes that val has at most one bit set, ie either equals 0 or equals 1<<r for some value of r from 0 to 63?

Answer 1

The second part of the question, about the special case of a 1-bit deposit, requires two steps. In the first step, we need to determine the bit index r of the single 1-bit in val, with a suitable response in case val is zero. This can easily be accomplished via the POSIX function ffs, or if r is known by other means, as alluded to by the asker in comments. In the second step we need to identify bit index i of the r-th 1-bit in mask, if it exists. We can then deposit the r-th bit of val at bit i.

One way of finding the index of the r-th 1-bit in mask is to tally the 1-bits using a classical population count algorithm based on binary partitioning, and record all of the intermediate group-wise bit counts. We then perform a binary search on the recorded bit-count data to identify the position of the desired bit.

The following C-code demonstrates this using 64-bit data. Whether this is actually faster than the iterative method will very much depend on typical values of mask and val.

#include <stdint.h>

/* Find the index of the n-th 1-bit in mask, n >= 0
   The index of the least significant bit is 0 
   Return -1 if there is no such bit
*/
int find_nth_set_bit (uint64_t mask, int n)
{
    int t, i = n, r = 0;
    const uint64_t m1 = 0x5555555555555555ULL; // even bits
    const uint64_t m2 = 0x3333333333333333ULL; // even 2-bit groups
    const uint64_t m4 = 0x0f0f0f0f0f0f0f0fULL; // even nibbles
    const uint64_t m8 = 0x00ff00ff00ff00ffULL; // even bytes
    uint64_t c1 = mask;
    uint64_t c2 = c1 - ((c1 >> 1) & m1);
    uint64_t c4 = ((c2 >> 2) & m2) + (c2 & m2);
    uint64_t c8 = ((c4 >> 4) + c4) & m4;
    uint64_t c16 = ((c8 >> 8) + c8) & m8;
    uint64_t c32 = (c16 >> 16) + c16;
    int c64 = (int)(((c32 >> 32) + c32) & 0x7f);
    t = (c32    ) & 0x3f; if (i >= t) { r += 32; i -= t; }
    t = (c16>> r) & 0x1f; if (i >= t) { r += 16; i -= t; }
    t = (c8 >> r) & 0x0f; if (i >= t) { r +=  8; i -= t; }
    t = (c4 >> r) & 0x07; if (i >= t) { r +=  4; i -= t; }
    t = (c2 >> r) & 0x03; if (i >= t) { r +=  2; i -= t; }
    t = (c1 >> r) & 0x01; if (i >= t) { r +=  1;         }
    if (n >= c64) r = -1;
    return r; 
}

/* val is either zero or has a single 1-bit.
   Return -1 if val is zero, otherwise the index of the 1-bit
   The index of the least significant bit is 0
*/
int find_bit_index (uint64_t val)
{
    return ffsll (val) - 1;
}

uint64_t deposit_single_bit (uint64_t val, uint64_t mask)
{
    uint64_t res = (uint64_t)0;
    int r = find_bit_index (val);
    if (r >= 0) {
        int i = find_nth_set_bit (mask, r);
        if (i >= 0) res = (uint64_t)1 << i;
    } 
    return res;
}

Answer 2

Just as an additional note, since this came up for me again, I found Sebastiano Vigna's method for finding the nth set bit to be faster in practice. It contains no branches or conditional moves.

Note that leq_bytes and gt_zero_bytes might be implementable using SSE instructions, but this version the advantage of being completely portable.

const uint64_t sMSBs8 = 0x8080808080808080ull;
const uint64_t sLSBs8 = 0x0101010101010101ull;
    
inline uint64_t
leq_bytes(uint64_t pX, uint64_t pY)
{
    return ((((pY | sMSBs8) - (pX & ~sMSBs8)) ^ pX ^ pY) & sMSBs8) >> 7;
}
    
    
inline uint64_t
gt_zero_bytes(uint64_t pX)
{
    return ((pX | ((pX | sMSBs8) - sLSBs8)) & sMSBs8) >> 7;
}
    
    
inline uint64_t find_nth_set_bit(uint64_t pWord, uint64_t pR)
{
    const uint64_t sOnesStep4  = 0x1111111111111111ull;
    const uint64_t sIncrStep8  = 0x8040201008040201ull;

    uint64_t byte_sums = pWord - ((pWord & 0xA*sOnesStep4) >> 1);
    byte_sums = (byte_sums & 3*sOnesStep4) + ((byte_sums >> 2) & 3*sOnesStep4);
    byte_sums = (byte_sums + (byte_sums >> 4)) & 0xF*sLSBs8;
    byte_sums *= sLSBs8;
    
    const uint64_t k_step_8 = pR * sLSBs8;
    const uint64_t place
        = (leq_bytes( byte_sums, k_step_8 ) * sLSBs8 >> 53) & ~0x7;
    const int byte_rank = pR - (((byte_sums << 8) >> place) & 0xFF);
    const uint64_t spread_bits = (pWord >> place & 0xFF) * sLSBs8 & sIncrStep8;
    const uint64_t bit_sums = gt_zero_bytes(spread_bits) * sLSBs8;
    const uint64_t byte_rank_step_8 = byte_rank * sLSBs8;
    return place + (leq_bytes( bit_sums, byte_rank_step_8 ) * sLSBs8 >> 56);
}