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.

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.

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 <noreply@anthropic.com>

Author: DiamonDinoia

.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"

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                          |
++---------------------------------------+----------------------------------------------------+
 
 ----
 

include/xsimd/arch/common/xsimd_common_bitwise.hpp

@@ -0,0 +1,70 @@
+/***************************************************************************
+ * 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 <climits>
+#include <cstddef>
+
+namespace xsimd
+{
+
+    namespace kernel
+    {
+
+        using namespace types;
+
+        namespace detail
+        {
+            // bit i is set when i / s is even: 0x55.. for s == 1, 0x33.. for
+            // s == 2, 0x0f0f.. for s == 4, and so on up to s == bits / 2
+            template <class U>
+            constexpr U alternating_mask(unsigned s) noexcept
+            {
+                U m = 0;
+                for (unsigned i = 0; i < sizeof(U) * CHAR_BIT; ++i)
+                    if (((i / s) & 1u) == 0u)
+                        m = U(m | U(U(1) << i));
+                return m;
+            }
+        }
+
+        // popcount
+        // SWAR fold, popcount64b from Hacker's Delight, listed in
+        // https://en.wikipedia.org/wiki/Hamming_weight#Efficient_implementation
+        template <class A, class T, class /*=std::enable_if_t<std::is_integral_v<T>>*/>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<common>) noexcept
+        {
+            using U = as_unsigned_integer_t<T>;
+            using b_type = batch<U, A>;
+            constexpr std::size_t bits = sizeof(T) * CHAR_BIT;
+
+            b_type x = bitwise_cast<U>(self);
+            x = x - ((x >> 1) & b_type(detail::alternating_mask<U>(1)));
+            b_type const m2(detail::alternating_mask<U>(2));
+            x = (x & m2) + ((x >> 2) & m2);
+            x = (x + (x >> 4)) & b_type(detail::alternating_mask<U>(4));
+            if constexpr (bits >= 16)
+                x = (x + (x >> 8)) & b_type(detail::alternating_mask<U>(8));
+            if constexpr (bits >= 32)
+                x = (x + (x >> 16)) & b_type(detail::alternating_mask<U>(16));
+            if constexpr (bits >= 64)
+                x = (x + (x >> 32)) & b_type(detail::alternating_mask<U>(32));
+            return bitwise_cast<T>(x);
+        }
+    }
+}
+
+#endif

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 <class A, class T, class = std::enable_if_t<std::is_integral_v<T>>>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<avx>) noexcept
+        {
+            return detail::fwd_to_sse([](__m128i s) noexcept
+                                      { return popcount(batch<T, ssse3>(s), ssse3 {}); },
+                                      self);
+        }
+
         // reciprocal
         template <class A>
         XSIMD_INLINE batch<float, A> reciprocal(batch<float, A> const& self,

include/xsimd/arch/xsimd_avx2.hpp

@@ -962,6 +962,38 @@ namespace xsimd
                                                { return batch<uint64_t, A>(_mm256_mul_epu32(a, b)); });
         }
 
+        // popcount
+        template <class A, class T, class = std::enable_if_t<std::is_integral_v<T>>>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<avx2>) noexcept
+        {
+            // per-byte counts from a nibble lookup indexed by VPSHUFB, after
+            // Wojciech Muła, http://0x80.pl/articles/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 <class A, class T, class = std::enable_if_t<std::is_integral_v<T>>>
         XSIMD_INLINE T reduce_add(batch<T, A> const& self, requires_arch<avx2>) noexcept

include/xsimd/arch/xsimd_avx512bw.hpp

@@ -561,6 +561,49 @@ namespace xsimd
             return detail::compare_int_avx512bw<A, T, _MM_CMPINT_NE>(self, other);
         }
 
+        namespace detail
+        {
+            // per-byte counts from a nibble lookup indexed by VPSHUFB, after
+            // Wojciech Muła, http://0x80.pl/articles/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 <class A, class T, class = std::enable_if_t<std::is_integral_v<T>>>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<avx512bw>) 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 <class A, class T, class = std::enable_if_t<std::is_integral_v<T>>>
         XSIMD_INLINE batch<T, A> sadd(batch<T, A> const& self, batch<T, A> const& other, requires_arch<avx512bw>) noexcept

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 <class A, class T, class = std::enable_if_t<std::is_integral_v<T> && sizeof(T) == 4>>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<avx512vnni<avx512bw>>) 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

include/xsimd/arch/xsimd_avx512vnni_avx512vbmi2.hpp

@@ -13,5 +13,27 @@
 #define XSIMD_AVX512VNNI_AVX512VBMI2_HPP
 
 #include "../types/xsimd_avx512vnni_avx512vbmi2_register.hpp"
+#include "./xsimd_avx512bw.hpp"
+#include "./xsimd_avx512vbmi2.hpp"
+
+namespace xsimd
+{
+
+    namespace kernel
+    {
+
+        using namespace types;
+
+        // popcount
+        template <class A, class T, class = std::enable_if_t<std::is_integral_v<T> && sizeof(T) == 4>>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<avx512vnni<avx512vbmi2>>) 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

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"

include/xsimd/arch/xsimd_common_fwd.hpp

@@ -85,6 +85,8 @@ namespace xsimd
         XSIMD_INLINE batch<T, A> rotr(batch<T, A> const& self, STy other, requires_arch<common>) noexcept;
         template <size_t count, class A, class T>
         XSIMD_INLINE batch<T, A> rotr(batch<T, A> const& self, requires_arch<common>) noexcept;
+        template <class A, class T, class = std::enable_if_t<std::is_integral_v<T>>>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<common>) noexcept;
         template <class A, class T>
         XSIMD_INLINE batch<T, A> load(T const* mem, aligned_mode, requires_arch<A>) noexcept;
         template <class A, class T>