c++ 20 std:: popcount

以下建议已合并http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p0553r4.html，并应将其添加到<bit>头。

我希望用法是这样的:

#include <bit>
#include <iostream>

int main() {
    std::cout << std::popcount(0x55) << std::endl;
}

当支持GCC时，我会尝试一下，GCC 9.1.0带有g++-9 -std=c++2a仍然不支持它。

提案说:

标题:< > 命名空间STD { // 25.5.6，计数 模板类T > < conexpr int popcount(T x) noexcept;

and:

模板类T > < conexpr int popcount(T x) noexcept; 约束:T是无符号整数类型(3.9.1 [basic.fundamental])。 返回:x值中的1位数。

std::rotl和std::rotr也被添加来执行循环位旋转:c++中循环移位(旋转)操作的最佳实践

这是一个有助于了解您的微架构的问题。我只是在gcc 4.3.3下用-O3编译的两个变量使用c++内联来计时，以消除函数调用开销，十亿次迭代，保持所有计数的运行总和，以确保编译器不删除任何重要的东西，使用rdtsc计时(精确的时钟周期)。

inline int pop2(unsigned x, unsigned y)
{
    x = x - ((x >> 1) & 0x55555555);
    y = y - ((y >> 1) & 0x55555555);
    x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
    y = (y & 0x33333333) + ((y >> 2) & 0x33333333);
    x = (x + (x >> 4)) & 0x0F0F0F0F;
    y = (y + (y >> 4)) & 0x0F0F0F0F;
    x = x + (x >> 8);
    y = y + (y >> 8);
    x = x + (x >> 16);
    y = y + (y >> 16);
    return (x+y) & 0x000000FF;
}

未经修改的黑客喜悦需要122亿周期。我的并行版本(计算的比特数是它的两倍)的运行周期为13.0千兆周期。在2.4GHz的酷睿双核上，两者总共消耗了10.5秒。在这个时钟频率下，25千兆周期= 10秒多一点，所以我相信我的计时是正确的。

这与指令依赖链有关，这对算法非常不利。通过使用一对64位寄存器，我几乎可以再次将速度提高一倍。事实上，如果我聪明一点，早点加上x+y，我就可以减少一些移位。64位版本做了一些小的调整，结果是相同的，但又增加了一倍的比特数。

对于128位SIMD寄存器，这是另一个因素，SSE指令集通常也有聪明的快捷方式。

没有理由让代码特别透明。该算法界面简单，可在多处在线引用，并能通过全面的单元测试。偶然发现它的程序员甚至可能学到一些东西。这些位操作在机器级别上是非常自然的。

好吧，我决定搁置调整后的64位版本。对于这个sizeof(unsigned long) == 8

inline int pop2(unsigned long x, unsigned long y)
{
    x = x - ((x >> 1) & 0x5555555555555555);
    y = y - ((y >> 1) & 0x5555555555555555);
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333);
    y = (y & 0x3333333333333333) + ((y >> 2) & 0x3333333333333333);
    x = (x + (x >> 4)) & 0x0F0F0F0F0F0F0F0F;
    y = (y + (y >> 4)) & 0x0F0F0F0F0F0F0F0F;
    x = x + y; 
    x = x + (x >> 8);
    x = x + (x >> 16);
    x = x + (x >> 32); 
    return x & 0xFF;
}

这看起来是对的(不过我没有仔细测试)。现在计时结果是10.70亿周期/ 14.1亿周期。后面的数字加起来是1280亿比特，相当于这台机器运行了5.9秒。非并行版本稍微加快了一点，因为我在64位模式下运行，它更喜欢64位寄存器，而不是32位寄存器。

让我们看看这里是否有更多的OOO管道。这有点复杂，所以我实际上测试了一些。每一项单独加起来是64，所有项加起来是256。

inline int pop4(unsigned long x, unsigned long y, 
                unsigned long u, unsigned long v)
{
  enum { m1 = 0x5555555555555555, 
         m2 = 0x3333333333333333, 
         m3 = 0x0F0F0F0F0F0F0F0F, 
         m4 = 0x000000FF000000FF };

    x = x - ((x >> 1) & m1);
    y = y - ((y >> 1) & m1);
    u = u - ((u >> 1) & m1);
    v = v - ((v >> 1) & m1);
    x = (x & m2) + ((x >> 2) & m2);
    y = (y & m2) + ((y >> 2) & m2);
    u = (u & m2) + ((u >> 2) & m2);
    v = (v & m2) + ((v >> 2) & m2);
    x = x + y; 
    u = u + v; 
    x = (x & m3) + ((x >> 4) & m3);
    u = (u & m3) + ((u >> 4) & m3);
    x = x + u; 
    x = x + (x >> 8);
    x = x + (x >> 16);
    x = x & m4; 
    x = x + (x >> 32);
    return x & 0x000001FF;
}

我兴奋了一会儿，但结果是gcc在-O3上玩内联的把戏，尽管我在一些测试中没有使用内联关键字。当我让gcc玩把戏时，对pop4()的十亿次调用需要12.56 gigacycles，但我确定它是将参数折叠为常量表达式。更实际的数字似乎是19.6gc，以实现30%的加速。我的测试循环现在看起来像这样，确保每个参数足够不同，以阻止gcc耍花招。

   hitime b4 = rdtsc(); 
   for (unsigned long i = 10L * 1000*1000*1000; i < 11L * 1000*1000*1000; ++i) 
      sum += pop4 (i,  i^1, ~i, i|1); 
   hitime e4 = rdtsc(); 

