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.

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.

The NEON kernel takes UDOT when __ARM_FEATURE_DOTPROD is defined: one
instruction replaces the two pairwise widening adds that fold byte counts into
32-bit elements. The Arm cross-compilation matrix gains an armv8.2-a+dotprod
row, so qemu covers both paths.

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/cross-arm.yml

@@ -13,6 +13,7 @@ jobs:
           - { platform: 'arm',     arch: 'armv7-a', dir: 'arm-linux-gnueabihf', flags: '-mfpu=neon', full: 'ON'}
           - { platform: 'arm',     arch: 'armv7-a', dir: 'arm-linux-gnueabihf', flags: '-mfpu=vfpv3-d16', full: 'OFF' } # no neon
           - { platform: 'aarch64', arch: 'armv8-a', dir: 'aarch64-linux-gnu', flags: '', full: 'ON' }
+          - { platform: 'aarch64', arch: 'armv8.2-a+dotprod', dir: 'aarch64-linux-gnu', flags: '', full: 'OFF' }
         sys:
           - { compiler: 'gcc',   version: '10' }
           - { compiler: 'gcc',   version: '14' }

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