英特尔位操作指令集 2 (BMI2) 中的并行存款指令 (
PDEPU64 _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;
}
pdep这个算法是 O(n)，其中 n 是
maskmask是否可能有更有效的最坏情况算法？
是否可以制作一个更快的版本，假设
val1<<rr问题的第二部分，关于 1 位存款的特殊情况，需要两个步骤。第一步，我们需要确定
rvalvalffsrirmaskrvali查找
rmask以下
Cmaskval#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;
}
作为补充说明，由于这个问题再次出现在我身上，我发现Sebastiano Vigna 的方法在实践中更快地找到第n组位。它不包含分支或条件移动。
请注意，
leq_bytesgt_zero_bytesconst 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);
}
确实有。
我刚刚完成了以下例程的编码。有多个优化点，但它在 O(lg^2 N) 渐近时间内完成了任务。
#include <cinttypes>
uint64_t pdep64r (uint64_t src, uint64_t mask) {
static const int W = 8*sizeof(uint64_t); // bits in word
static const int LGW = 6; // log2(64)
static const uint64_t UMASK[]= {0x0001, 0x0000000100000001, 0x0001000100010001, 0x0101010101010101,
0x1111111111111111, 0x5555555555555555};
static const uint64_t XMASK[]= {0xFFFFFFFFFFFFFFFF, 0xFFFFFFFF, 0xFFFF, 0xFF, 0x0F, 0x03};
static const uint64_t MXOR[] = {0x00000000FFFFFFFF, 0x0000FFFF0000FFFF,
0x00FF00FF00FF00FF, 0x0F0F0F0F0F0F0F0F,
0x3333333333333333, 0x5555555555555555};
// recursive popcount every 2^k subblock of the mask and store the counts
uint64_t pcmask[LGW+1]; {
pcmask[LGW] = mask;
for (int k=LGW; k>0; --k) {
pcmask[k-1] = (pcmask[k] & MXOR[k-1]) +
((pcmask[k] & ~MXOR[k-1]) >> (1<<(LGW-k) ) );
}
}
uint64_t dst = src;
// solve by recursively splitting the problems into 2^k x 2^(N-k) blocks
for (int k=0; k<LGW; ++k) {
uint64_t pclo = pcmask[k+1] & MXOR[k]; // popcount of half
const uint64_t UNIT = UMASK[k];
const uint64_t EXT = XMASK[k];
uint64_t m = UNIT;
for (int n=0; n<LGW-k; ++n) {
uint64_t d = ( (pclo>>n) & UNIT)*EXT; // select blocks with d0=1
m = ( (m<<(1<<n)) & d ) | (m & ~d); // expand the mask m by a power of 2 if that digit is set
}
m -= UNIT;
uint64_t l = dst & m; // keep in lower half
uint64_t h = dst & ~m; // must be shifted to upper half
for (int n=0; n<LGW-k; ++n) {
uint64_t d = ( (pclo>>n) & UNIT)*EXT; // select blocks with d0=1
h = ( (h>>(1<<n)) & d ) | (h & ~d); // shift <<(2^d) only for blocks with d set
}
h <<= (1<<(LGW-k-1));
h &= ~MXOR[k]; // remove spillovers to next block
dst = l | h;
}
// as we only introduced zeros by shifting left there may be leftover bits that need cancelling
return dst & mask;
}