2560亿比特加起来在8.17秒内过去了。根据16位表查找的基准测试，3200万比特的计算结果为1.02秒。不能直接比较，因为另一个工作台没有给出时钟速度，但看起来我已经把64KB表版本的鼻涕打出来了，这首先是L1缓存的悲惨使用。

更新:决定做明显的和创建pop6()通过增加四个重复的行。结果是22.8gc, 3840亿比特在9.5秒内加起来。所以还有20%现在是800毫秒，320亿比特。

下面是功能优等递归解决方案，它是迄今为止最纯粹的一个(并且可以用于任何位长度!):

template<typename T>
int popcnt(T n)
{
  if (n>0)
    return n&1 + popcnt(n>>1);
  return 0; 
}

一个简单的算法来计算设置位的数量:

int countbits(n) {
    int count = 0;
    while(n != 0) {
        n = n & (n-1);
        count++;
    }
    return count;
}

以11(1011)为例，尝试手动运行该算法。它应该对你有很大帮助!

对于那些想要在c++ 11中为任何无符号整数类型作为consexpr函数的人(tacklelib/include/tacklelib/utility/math.hpp):

#include <stdint.h>
#include <limits>
#include <type_traits>

const constexpr uint32_t uint32_max = (std::numeric_limits<uint32_t>::max)();

namespace detail
{
    template <typename T>
    inline constexpr T _count_bits_0(const T & v)
    {
        return v - ((v >> 1) & 0x55555555);
    }

    template <typename T>
    inline constexpr T _count_bits_1(const T & v)
    {
        return (v & 0x33333333) + ((v >> 2) & 0x33333333);
    }

    template <typename T>
    inline constexpr T _count_bits_2(const T & v)
    {
        return (v + (v >> 4)) & 0x0F0F0F0F;
    }

    template <typename T>
    inline constexpr T _count_bits_3(const T & v)
    {
        return v + (v >> 8);
    }

    template <typename T>
    inline constexpr T _count_bits_4(const T & v)
    {
        return v + (v >> 16);
    }

    template <typename T>
    inline constexpr T _count_bits_5(const T & v)
    {
        return v & 0x0000003F;
    }

    template <typename T, bool greater_than_uint32>
    struct _impl
    {
        static inline constexpr T _count_bits_with_shift(const T & v)
        {
            return
                detail::_count_bits_5(
                    detail::_count_bits_4(
                        detail::_count_bits_3(
                            detail::_count_bits_2(
                                detail::_count_bits_1(
                                    detail::_count_bits_0(v)))))) + count_bits(v >> 32);
        }
    };

    template <typename T>
    struct _impl<T, false>
    {
        static inline constexpr T _count_bits_with_shift(const T & v)
        {
            return 0;
        }
    };
}

template <typename T>
inline constexpr T count_bits(const T & v)
{
    static_assert(std::is_integral<T>::value, "type T must be an integer");
    static_assert(!std::is_signed<T>::value, "type T must be not signed");

    return uint32_max >= v ?
        detail::_count_bits_5(
            detail::_count_bits_4(
                detail::_count_bits_3(
                    detail::_count_bits_2(
                        detail::_count_bits_1(
                            detail::_count_bits_0(v)))))) :
        detail::_impl<T, sizeof(uint32_t) < sizeof(v)>::_count_bits_with_shift(v);
}

谷歌测试库中的附加测试:

#include <stdlib.h>
#include <time.h>

namespace {
    template <typename T>
    inline uint32_t _test_count_bits(const T & v)
    {
        uint32_t count = 0;
        T n = v;
        while (n > 0) {
            if (n % 2) {
                count += 1;
            }
            n /= 2;
        }
        return count;
    }
}

TEST(FunctionsTest, random_count_bits_uint32_100K)
{
    srand(uint_t(time(NULL)));
    for (uint32_t i = 0; i < 100000; i++) {
        const uint32_t r = uint32_t(rand()) + (uint32_t(rand()) << 16);
        ASSERT_EQ(_test_count_bits(r), count_bits(r));
    }
}

TEST(FunctionsTest, random_count_bits_uint64_100K)
{
    srand(uint_t(time(NULL)));
    for (uint32_t i = 0; i < 100000; i++) {
        const uint64_t r = uint64_t(rand()) + (uint64_t(rand()) << 16) + (uint64_t(rand()) << 32) + (uint64_t(rand()) << 48);
        ASSERT_EQ(_test_count_bits(r), count_bits(r));
    }
}

我在任何地方都没见过这种方法:

int nbits(unsigned char v) {
    return ((((v - ((v >> 1) & 0x55)) * 0x1010101) & 0x30c00c03) * 0x10040041) >> 0x1c;
}

它每字节工作一次，所以对于一个32位整数，它必须被调用四次。它源于横向加法，但它使用两个32位乘法将指令数量减少到只有7条。

大多数当前的C编译器将使用SIMD (SSE2)指令优化这个函数，当请求的数量是4的倍数时，它变得非常有竞争力。它是可移植的，可以定义为宏或内联函数，并且不需要数据表。

这种方法可以扩展为一次处理16位，使用64位乘法。但是，当所有16位都被设置时，它会失败，返回0，所以它只能在0xFFFF输入值不存在时使用。由于64位操作，它也比较慢，并且没有很好地优化。