feat: add per-lane popcount on batch

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

include/xsimd/arch/common/xsimd_common_bit.hpp

@@ -7,33 +7,19 @@
 
 #include "../../config/xsimd_config.hpp"
 
-#if XSIMD_CPP_VERSION > 202002L
+#include <climits>
+#include <cstddef>
+#include <type_traits>
 
+#if XSIMD_CPP_VERSION >= 202002L
 #include <version>
+#endif
 
-#if __cpp_lib_bitops >= 201907L
-
+#if XSIMD_CPP_VERSION >= 202002L && __cpp_lib_bitops >= 201907L
 #include <bit>
-
-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 <climits>
-#include <type_traits>
-
 #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 <class T>
+        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 <class T, class = std::enable_if_t<std::is_unsigned_v<T>>>
+        // 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 <class T>
+        XSIMD_INLINE constexpr int popcount_swar(T v) noexcept
+        {
+            static_assert(std::is_unsigned_v<T>, "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<W>(0x55));
+            x = (x & repeat_pattern<W>(0x33)) + ((x >> 2) & repeat_pattern<W>(0x33));
+            x = (x + (x >> 4)) & repeat_pattern<W>(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 <class T>
         XSIMD_INLINE int popcount(T x) noexcept
+#else
+        template <class T>
+        XSIMD_INLINE constexpr int popcount(T x) noexcept
+#endif
         {
-#if XSIMD_HAS_BUILTIN(__builtin_popcountg)
+            static_assert(std::is_unsigned_v<T>, "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 <class T, class = std::enable_if_t<std::is_unsigned_v<T>>>
         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

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 <climits>
+#include <cstddef>
+#include <cstdint>
+
+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 <class A, class T, class>
+        XSIMD_INLINE batch<T, A> popcount(batch<T, A> const& self, requires_arch<common>) noexcept
+        {
+            using w_type = batch<uint64_t, A>;
+            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<uint64_t>(self);
+            x = x - ((x >> 1) & w_type(xsimd::detail::repeat_pattern<uint64_t>(0x55)));
+            w_type const m2(xsimd::detail::repeat_pattern<uint64_t>(0x33));
+            x = (x & m2) + ((x >> 2) & m2);
+            x = (x + (x >> 4)) & w_type(xsimd::detail::repeat_pattern<uint64_t>(0x0f));
+            if constexpr (bits == 8)
+                return bitwise_cast<T>(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<T>(x) & batch<T, A>(T(0xff));
+        }
+    }
+}
+
+#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,40 @@ 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: |(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 <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,51 @@ 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: |(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 <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