From 55d2d40a05261ce8a89d9373cc8749bba73319a7 Mon Sep 17 00:00:00 2001 From: Marco Barbone Date: Thu, 27 Aug 2026 10:03:44 -0400 Subject: [PATCH 1/2] fix: use std::popcount in C++20, not only in C++23 xsimd_common_bit.hpp gated the standard path on XSIMD_CPP_VERSION > 202002L, so exactly C++20 skipped std::popcount and compiled the backport instead. Every other XSIMD_CPP_VERSION test in the tree uses >=. The static_assert in test_bit.cpp is the check: std::popcount is constexpr and the backport is not, so a C++20 build that stops reaching the standard library fails to compile. Assisted-by: Claude Opus 5 --- include/xsimd/arch/common/xsimd_common_bit.hpp | 2 +- test/test_bit.cpp | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/include/xsimd/arch/common/xsimd_common_bit.hpp b/include/xsimd/arch/common/xsimd_common_bit.hpp index 607bddf1c..b811da869 100644 --- a/include/xsimd/arch/common/xsimd_common_bit.hpp +++ b/include/xsimd/arch/common/xsimd_common_bit.hpp @@ -7,7 +7,7 @@ #include "../../config/xsimd_config.hpp" -#if XSIMD_CPP_VERSION > 202002L +#if XSIMD_CPP_VERSION >= 202002L #include diff --git a/test/test_bit.cpp b/test/test_bit.cpp index de3b60c27..e00fa433e 100644 --- a/test/test_bit.cpp +++ b/test/test_bit.cpp @@ -14,6 +14,15 @@ #include "test_utils.hpp" +// std::popcount is constexpr and the backport in xsimd_common_bit.hpp is not, +// so this fails to compile if the C++20 build stops reaching the standard one. +#if XSIMD_CPP_VERSION >= 202002L +#include +#if __cpp_lib_bitops >= 201907L +static_assert(xsimd::detail::popcount(0xffu) == 8, "C++20 must use std::popcount"); +#endif +#endif + template struct bit_test { From 15211880b11ae95b384488ea4e7eaf6c4f8959c7 Mon Sep 17 00:00:00 2001 From: Marco Barbone Date: Thu, 13 Aug 2026 10:43:25 -0400 Subject: [PATCH 2/2] feat: add per-lane popcount on batch Count the bits set in each element of an integer batch. The common kernel is the SWAR fold; x86 uses the PSHUFB nibble lookup from SSSE3 up, NEON uses CNT with pairwise widening adds, SVE uses svcnt_x, and WASM uses i8x16.popcnt with pairwise widening extends, folding 32-bit counts with a shift-and-add pair for 64-bit elements. VSX and VXE use vec_popcnt, which the compiler maps to a single VPOPCNTB/H/W/D or VPOPCT. The scalar overload forwards to detail::popcount and returns T, so it agrees with the batch overload. The common kernel runs the whole fold on 64-bit lanes whatever T is, and masks the surviving byte at the end. Every step of the fold confines its own carries, so a bit that crosses a T boundary is masked off again, and a target with no narrow shift does not pay for one to be emulated. This drops emulated<128> u8 from 356 instructions to 20, and measures 35.4x there, 12.6x on emulated<128> u16 and 2.60x on avx2 u8. repeat_pattern() builds the fold's masks from the pattern and the type of the repeat unit, so 0x5555555555555555 reads as repeat_pattern() rather than as a literal. For 64-bit elements the two x86 nibble tables carry a +4 and a -4 bias, after libpopcnt, so PSADBW yields the byte count and the 8-byte sum in one instruction. This drops the VPADDB, and measures 1.09x on SSE and AVX2 and 1.05x on AVX-512. sse2 gets its own kernel, since the SWAR fold is what a default x86-64 build runs and it leaves PSADBW and PMADDWD unused. PSADBW sums the eight byte counts of a 64-bit element in one instruction, PMADDWD adds the two halves of a 32-bit one, and the byte shifts drop the mask that the fold already applies. This cuts 64-bit elements from 25 instructions to 16, 32-bit ones from 22 to 19 and 8-bit ones from 15 to 14, and measures 1.5x, 1.09x and 1.08x. 16-bit elements keep the fold's instruction count, since SSE2 has no PMADDUBSW to combine the two byte counts of an element. avx512vnni gains a 32-bit kernel: VPDPBUSD does in one uop what the VPMADDUBSW and VPMADDWD pair does in two, and the zero accumulator is free because the register copy is eliminated at rename. This measures 1.14x. The same substitution on 256-bit vectors is neutral, since three ports serve them, so avxvnni gets no kernel. The CI job labelled avx512vnni built for knm, which enables avx5124vnniw rather than avx512vnni and selected the avx512pf arch, so it covered neither kernel. It now builds for cascadelake. Assisted-by: Claude Opus 5 --- .github/workflows/linux.yml | 3 +- docs/source/api/bitwise_operators_index.rst | 2 + .../arch/common/xsimd_common_bitwise.hpp | 87 +++++++++++++++++++ include/xsimd/arch/xsimd_avx.hpp | 9 ++ include/xsimd/arch/xsimd_avx2.hpp | 32 +++++++ include/xsimd/arch/xsimd_avx512bw.hpp | 43 +++++++++ .../xsimd/arch/xsimd_avx512vnni_avx512bw.hpp | 21 +++++ .../arch/xsimd_avx512vnni_avx512vbmi2.hpp | 21 +++++ 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 | 32 +++++++ 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 | 14 +++ test/test_batch_int.cpp | 65 ++++++++++++++ test/test_bit.cpp | 11 +++ test/test_xsimd_api.cpp | 37 ++++++++ 22 files changed, 492 insertions(+), 1 deletion(-) create mode 100644 include/xsimd/arch/common/xsimd_common_bitwise.hpp diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 6fe7779a7..2a6893f50 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -116,7 +116,8 @@ jobs: CMAKE_EXTRA_ARGS="$CMAKE_EXTRA_ARGS -DTARGET_ARCH=icelake-server" fi if [[ '${{ matrix.sys.flags }}' == 'avx512vnni' ]]; then - CMAKE_EXTRA_ARGS="$CMAKE_EXTRA_ARGS -DTARGET_ARCH=knm" + # knm enables avx5124vnniw, not avx512vnni: it selects the avx512pf arch + CMAKE_EXTRA_ARGS="$CMAKE_EXTRA_ARGS -DTARGET_ARCH=cascadelake" fi if [[ '${{ matrix.sys.flags }}' == 'i386' ]]; then CXX_FLAGS="$CXX_FLAGS -m32" diff --git a/docs/source/api/bitwise_operators_index.rst b/docs/source/api/bitwise_operators_index.rst index e4631dcb3..0dc654879 100644 --- a/docs/source/api/bitwise_operators_index.rst +++ b/docs/source/api/bitwise_operators_index.rst @@ -48,6 +48,8 @@ Bitwise Operators +---------------------------------------+----------------------------------------------------+ | :cpp:func:`rotl` | per slot rotate left | +---------------------------------------+----------------------------------------------------+ +| :cpp:func:`popcount` | per slot population count | ++---------------------------------------+----------------------------------------------------+ ---- 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..834fc212e --- /dev/null +++ b/include/xsimd/arch/common/xsimd_common_bitwise.hpp @@ -0,0 +1,87 @@ +/*************************************************************************** + * 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_details.hpp" + +#include +#include +#include +#include + +namespace xsimd +{ + + namespace kernel + { + + using namespace types; + + namespace detail + { + // Pattern P repeated across U from bit 0 up. The type of P is the + // repeat unit, so repeat_pattern() is + // 0x5555555555555555 and repeat_pattern() + // is 0x0f0f0f0f. + template + constexpr U repeat_pattern() noexcept + { + using unit = decltype(P); + static_assert(std::is_unsigned::value, "the repeat unit must be unsigned"); + static_assert(P != 0, "the pattern must be non-zero"); + static_assert(sizeof(U) % sizeof(unit) == 0, "the repeat unit must divide U"); + U result = 0; + for (std::size_t i = 0; i < sizeof(U) / sizeof(unit); ++i) + result = U(result | U(U(P) << (i * sizeof(unit) * CHAR_BIT))); + return result; + } + } + + // popcount + // SWAR fold, popcount64b from Hacker's Delight, listed in + // https://en.wikipedia.org/wiki/Hamming_weight#Efficient_implementation + template >*/> + XSIMD_INLINE batch popcount(batch const& self, requires_arch) noexcept + { + using U = as_unsigned_integer_t; + 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(detail::repeat_pattern())); + w_type const m2(detail::repeat_pattern()); + x = (x & m2) + ((x >> 2) & m2); + x = (x + (x >> 4)) & w_type(detail::repeat_pattern()); + 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(U(0xff))); + } + } +} + +#endif diff --git a/include/xsimd/arch/xsimd_avx.hpp b/include/xsimd/arch/xsimd_avx.hpp index 84b1ba3a0..8ae82c24f 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..e10c84010 100644 --- a/include/xsimd/arch/xsimd_avx2.hpp +++ b/include/xsimd/arch/xsimd_avx2.hpp @@ -962,6 +962,38 @@ 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, 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..f44d8742f 100644 --- a/include/xsimd/arch/xsimd_avx512bw.hpp +++ b/include/xsimd/arch/xsimd_avx512bw.hpp @@ -561,6 +561,49 @@ 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, 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..64295f593 100644 --- a/include/xsimd/arch/xsimd_avx512vnni_avx512bw.hpp +++ b/include/xsimd/arch/xsimd_avx512vnni_avx512bw.hpp @@ -13,5 +13,26 @@ #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 it needs + // costs nothing, since the register copy is eliminated at rename + 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..9c9ce45cb 100644 --- a/include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp +++ b/include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp @@ -13,5 +13,26 @@ #define XSIMD_AVX512VNNI_AVX512VBMI2_HPP #include "../types/xsimd_avx512vnni_avx512vbmi2_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 it needs + // costs nothing, since the register copy is eliminated at rename + return _mm512_dpbusd_epi32(_mm512_setzero_si512(), detail::popcount_bytes(self), _mm512_set1_epi8(1)); + } + } +} #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 b247b2bd6..746b6aa51 100644 --- a/include/xsimd/arch/xsimd_common_fwd.hpp +++ b/include/xsimd/arch/xsimd_common_fwd.hpp @@ -85,6 +85,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 c0aff00bb..84f576478 100644 --- a/include/xsimd/arch/xsimd_neon.hpp +++ b/include/xsimd/arch/xsimd_neon.hpp @@ -3680,6 +3680,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..5c6d0e3df 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 std::enable_if_t, T> + popcount(T x) noexcept + { + return T(detail::popcount(std::make_unsigned_t(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 eccba3b36..1682cb477 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); + if constexpr (sizeof(T) == 1) + return counts; + // PSADBW sums the eight byte counts of a 64-bit element in one + // instruction, which the shift-and-mask fold needs nine to do + 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..5ef870a4c 100644 --- a/include/xsimd/arch/xsimd_ssse3.hpp +++ b/include/xsimd/arch/xsimd_ssse3.hpp @@ -82,6 +82,38 @@ 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, 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 dc36bcf00..1bd5f31d5 100644 --- a/include/xsimd/arch/xsimd_sve.hpp +++ b/include/xsimd/arch/xsimd_sve.hpp @@ -1206,6 +1206,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..b3686e279 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_i64x2_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..a5f8a3a8e 100644 --- a/include/xsimd/types/xsimd_api.hpp +++ b/include/xsimd/types/xsimd_api.hpp @@ -1946,6 +1946,20 @@ 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 integer values. + * @return per slot population count. + */ + template + XSIMD_INLINE batch popcount(batch const& x) noexcept + { + 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..6aec94fc1 100644 --- a/test/test_batch_int.cpp +++ b/test/test_batch_int.cpp @@ -303,6 +303,66 @@ struct batch_int_test t.run(); } + // 0, ~0, single bits, prefixes, suffixes and pseudo-random words: the + // patterns the bit ops actually branch on + 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); + CHECK_BATCH_EQ(xsimd::popcount(batch_type::load_unaligned(in.data())), expected); + } + } + void test_less_than_underflow() const { batch_type test_negative_compare = batch_type(5) - 6; @@ -361,5 +421,10 @@ TEST_CASE_TEMPLATE("[batch int tests]", B, BATCH_INT_TYPES) { Test.test_less_than_underflow(); } + + SUBCASE("popcount") + { + Test.test_popcount(); + } } #endif diff --git a/test/test_bit.cpp b/test/test_bit.cpp index e00fa433e..b86191207 100644 --- a/test/test_bit.cpp +++ b/test/test_bit.cpp @@ -223,6 +223,17 @@ struct bit_test } }; +// repeat_pattern is a constexpr mask table, so the check is the compile +// itself: a wrong value or a wrong repeat stride fails to build. +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x55, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x3333, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x0f0f0f0f, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x0f0f0f0f, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x5555555555555555ULL, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x3333333333333333ULL, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0x0f0f0f0f0f0f0f0fULL, "repeat_pattern"); +static_assert(xsimd::kernel::detail::repeat_pattern() == 0xffffffffffffffffULL, "repeat_pattern"); + TEST_CASE_TEMPLATE("[bit operations]", T, uint8_t, uint16_t, uint32_t, uint64_t) { diff --git a/test/test_xsimd_api.cpp b/test/test_xsimd_api.cpp index bc8d8f6bd..fd133b348 100644 --- a/test/test_xsimd_api.cpp +++ b/test/test_xsimd_api.cpp @@ -413,6 +413,38 @@ struct xsimd_api_integral_types_functions CHECK_EQ(extract(xsimd::mod(T(val0), T(val1))), val0 % val1); } + // 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() + { + constexpr int bits = std::numeric_limits::digits + std::numeric_limits::is_signed; + using U = std::make_unsigned_t; + for (int i = 0; i < bits; ++i) + { + value_type const single = value_type(U(U(1) << i)); + value_type const prefix = value_type(U(U(~U(0)) << i)); + value_type const suffix = value_type(U(~U(U(~U(0)) << i))); + INFO("popcount, bit " << i); + CHECK_EQ(extract(xsimd::popcount(T(single))), naive_popcount(single)); + CHECK_EQ(extract(xsimd::popcount(T(prefix))), naive_popcount(prefix)); + CHECK_EQ(extract(xsimd::popcount(T(suffix))), naive_popcount(suffix)); + } + for (value_type v : { value_type(0), value_type(U(~U(0))), value_type(0x5a), value_type(0x3c) }) + { + INFO("popcount, value " << int64_t(v)); + CHECK_EQ(extract(xsimd::popcount(T(v))), naive_popcount(v)); + } + } + void test_rotl() { constexpr auto N = std::numeric_limits::digits + std::numeric_limits::is_signed; @@ -479,6 +511,11 @@ TEST_CASE_TEMPLATE("[xsimd api | integral types functions]", B, INTEGRAL_TYPES) Test.test_mod(); } + SUBCASE("popcount") + { + Test.test_popcount(); + } + SUBCASE("rotl") { Test.test_rotl();