From a97b94e1ab259ea64e5cadfa1a333d42ef029b19 Mon Sep 17 00:00:00 2001 From: Marco Barbone Date: Fri, 9 Oct 2026 15:02:24 -0400 Subject: [PATCH] feat: add per-lane popcount for unsigned integral batches SWAR fold via detail::repeat_pattern() masks (Wikipedia Hamming_weight), with per-arch dispatch (SSE2/SSSE3/AVX2/AVX512BW/VNNI, NEON vcnt, WASM via unsigned shift), scalar fallback std::popcount / __builtin_popcount. Batch popcount restricted to unsigned integral, bool excluded. Tests: per-bit oracle, all-lanes-to-N-bits, batch-vs-scalar consistency. --- docs/source/api/bitwise_operators_index.rst | 3 + .../xsimd/arch/common/xsimd_common_bit.hpp | 143 +++++++++--------- .../arch/common/xsimd_common_bitwise.hpp | 66 ++++++++ include/xsimd/arch/xsimd_avx.hpp | 9 ++ include/xsimd/arch/xsimd_avx2.hpp | 34 +++++ include/xsimd/arch/xsimd_avx512bw.hpp | 45 ++++++ .../xsimd/arch/xsimd_avx512vnni_avx512bw.hpp | 22 +++ .../arch/xsimd_avx512vnni_avx512vbmi2.hpp | 20 +++ include/xsimd/arch/xsimd_common.hpp | 1 + include/xsimd/arch/xsimd_common_fwd.hpp | 2 + include/xsimd/arch/xsimd_neon.hpp | 25 +++ include/xsimd/arch/xsimd_scalar.hpp | 9 ++ include/xsimd/arch/xsimd_sse2.hpp | 33 ++++ include/xsimd/arch/xsimd_ssse3.hpp | 34 +++++ include/xsimd/arch/xsimd_sve.hpp | 8 + include/xsimd/arch/xsimd_vsx.hpp | 9 ++ include/xsimd/arch/xsimd_vxe.hpp | 9 ++ include/xsimd/arch/xsimd_wasm.hpp | 20 +++ include/xsimd/types/xsimd_api.hpp | 15 ++ test/test_batch_int.cpp | 111 ++++++++++++++ test/test_bit.cpp | 56 ++++--- test/test_xsimd_api.cpp | 55 +++++++ 22 files changed, 632 insertions(+), 97 deletions(-) create mode 100644 include/xsimd/arch/common/xsimd_common_bitwise.hpp diff --git a/docs/source/api/bitwise_operators_index.rst b/docs/source/api/bitwise_operators_index.rst index e4631dcb3..4934c3bec 100644 --- a/docs/source/api/bitwise_operators_index.rst +++ b/docs/source/api/bitwise_operators_index.rst @@ -48,6 +48,9 @@ Bitwise Operators +---------------------------------------+----------------------------------------------------+ | :cpp:func:`rotl` | per slot rotate left | +---------------------------------------+----------------------------------------------------+ +| :cpp:func:`popcount` | unsigned integral types only; per slot population | +| | count, returned as the same unsigned type | ++---------------------------------------+----------------------------------------------------+ ---- diff --git a/include/xsimd/arch/common/xsimd_common_bit.hpp b/include/xsimd/arch/common/xsimd_common_bit.hpp index 607bddf1c..e559d03d0 100644 --- a/include/xsimd/arch/common/xsimd_common_bit.hpp +++ b/include/xsimd/arch/common/xsimd_common_bit.hpp @@ -7,33 +7,19 @@ #include "../../config/xsimd_config.hpp" -#if XSIMD_CPP_VERSION > 202002L +#include +#include +#include +#if XSIMD_CPP_VERSION >= 202002L #include +#endif -#if __cpp_lib_bitops >= 201907L - +#if XSIMD_CPP_VERSION >= 202002L && __cpp_lib_bitops >= 201907L #include - -namespace xsimd -{ - namespace detail - { - using std::countl_one; - using std::countl_zero; - using std::countr_one; - using std::countr_zero; - using std::popcount; - } -} - +#define XSIMD_HAS_STD_BITOPS 1 #endif -#else - -#include -#include - #ifdef __has_builtin #define XSIMD_HAS_BUILTIN(x) __has_builtin(x) #else @@ -48,81 +34,92 @@ namespace xsimd { namespace detail { - // FIXME: We could do better by dispatching to the appropriate popcount instruction - // depending on the arch. + // Lane value made of \c byte repeated across every byte of \c T. + template + XSIMD_INLINE constexpr T repeat_pattern(unsigned char byte) noexcept + { + T out = 0; + for (std::size_t i = 0; i < sizeof(T); ++i) + out = T((out << CHAR_BIT) | byte); + return out; + } - template >> + // Portable word-parallel popcount. Bits fold to per-byte counts, then + // a single multiply sums the byte counts into the top byte. W is at + // least unsigned int, otherwise the final multiply loses the sum for + // narrow T. + // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel + template + XSIMD_INLINE constexpr int popcount_swar(T v) noexcept + { + static_assert(std::is_unsigned_v, "popcount requires an unsigned integral type"); + using W = std::conditional_t<(sizeof(T) < sizeof(unsigned int)), unsigned int, T>; + W x = W(v); + x = x - ((x >> 1) & repeat_pattern(0x55)); + x = (x & repeat_pattern(0x33)) + ((x >> 2) & repeat_pattern(0x33)); + x = (x + (x >> 4)) & repeat_pattern(0x0f); + return int((x * W(~W(0) / 255)) >> ((sizeof(W) - 1) * CHAR_BIT)); + } + + // Dispatch order, first match wins: SWAR for GCC on x86 without + // POPCNT, std::popcount (C++20), __builtin_popcountg (GCC 14, + // Clang 19), sized builtins, MSVC intrinsics, SWAR. The first rung + // beats std::popcount because GCC lowers it, like every builtin, to + // a libgcc call that is slower than the inline fold. The builtins + // fold to constant expressions, the MSVC intrinsics do not. +#if defined(_MSC_VER) && !defined(__clang__) && !defined(XSIMD_HAS_STD_BITOPS) + template XSIMD_INLINE int popcount(T x) noexcept +#else + template + XSIMD_INLINE constexpr int popcount(T x) noexcept +#endif { -#if XSIMD_HAS_BUILTIN(__builtin_popcountg) + static_assert(std::is_unsigned_v, "popcount requires an unsigned integral type"); +#if defined(__GNUC__) && !defined(__clang__) && (defined(__i386__) || defined(__x86_64__)) && !defined(__POPCNT__) + return popcount_swar(x); +#elif defined(XSIMD_HAS_STD_BITOPS) + return std::popcount(x); +#elif XSIMD_HAS_BUILTIN(__builtin_popcountg) return __builtin_popcountg(x); -#else - if constexpr (sizeof(T) == 1) +#elif XSIMD_HAS_BUILTIN(__builtin_popcount) && XSIMD_HAS_BUILTIN(__builtin_popcountll) + if constexpr (sizeof(T) <= 4) { -#if XSIMD_HAS_BUILTIN(__builtin_popcount) return __builtin_popcount(x); -#elif defined(_MSC_VER) - return __popcnt(x); -#else - // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSet64 - return ((uint64_t)x * 0x200040008001ULL & 0x111111111111111ULL) % 0xf; -#endif } - else if constexpr (sizeof(T) == 2) + else { -#if XSIMD_HAS_BUILTIN(__builtin_popcount) - return __builtin_popcount(x); + return __builtin_popcountll(x); + } #elif defined(_MSC_VER) + if constexpr (sizeof(T) == 2) + { return __popcnt16(x); -#else - // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSet64 - constexpr unsigned long long msb12 = 0x1001001001001ULL; - constexpr unsigned long long mask5 = 0x84210842108421ULL; - - unsigned int v = (unsigned int)x; - - return ((v & 0xfff) * msb12 & mask5) % 0x1f - + (((v & 0xfff000) >> 12) * msb12 & mask5) % 0x1f; -#endif } - else if constexpr (sizeof(T) == 4) + else if constexpr (sizeof(T) == 1 || sizeof(T) == 4) { -#if XSIMD_HAS_BUILTIN(__builtin_popcount) - return __builtin_popcount(x); -#elif defined(_MSC_VER) return __popcnt(x); -#else - // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel - x = x - ((x >> 1) & (T) ~(T)0 / 3); - x = (x & (T) ~(T)0 / 15 * 3) + ((x >> 2) & (T) ~(T)0 / 15 * 3); - x = (x + (x >> 4)) & (T) ~(T)0 / 255 * 15; - return (x * ((T) ~(T)0 / 255)) >> (sizeof(T) - 1) * CHAR_BIT; -#endif } else { // sizeof(T) == 8 -#if XSIMD_HAS_BUILTIN(__builtin_popcountll) - return __builtin_popcountll(x); -#elif XSIMD_HAS_BUILTIN(__builtin_popcount) - return __builtin_popcount((unsigned int)x) + __builtin_popcount((unsigned int)(x >> 32)); -#elif defined(_MSC_VER) #ifdef _M_X64 return (int)__popcnt64(x); #else return (int)(__popcnt((unsigned int)x) + __popcnt((unsigned int)(x >> 32))); -#endif -#else - // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel - x = x - ((x >> 1) & (T) ~(T)0 / 3); - x = (x & (T) ~(T)0 / 15 * 3) + ((x >> 2) & (T) ~(T)0 / 15 * 3); - x = (x + (x >> 4)) & (T) ~(T)0 / 255 * 15; - return (x * ((T) ~(T)0 / 255)) >> (sizeof(T) - 1) * CHAR_BIT; #endif } +#else + return popcount_swar(x); #endif } +#ifdef XSIMD_HAS_STD_BITOPS + using std::countl_one; + using std::countl_zero; + using std::countr_one; + using std::countr_zero; +#else template >> XSIMD_INLINE int countl_zero(T x) noexcept { @@ -224,9 +221,11 @@ namespace xsimd { return countr_zero(T(~x)); } - +#endif } } -#endif +#undef XSIMD_HAS_STD_BITOPS +#undef XSIMD_HAS_BUILTIN + #endif diff --git a/include/xsimd/arch/common/xsimd_common_bitwise.hpp b/include/xsimd/arch/common/xsimd_common_bitwise.hpp new file mode 100644 index 000000000..11c4832a8 --- /dev/null +++ b/include/xsimd/arch/common/xsimd_common_bitwise.hpp @@ -0,0 +1,66 @@ +/*************************************************************************** + * Copyright (c) Johan Mabille, Sylvain Corlay, Wolf Vollprecht and * + * Martin Renou * + * Copyright (c) QuantStack * + * Copyright (c) Serge Guelton * + * Copyright (c) Marco Barbone * + * * + * Distributed under the terms of the BSD 3-Clause License. * + * * + * The full license is in the file LICENSE, distributed with this software. * + ****************************************************************************/ + +#ifndef XSIMD_COMMON_BITWISE_HPP +#define XSIMD_COMMON_BITWISE_HPP + +#include "./xsimd_common_bit.hpp" +#include "./xsimd_common_details.hpp" + +#include +#include +#include + +namespace xsimd +{ + + namespace kernel + { + + using namespace types; + + // popcount + // SWAR fold on 64-bit lanes, the efficient implementation from + // https://en.wikipedia.org/wiki/Hamming_weight#Efficient_implementation + template + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + using w_type = batch; + constexpr std::size_t bits = sizeof(T) * CHAR_BIT; + + // Every step runs on 64-bit lanes whatever T is. Each step confines + // its own carries, so a bit that crosses a T boundary is always + // masked off again, and targets with no narrow shift do not pay for + // one to be emulated. + w_type x = bitwise_cast(self); + x = x - ((x >> 1) & w_type(xsimd::detail::repeat_pattern(0x55))); + w_type const m2(xsimd::detail::repeat_pattern(0x33)); + x = (x & m2) + ((x >> 2) & m2); + x = (x + (x >> 4)) & w_type(xsimd::detail::repeat_pattern(0x0f)); + if constexpr (bits == 8) + return bitwise_cast(x); + + // Byte counts are at most 8, so a per-lane sum is at most 64 and no + // byte carries into the next. The final mask keeps the one byte per + // lane that holds the whole count and drops the bytes that summed + // across a lane boundary. + x = x + (x >> 8); + if constexpr (bits >= 32) + x = x + (x >> 16); + if constexpr (bits >= 64) + x = x + (x >> 32); + return bitwise_cast(x) & batch(T(0xff)); + } + } +} + +#endif diff --git a/include/xsimd/arch/xsimd_avx.hpp b/include/xsimd/arch/xsimd_avx.hpp index 49bd12b90..28fde929a 100644 --- a/include/xsimd/arch/xsimd_avx.hpp +++ b/include/xsimd/arch/xsimd_avx.hpp @@ -1411,6 +1411,15 @@ namespace xsimd return _mm256_castps_si256(_mm256_xor_ps(_mm256_castsi256_ps(self.data), _mm256_castsi256_ps(other.data))); } + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + return detail::fwd_to_sse([](__m128i s) noexcept + { return popcount(batch(s), ssse3 {}); }, + self); + } + // reciprocal template XSIMD_INLINE batch reciprocal(batch const& self, diff --git a/include/xsimd/arch/xsimd_avx2.hpp b/include/xsimd/arch/xsimd_avx2.hpp index ba6825cb8..75e7a9f2b 100644 --- a/include/xsimd/arch/xsimd_avx2.hpp +++ b/include/xsimd/arch/xsimd_avx2.hpp @@ -962,6 +962,40 @@ namespace xsimd { return batch(_mm256_mul_epu32(a, b)); }); } + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // per-byte counts from a nibble lookup indexed by VPSHUFB, after + // Wojciech Muła, http://0x80.pl/notesen/2008-05-24-sse-popcount.html + __m256i const low_mask = _mm256_set1_epi8(0x0f); + __m256i const lo = _mm256_and_si256(self, low_mask); + __m256i const hi = _mm256_and_si256(_mm256_srli_epi16(self, 4), low_mask); + if constexpr (sizeof(T) == 8) + { + // tables biased by +4 and -4 turn the VPSADBW difference into the + // per-byte count: |(c_lo+4) - (4-c_hi)| = c_lo + c_hi, and VPSADBW + // sums |a-b| over each block of eight bytes, so one instruction + // adds the nibble counts and sums the eight bytes, + // after https://github.com/kimwalisch/libpopcnt + __m256i const lookup_lo = _mm256_broadcastsi128_si256(_mm_setr_epi8(4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8)); + __m256i const lookup_hi = _mm256_broadcastsi128_si256(_mm_setr_epi8(4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0)); + return _mm256_sad_epu8(_mm256_shuffle_epi8(lookup_lo, lo), _mm256_shuffle_epi8(lookup_hi, hi)); + } + else + { + __m256i const lookup = _mm256_broadcastsi128_si256(_mm_setr_epi8(0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4)); + __m256i const counts = _mm256_add_epi8(_mm256_shuffle_epi8(lookup, lo), _mm256_shuffle_epi8(lookup, hi)); + // wider elements sum their byte counts + if constexpr (sizeof(T) == 1) + return counts; + else if constexpr (sizeof(T) == 2) + return _mm256_maddubs_epi16(counts, _mm256_set1_epi8(1)); + else + return _mm256_madd_epi16(_mm256_maddubs_epi16(counts, _mm256_set1_epi8(1)), _mm256_set1_epi16(1)); + } + } + // reduce_add template >> XSIMD_INLINE T reduce_add(batch const& self, requires_arch) noexcept diff --git a/include/xsimd/arch/xsimd_avx512bw.hpp b/include/xsimd/arch/xsimd_avx512bw.hpp index 2d32002b9..584274bf5 100644 --- a/include/xsimd/arch/xsimd_avx512bw.hpp +++ b/include/xsimd/arch/xsimd_avx512bw.hpp @@ -561,6 +561,51 @@ namespace xsimd return detail::compare_int_avx512bw(self, other); } + namespace detail + { + // per-byte counts from a nibble lookup indexed by VPSHUFB, after + // Wojciech Muła, http://0x80.pl/notesen/2008-05-24-sse-popcount.html + XSIMD_INLINE __m512i popcount_bytes(__m512i self) noexcept + { + __m512i const low_mask = _mm512_set1_epi8(0x0f); + __m512i const lookup = _mm512_broadcast_i32x4(_mm_setr_epi8(0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4)); + __m512i const lo = _mm512_and_si512(self, low_mask); + __m512i const hi = _mm512_and_si512(_mm512_srli_epi16(self, 4), low_mask); + return _mm512_add_epi8(_mm512_shuffle_epi8(lookup, lo), _mm512_shuffle_epi8(lookup, hi)); + } + } + + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + if constexpr (sizeof(T) == 8) + { + __m512i const low_mask = _mm512_set1_epi8(0x0f); + __m512i const lo = _mm512_and_si512(self, low_mask); + __m512i const hi = _mm512_and_si512(_mm512_srli_epi16(self, 4), low_mask); + // tables biased by +4 and -4 turn the VPSADBW difference into the + // per-byte count: |(c_lo+4) - (4-c_hi)| = c_lo + c_hi, and VPSADBW + // sums |a-b| over each block of eight bytes, so one instruction + // adds the nibble counts and sums the eight bytes, + // after https://github.com/kimwalisch/libpopcnt + __m512i const lookup_lo = _mm512_broadcast_i32x4(_mm_setr_epi8(4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8)); + __m512i const lookup_hi = _mm512_broadcast_i32x4(_mm_setr_epi8(4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0)); + return _mm512_sad_epu8(_mm512_shuffle_epi8(lookup_lo, lo), _mm512_shuffle_epi8(lookup_hi, hi)); + } + else + { + __m512i const counts = detail::popcount_bytes(self); + // wider elements sum their byte counts + if constexpr (sizeof(T) == 1) + return counts; + else if constexpr (sizeof(T) == 2) + return _mm512_maddubs_epi16(counts, _mm512_set1_epi8(1)); + else + return _mm512_madd_epi16(_mm512_maddubs_epi16(counts, _mm512_set1_epi8(1)), _mm512_set1_epi16(1)); + } + } + // sadd template >> XSIMD_INLINE batch sadd(batch const& self, batch const& other, requires_arch) noexcept diff --git a/include/xsimd/arch/xsimd_avx512vnni_avx512bw.hpp b/include/xsimd/arch/xsimd_avx512vnni_avx512bw.hpp index c95069df1..ec3f82817 100644 --- a/include/xsimd/arch/xsimd_avx512vnni_avx512bw.hpp +++ b/include/xsimd/arch/xsimd_avx512vnni_avx512bw.hpp @@ -13,5 +13,27 @@ #define XSIMD_AVX512VNNI_AVX512_BW_HPP #include "../types/xsimd_avx512vnni_avx512bw_register.hpp" +#include "./xsimd_avx512bw.hpp" + +namespace xsimd +{ + + namespace kernel + { + + using namespace types; + + // popcount + template && sizeof(T) == 4>> + XSIMD_INLINE batch popcount(batch const& self, requires_arch>) noexcept + { + // VPDPBUSD does in one uop what the VPMADDUBSW + VPMADDWD pair of + // the avx512bw kernel does in two; the zero accumulator compiles + // to the vpxor zero idiom, which clears dependencies on the old + // register value + return _mm512_dpbusd_epi32(_mm512_setzero_si512(), detail::popcount_bytes(self), _mm512_set1_epi8(1)); + } + } +} #endif diff --git a/include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp b/include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp index 552869d25..8fb641707 100644 --- a/include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp +++ b/include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp @@ -13,5 +13,25 @@ #define XSIMD_AVX512VNNI_AVX512VBMI2_HPP #include "../types/xsimd_avx512vnni_avx512vbmi2_register.hpp" +#include "./xsimd_avx512vnni_avx512bw.hpp" + +namespace xsimd +{ + + namespace kernel + { + + using namespace types; + + // popcount + // avx512vnni derives from avx512vbmi2, not from + // avx512vnni, so dispatch needs this forwarder + template && sizeof(T) == 4>> + XSIMD_INLINE batch popcount(batch const& self, requires_arch>) noexcept + { + return popcount(self, avx512vnni {}); + } + } +} #endif diff --git a/include/xsimd/arch/xsimd_common.hpp b/include/xsimd/arch/xsimd_common.hpp index 1d800e349..37da93b52 100644 --- a/include/xsimd/arch/xsimd_common.hpp +++ b/include/xsimd/arch/xsimd_common.hpp @@ -14,6 +14,7 @@ #include "./common/xsimd_common_arithmetic.hpp" #include "./common/xsimd_common_bit.hpp" +#include "./common/xsimd_common_bitwise.hpp" #include "./common/xsimd_common_cast.hpp" #include "./common/xsimd_common_complex.hpp" #include "./common/xsimd_common_logical.hpp" diff --git a/include/xsimd/arch/xsimd_common_fwd.hpp b/include/xsimd/arch/xsimd_common_fwd.hpp index 139b0b282..4f7e25294 100644 --- a/include/xsimd/arch/xsimd_common_fwd.hpp +++ b/include/xsimd/arch/xsimd_common_fwd.hpp @@ -89,6 +89,8 @@ namespace xsimd XSIMD_INLINE batch rotr(batch const& self, STy other, requires_arch) noexcept; template XSIMD_INLINE batch rotr(batch const& self, requires_arch) noexcept; + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept; template XSIMD_INLINE batch load(T const* mem, aligned_mode, requires_arch) noexcept; template diff --git a/include/xsimd/arch/xsimd_neon.hpp b/include/xsimd/arch/xsimd_neon.hpp index 4521f38da..860cf001d 100644 --- a/include/xsimd/arch/xsimd_neon.hpp +++ b/include/xsimd/arch/xsimd_neon.hpp @@ -3682,6 +3682,31 @@ namespace xsimd WRAP_MASK_OP(countr_one) #undef WRAP_MASK_OP + + /************ + * popcount * + ************/ + + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // CNT counts bytes; wider elements fold them with pairwise widening + // adds, after Wojciech Muła's popcnt_neon_vcnt, + // https://github.com/WojciechMula/sse-popcount + uint8x16_t counts = vcntq_u8(bitwise_cast(self).data); + if constexpr (sizeof(T) == 1) + return bitwise_cast(batch(counts)); + else if constexpr (sizeof(T) == 2) + return bitwise_cast(batch(vpaddlq_u8(counts))); + else + { + uint32x4_t const wide = vpaddlq_u16(vpaddlq_u8(counts)); + if constexpr (sizeof(T) == 4) + return bitwise_cast(batch(wide)); + else + return bitwise_cast(batch(vpaddlq_u32(wide))); + } + } } } diff --git a/include/xsimd/arch/xsimd_scalar.hpp b/include/xsimd/arch/xsimd_scalar.hpp index 6adc9dce5..f76f837d1 100644 --- a/include/xsimd/arch/xsimd_scalar.hpp +++ b/include/xsimd/arch/xsimd_scalar.hpp @@ -13,6 +13,7 @@ #define XSIMD_SCALAR_HPP #include "../config/xsimd_macros.hpp" +#include "./common/xsimd_common_bit.hpp" #include #include @@ -408,6 +409,14 @@ namespace xsimd return +x; } + // returns T rather than int, so that the scalar and the batch overload agree + template + XSIMD_INLINE T popcount(T x) noexcept + { + static_assert(std::is_unsigned_v && !std::is_same_v, "popcount requires an unsigned integral type"); + return T(detail::popcount(x)); + } + XSIMD_INLINE float reciprocal(float const& x) noexcept { return 1.f / x; diff --git a/include/xsimd/arch/xsimd_sse2.hpp b/include/xsimd/arch/xsimd_sse2.hpp index 58e35254d..3fcc3acd4 100644 --- a/include/xsimd/arch/xsimd_sse2.hpp +++ b/include/xsimd/arch/xsimd_sse2.hpp @@ -1577,6 +1577,39 @@ namespace xsimd return _mm_xor_pd(self, other); } + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // SWAR fold to per-byte counts, after popcount64b from Hacker's + // Delight, vectorized as in Wojciech Muła's popcnt_SSE_bit_parallel, + // https://github.com/WojciechMula/sse-popcount. The byte shifts need + // no mask of their own, since every bit they leak across a byte + // boundary falls under one of the three masks the fold already + // applies. + __m128i const m1 = _mm_set1_epi8(0x55); + __m128i const m2 = _mm_set1_epi8(0x33); + __m128i const m4 = _mm_set1_epi8(0x0f); + __m128i x = _mm_sub_epi8(self, _mm_and_si128(_mm_srli_epi16(self, 1), m1)); + x = _mm_add_epi8(_mm_and_si128(x, m2), _mm_and_si128(_mm_srli_epi16(x, 2), m2)); + __m128i const counts = _mm_and_si128(_mm_add_epi8(x, _mm_srli_epi16(x, 4)), m4); + // PSADBW sums the eight byte counts of a 64-bit element in one + // instruction, which the shift-and-mask fold needs nine to do + if constexpr (sizeof(T) == 1) + return counts; + else if constexpr (sizeof(T) == 8) + return _mm_sad_epu8(counts, _mm_setzero_si128()); + else + { + __m128i const pairs = _mm_and_si128(_mm_add_epi8(counts, _mm_srli_epi16(counts, 8)), _mm_set1_epi16(0xff)); + if constexpr (sizeof(T) == 2) + return pairs; + // PMADDWD adds the two halves of a 32-bit element in one instruction + else + return _mm_madd_epi16(pairs, _mm_set1_epi16(1)); + } + } + // reciprocal template XSIMD_INLINE batch reciprocal(batch const& self, diff --git a/include/xsimd/arch/xsimd_ssse3.hpp b/include/xsimd/arch/xsimd_ssse3.hpp index 4451d5bac..3f00c87d5 100644 --- a/include/xsimd/arch/xsimd_ssse3.hpp +++ b/include/xsimd/arch/xsimd_ssse3.hpp @@ -82,6 +82,40 @@ namespace xsimd return detail::extract_pair(self, other, i, std::make_index_sequence()); } + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // per-byte counts from a nibble lookup indexed by PSHUFB, after + // Wojciech Muła, http://0x80.pl/notesen/2008-05-24-sse-popcount.html + __m128i const low_mask = _mm_set1_epi8(0x0f); + __m128i const lo = _mm_and_si128(self, low_mask); + __m128i const hi = _mm_and_si128(_mm_srli_epi16(self, 4), low_mask); + if constexpr (sizeof(T) == 8) + { + // tables biased by +4 and -4 turn the PSADBW difference into the + // per-byte count: |(c_lo+4) - (4-c_hi)| = c_lo + c_hi, and PSADBW + // sums |a-b| over each block of eight bytes, so one instruction + // adds the nibble counts and sums the eight bytes, + // after https://github.com/kimwalisch/libpopcnt + __m128i const lookup_lo = _mm_setr_epi8(4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8); + __m128i const lookup_hi = _mm_setr_epi8(4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0); + return _mm_sad_epu8(_mm_shuffle_epi8(lookup_lo, lo), _mm_shuffle_epi8(lookup_hi, hi)); + } + else + { + __m128i const lookup = _mm_setr_epi8(0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4); + __m128i const counts = _mm_add_epi8(_mm_shuffle_epi8(lookup, lo), _mm_shuffle_epi8(lookup, hi)); + // wider elements sum their byte counts + if constexpr (sizeof(T) == 1) + return counts; + else if constexpr (sizeof(T) == 2) + return _mm_maddubs_epi16(counts, _mm_set1_epi8(1)); + else + return _mm_madd_epi16(_mm_maddubs_epi16(counts, _mm_set1_epi8(1)), _mm_set1_epi16(1)); + } + } + // reduce_add template >> XSIMD_INLINE T reduce_add(batch const& self, requires_arch) noexcept diff --git a/include/xsimd/arch/xsimd_sve.hpp b/include/xsimd/arch/xsimd_sve.hpp index b5a9c9612..a1e70aa6c 100644 --- a/include/xsimd/arch/xsimd_sve.hpp +++ b/include/xsimd/arch/xsimd_sve.hpp @@ -1207,6 +1207,14 @@ namespace xsimd return svscale_x(detail_sve::ptrue(), x, exp); } + // popcount + template = 0> + XSIMD_INLINE batch popcount(const batch& self, requires_arch) noexcept + { + using U = as_unsigned_integer_t; + return bitwise_cast(batch(svcnt_x(detail_sve::ptrue(), self))); + } + } // namespace kernel } // namespace xsimd diff --git a/include/xsimd/arch/xsimd_vsx.hpp b/include/xsimd/arch/xsimd_vsx.hpp index b2361b5ba..54fa88407 100644 --- a/include/xsimd/arch/xsimd_vsx.hpp +++ b/include/xsimd/arch/xsimd_vsx.hpp @@ -536,6 +536,15 @@ namespace xsimd return ~vec_cmpeq(self.data, other.data); } + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // VPOPCNTB/H/W/D count the bits of each element in one instruction + using U = as_unsigned_integer_t; + return bitwise_cast(batch(vec_popcnt(bitwise_cast(self).data))); + } + // reciprocal template XSIMD_INLINE batch reciprocal(batch const& self, diff --git a/include/xsimd/arch/xsimd_vxe.hpp b/include/xsimd/arch/xsimd_vxe.hpp index 61e9bf8a5..b3ac0fcf2 100644 --- a/include/xsimd/arch/xsimd_vxe.hpp +++ b/include/xsimd/arch/xsimd_vxe.hpp @@ -408,6 +408,15 @@ namespace xsimd return vec_mergeh(row[0].data, row[1].data) + vec_mergel(row[0].data, row[1].data); } + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // VPOPCT counts the bits of each element in one instruction + using U = as_unsigned_integer_t; + return bitwise_cast(batch(vec_popcnt(bitwise_cast(self).data))); + } + // reduce_add template XSIMD_INLINE float reduce_add(batch const& self, requires_arch) noexcept diff --git a/include/xsimd/arch/xsimd_wasm.hpp b/include/xsimd/arch/xsimd_wasm.hpp index 87120e307..1d62e29d6 100644 --- a/include/xsimd/arch/xsimd_wasm.hpp +++ b/include/xsimd/arch/xsimd_wasm.hpp @@ -1813,6 +1813,26 @@ namespace xsimd { return wasm_i64x2_shuffle(self, other, 0, 2); } + + // popcount + template >> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + // i8x16.popcnt counts bytes only, and pairwise widening addition + // stops at 32-bit elements; a shift-and-add folds the 32-bit + // counts into 64-bit ones, a mask drops the stray high halves + v128_t counts = wasm_i8x16_popcnt(self); + if constexpr (sizeof(T) == 1) + return counts; + counts = wasm_u16x8_extadd_pairwise_u8x16(counts); + if constexpr (sizeof(T) == 2) + return counts; + counts = wasm_u32x4_extadd_pairwise_u16x8(counts); + if constexpr (sizeof(T) == 4) + return counts; + return wasm_v128_and(wasm_i64x2_add(counts, wasm_u64x2_shr(counts, 32)), + wasm_i64x2_splat(0xffffffff)); + } } } diff --git a/include/xsimd/types/xsimd_api.hpp b/include/xsimd/types/xsimd_api.hpp index 4fbb86a99..340f0029d 100644 --- a/include/xsimd/types/xsimd_api.hpp +++ b/include/xsimd/types/xsimd_api.hpp @@ -1946,6 +1946,21 @@ namespace xsimd return kernel::polar(r, theta, A {}); } + /** + * @ingroup batch_bitwise + * + * Counts the bits set in each element of \c x. + * @param x batch of unsigned integer values. + * @return per slot population count. + */ + template + XSIMD_INLINE batch popcount(batch const& x) noexcept + { + static_assert(std::is_unsigned_v && !std::is_same_v, "popcount requires an unsigned integral type"); + detail::static_check_supported_config(); + return kernel::popcount(x, A {}); + } + /** * @ingroup batch_arithmetic * diff --git a/test/test_batch_int.cpp b/test/test_batch_int.cpp index e8e9b23b3..8d8a5dc11 100644 --- a/test/test_batch_int.cpp +++ b/test/test_batch_int.cpp @@ -15,9 +15,15 @@ #include "test_utils.hpp" #include +#include namespace xsimd { + static_assert(std::is_same_v { })), batch>, "popcount on uint8_t batch returns batch"); + static_assert(std::is_same_v { })), batch>, "popcount on uint16_t batch returns batch"); + static_assert(std::is_same_v { })), batch>, "popcount on uint32_t batch returns batch"); + static_assert(std::is_same_v { })), batch>, "popcount on uint64_t batch returns batch"); + template struct test_int_min_max { @@ -303,6 +309,102 @@ struct batch_int_test t.run(); } + // 0, ~0, single bits, prefixes, suffixes and pseudo-random words + static array_type bit_patterns(size_t seed) + { + constexpr size_t bits = sizeof(value_type) * CHAR_BIT; + using U = std::make_unsigned_t; + array_type a; + for (size_t i = 0; i < size; ++i) + { + size_t k = seed * size + i; + size_t sh = k % bits; + U u; + switch (k % 6) + { + case 0: + u = U(0); + break; + case 1: + u = U(~U(0)); + break; + case 2: + u = U(U(1) << sh); + break; + case 3: + u = U(~U(0)) << sh; + break; + case 4: + u = U(U(U(1) << sh) - U(1)); + break; + default: + u = U(k * 0x9E3779B9u + 0x7F4A7C15u); + break; + } + a[i] = value_type(u); + } + return a; + } + + // independent oracle: one shift per bit, so the check cannot pass by + // agreeing with the SWAR fold or the lookup tables it is testing + static value_type naive_popcount(value_type v) + { + using U = std::make_unsigned_t; + int n = 0; + for (U u(v); u; u = U(u >> 1)) + n += int(u & U(1)); + return value_type(n); + } + + void test_popcount() const + { + for (size_t s = 0; s < 6; ++s) + { + array_type in = bit_patterns(s), expected; + std::transform(in.cbegin(), in.cend(), expected.begin(), naive_popcount); + INFO("popcount, pattern " << s); + batch_type const in_batch = batch_type::load_unaligned(in.data()); + CHECK_BATCH_EQ(xsimd::popcount(in_batch), expected); + // run the common SWAR kernel directly, so arch overrides cannot hide it + INFO("popcount common kernel, pattern " << s); + CHECK_BATCH_EQ((xsimd::kernel::popcount(in_batch, xsimd::kernel::requires_arch {})), expected); + } + + // one isolated bit per bit position, lanes mixed with 0 and ~0 + constexpr size_t bits = sizeof(value_type) * CHAR_BIT; + using U = std::make_unsigned_t; + for (size_t b = 0; b < bits; ++b) + { + for (size_t r = 0; r < 3; ++r) + { + array_type in, expected; + U const lane[3] = { U(U(1) << b), U(0), U(~U(0)) }; + for (size_t i = 0; i < size; ++i) + in[i] = value_type(lane[(i + r) % 3]); + std::transform(in.cbegin(), in.cend(), expected.begin(), naive_popcount); + INFO("popcount, bit " << b << ", lane offset " << r); + batch_type const in_batch = batch_type::load_unaligned(in.data()); + CHECK_BATCH_EQ(xsimd::popcount(in_batch), expected); + INFO("popcount common kernel, bit " << b << ", lane offset " << r); + CHECK_BATCH_EQ((xsimd::kernel::popcount(in_batch, xsimd::kernel::requires_arch {})), expected); + } + } + + // every lane holds the same value with bits - k bits set: a count + // spilling across a lane boundary raises a neighbour above bits - k + for (size_t k = 0; k < bits; ++k) + { + array_type in, expected; + for (size_t i = 0; i < size; ++i) + in[i] = value_type(U(U(~U(0)) >> k)); + std::transform(in.cbegin(), in.cend(), expected.begin(), naive_popcount); + INFO("popcount, all lanes to " << bits - k << " bits"); + batch_type const in_batch = batch_type::load_unaligned(in.data()); + CHECK_BATCH_EQ(xsimd::popcount(in_batch), expected); + } + } + void test_less_than_underflow() const { batch_type test_negative_compare = batch_type(5) - 6; @@ -361,5 +463,14 @@ TEST_CASE_TEMPLATE("[batch int tests]", B, BATCH_INT_TYPES) { Test.test_less_than_underflow(); } + + // popcount is defined on unsigned types only, much like std::popcount. + if constexpr (std::is_unsigned_v) + { + SUBCASE("popcount") + { + Test.test_popcount(); + } + } } #endif diff --git a/test/test_bit.cpp b/test/test_bit.cpp index de3b60c27..c05bbf902 100644 --- a/test/test_bit.cpp +++ b/test/test_bit.cpp @@ -14,6 +14,17 @@ #include "test_utils.hpp" +// std::popcount is constexpr and so are __builtin_popcount and +// __builtin_popcountll on GCC and Clang, so the backport popcount is a +// constant expression on the MSVC and fallback paths too wherever a builtin +// is available. +#if XSIMD_CPP_VERSION >= 202002L +#include +#endif +#if (XSIMD_CPP_VERSION >= 202002L && __cpp_lib_bitops >= 201907L) || defined(__GNUC__) || defined(__clang__) +static_assert(xsimd::detail::popcount(0xffu) == 8, "popcount must be usable in a constant expression"); +#endif + template struct bit_test { @@ -22,26 +33,23 @@ struct bit_test void test_popcount() { - // Zero - CHECK_EQ(xsimd::detail::popcount(T(0)), 0); - - // All bits set - CHECK_EQ(xsimd::detail::popcount(T(~T(0))), bits::value); - - // Single bit patterns - all should have popcount of 1 - for (int i = 0; i < bits::value; ++i) + auto check = [&](T v, int ref) { - T value = T(T(1) << i); - INFO("popcount(1 << " << i << ")"); - CHECK_EQ(xsimd::detail::popcount(value), 1); - } + INFO("popcount(0x" << std::hex << (unsigned long long)v << std::dec << ")"); + CHECK_EQ(xsimd::detail::popcount(v), ref); +#if !(XSIMD_CPP_VERSION >= 202002L && __cpp_lib_bitops >= 201907L) + // Run the portable fallback directly: on a compiler with a + // builtin the dispatch never reaches it. + CHECK_EQ(xsimd::detail::popcount_swar(v), ref); +#endif + }; - // Powers of 2 minus 1 - known popcounts - for (int i = 1; i < bits::value; ++i) + check(T(0), 0); + check(T(~T(0)), bits::value); + for (int i = 0; i < bits::value; ++i) { - T value = T((T(1) << i) - 1); - INFO("popcount((1 << " << i << ") - 1)"); - CHECK_EQ(xsimd::detail::popcount(value), i); + check(T(T(1) << i), 1); + check(T(T(~T(0)) >> i), bits::value - i); } // Alternating patterns @@ -54,17 +62,15 @@ struct bit_test pattern_aa |= T(0xAA) << (i * 8); pattern_55 |= T(0x55) << (i * 8); } - INFO("popcount(0xAA...)"); - CHECK_EQ(xsimd::detail::popcount(pattern_aa), bits::value / 2); - INFO("popcount(0x55...)"); - CHECK_EQ(xsimd::detail::popcount(pattern_55), bits::value / 2); + check(pattern_aa, bits::value / 2); + check(pattern_55, bits::value / 2); } // Specific test cases - CHECK_EQ(xsimd::detail::popcount(T(1)), 1); - CHECK_EQ(xsimd::detail::popcount(T(3)), 2); - CHECK_EQ(xsimd::detail::popcount(T(7)), 3); - CHECK_EQ(xsimd::detail::popcount(T(15)), 4); + check(T(1), 1); + check(T(3), 2); + check(T(7), 3); + check(T(15), 4); } void test_countl_zero() diff --git a/test/test_xsimd_api.cpp b/test/test_xsimd_api.cpp index bc8d8f6bd..ba3a04d18 100644 --- a/test/test_xsimd_api.cpp +++ b/test/test_xsimd_api.cpp @@ -14,6 +14,8 @@ #include +#include + template struct scalar_type { @@ -41,11 +43,14 @@ bool extract(xsimd::batch_bool const& batch) { return batch.get(0); } #define INTEGRAL_TYPES_HEAD char, unsigned char, signed char, short, unsigned short, int, unsigned int, long, unsigned long #ifdef XSIMD_NO_SUPPORTED_ARCHITECTURE #define INTEGRAL_TYPES_TAIL +#define UNSIGNED_TYPES #else #define INTEGRAL_TYPES_TAIL , xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch +#define UNSIGNED_TYPES , xsimd::batch, xsimd::batch, xsimd::batch, xsimd::batch #endif #define INTEGRAL_TYPES INTEGRAL_TYPES_HEAD INTEGRAL_TYPES_TAIL +#define UNSIGNED_INTEGRAL_TYPES unsigned char, unsigned short, unsigned int, unsigned long UNSIGNED_TYPES // @@ -413,6 +418,43 @@ struct xsimd_api_integral_types_functions CHECK_EQ(extract(xsimd::mod(T(val0), T(val1))), val0 % val1); } + void test_popcount() + { + constexpr int bits = std::numeric_limits::digits + std::numeric_limits::is_signed; + using U = std::make_unsigned_t; + + // Check every lane of the result against std::bitset; the scalar case + // must not build a batch, some CI jobs have no SIMD architecture. + auto check = [&](value_type v) + { + INFO("popcount, value " << int64_t(v)); + const auto expected = value_type(std::bitset(U(v)).count()); + if constexpr (std::is_integral_v) + { + CHECK_EQ(xsimd::popcount(v), expected); + } + else + { + auto got = xsimd::popcount(T(v)); + std::array lanes; + got.store_unaligned(lanes.data()); + for (std::size_t l = 0; l < lanes.size(); ++l) + CHECK_EQ(lanes[l], expected); + } + }; + + for (int i = 0; i < bits; ++i) + { + check(value_type(U(U(1) << i))); + check(value_type(U(U(~U(0)) << i))); + check(value_type(U(~U(U(~U(0)) << i)))); + } + check(value_type(0)); + check(value_type(U(~U(0)))); + check(value_type(0x5a)); + check(value_type(0x3c)); + } + void test_rotl() { constexpr auto N = std::numeric_limits::digits + std::numeric_limits::is_signed; @@ -500,6 +542,19 @@ TEST_CASE_TEMPLATE("[xsimd api | integral types functions]", B, INTEGRAL_TYPES) } } +// The batch popcount API accepts unsigned types only, much like std::popcount. +TEST_CASE_TEMPLATE("[xsimd api | unsigned integral types functions]", B, UNSIGNED_INTEGRAL_TYPES) +{ + using test_type = xsimd_api_integral_types_functions; + + test_type Test; + + SUBCASE("popcount") + { + Test.test_popcount(); + } +} + // Subtracting the type minimum must saturate, not wrap. Previously this was // computed as sadd(x, -min), and -min is not representable. TEST_CASE_TEMPLATE("[xsimd api | ssub at type minimum]", B, INTEGRAL_TYPES)