What's the fastest method to pack an array of booleans to integers?

mul_vs_bitset.cpp

#include <vector>
#include <bitset>
#include <array>
#include <random>
#include <algorithm>
#include <cstring>
#include <chrono>
#include <iostream>

void pack_booleans(std::vector<uint8_t>& dest, const std::vector<bool>& src)
{
  // How many bits to process at a time.
  // Tests on multiple platforms suggest that 1024 provides the most speed up.
  constexpr uint64_t bit_stride = 1024;
  constexpr uint64_t byte_stride = bit_stride / 8;
  constexpr uint64_t magic = 0x8040201008040201;

  auto a = std::array<uint8_t, bit_stride>();
  auto t = std::array<uint64_t, byte_stride>();
  size_t dest_idx = 0;
  auto itr_finish = src.cbegin() + (src.size() / bit_stride) * bit_stride;

  for (auto itr = src.cbegin(); itr != itr_finish; itr += bit_stride) {
    std::copy(itr, itr + bit_stride, a.begin());
    std::memcpy(t.data(), a.data(), a.size());
    std::transform(t.cbegin(), t.cend(), dest.begin() + dest_idx,
                   [](auto e) { return (magic * e) >> 56; }); 
    dest_idx += byte_stride;
  }
}

void use_bitset(std::vector<unsigned long>& dest, const std::vector<bool>& src)
{
  auto bset = std::bitset<64>();
  size_t idx = 0;

  for (size_t i = 0; i < src.size(); i += 64) {
    for (size_t j = 0; j < 64; j++)
      bset[j] = src[i + j];
    dest[idx++] = bset.to_ulong();
  }
}

int main(int argc, char* argv[])
{
  const size_t N = 2048ul * 1024 * 1024;
  auto src = std::vector<bool>(N, false);

  std::random_device rd;  //Will be used to obtain a seed for the random number engine
  std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd()
  std::uniform_int_distribution<uint64_t> distrib(0, N-1);
 
  for (size_t n=0; n<1000; ++n)
    src[distrib(gen)] = true;

  auto dest1 = std::vector<uint8_t>(N / 8); 
  auto start1 = std::chrono::steady_clock::now();
  pack_booleans(dest1, src);  
  auto end1 = std::chrono::steady_clock::now();
  auto millis1 = std::chrono::duration_cast<std::chrono::milliseconds>(end1 - start1).count();

  auto dest2 = std::vector<unsigned long>(N / 64);
  auto start2 = std::chrono::steady_clock::now();
  use_bitset(dest2, src);
  auto end2 = std::chrono::steady_clock::now();
  auto millis2 = std::chrono::duration_cast<std::chrono::milliseconds>(end2 - start2).count();

  std::cout << "multiplication takes time " << millis1 << ", bitset takes time " << millis2 << std::endl;

  std::cout << "random output 1: " << +dest1[distrib(gen) / 8] << std::endl;
  std::cout << "random output 2: " << dest2[distrib(gen) / 64] << std::endl;
}