First actually understand what the problem is. Population count (also called Hamming weights) is the count of 1 bits in a string. By count of 1 bits you might have guessed it right that we are talking about a binary string. For instance, 11000101 -> the population count is 4.
There is a compiler extension for this __builtin_popcount(int n) for this.
So the problem which we are solving is that we are given an array of binary strings and we have
to calculate the sum of all 1's in all the strings.
int population_count(int *arr, int n){
int sum = 0;
for(int i = 0; i < n; i++){
sum += __builtin_popcount(arr[i]);
}
return sum;
}
This assumes a x86-64 architecture, which has a seaprate instruction dedicated for this.
Now since we have simd, we can try solving this in a different way. I came across this while
reading the Algorithmica Book[1].
This approach was introduced by Wojciech Muła as mentioned in the book.
For simplification, we take a array of 8 integers, since they fit on a single 256 bit wide register(ymm registers).
PSHUFB stands for packed shuffle bytes. The intrinsic that corresponds to it is _mm256_shuffle_epi8(arr, mask). It returns the values, shuffled as the indexes that are passed as mask.
For example: arr is [1, 2, 3, 4, 5] and the mask is [3, 2, 1, 4, 0], then the result
will be [4, 3, 2, 5, 1].
The thing with PSHUFB is that even if you use it on a 256 wide register, it actually operates on 2 lanes of 128 bits. Also the indexes that it takes as input have to be 4 bits.
Since the shuffle intrinsic takes 4 bit indexes, we can make the table for 0 -> 15, each index corresponding to its population count. Since it operates on 2 separate lanes, the lookup table will also be used over 2 lanes. For this reason we make the table twice in the lookup, the first one works on the lower lane and the second one on the higher lane.
const __m256i lookup = _mm256_setr_epi8(
/* 0 */ 0, /* 1 */ 1, /* 2 */ 1, /* 3 */ 2,
/* 4 */ 1, /* 5 */ 2, /* 6 */ 2, /* 7 */ 3,
/* 8 */ 1, /* 9 */ 2, /* A */ 2, /* B */ 3,
/* C */ 2, /* D */ 3, /* E */ 3, /* F */ 4,
/* 0 */ 0, /* 1 */ 1, /* 2 */ 1, /* 3 */ 2,
/* 4 */ 1, /* 5 */ 2, /* 6 */ 2, /* 7 */ 3,
/* 8 */ 1, /* 9 */ 2, /* A */ 2, /* B */ 3,
/* C */ 2, /* D */ 3, /* E */ 3, /* F */ 4
);
Now since pshufb works on 4 bit indexes, we will do it in 2 stages. We will first calculate the pop count of the lower 4 bits and then the higher 4 bits.
const __m256i lower_mask = _mm256_set1_epi8(0x0f);
__m256i low_m = _mm256_and_si256(a, lower_mask); // here a is the register contains 8 ints;
// low_m now contains the bytes as their lower 4 bits only.
// high_m contains the higher 4 bits of all bytes
__m256i high_m = _mm256_and_si256(
_mm256_srli_epi16(a, 4),
lower_mask
);
We are calculating the population count of each byte, by calculating for the upper and lower 4 bits
separately since pshufb takes 4 bit indexes.
So low_m here has the lower 4 bits for each byte and high_m has the upper 4 bits.
Now we calculate the actual number of 1's of each 4 bit number through the lookup table.
__m256i pl = _mm256_shuffle_epi8(lookup, low_m);
__m256i ph = _mm256_shuffle_epi8(lookup, high_m);
The low_m and high_m act as indices and lookup table returns their population count.
The pl and ph now contain the population count of the 4-bit indices.
The sum can be found as:
__m256i s = _mm256_add_epi8(pl, ph);
Now we divide the 32 bytes into 4 groups of 8 byte each.
__m256i sum = _mm256_sad_epu8(
s, _mm256_setzero_si256()
);
uint64_t *arr = (uint64_t *)∑
int ans = arr[0] + arr[1] + arr[2] + arr[3];
printf("%d\n", ans);
Now we have our final answer in ans variable. You can extend this to more than 8 binary strings by looping over them, taking 8 at a time.