Set bits randomly

makefile

all: a.out b.out

a.out: randbit.cpp walker.hpp
	g++ -O3  $< -o $@

b.out: randbit.cpp walker.hpp
	g++ -O3 -DELIMINATE_BRANCH $< -o $@

clean:
	rm -f a.out b.out test.dat result.png

randbit.cpp

//------------------------------------------------------------------------
#include <cstdio>
#include <random>
#include <iostream>
#include <algorithm>
#include <functional>
#include <chrono>
#include "walker.hpp"
//#define ELIMINATE_BRANCH
//------------------------------------------------------------------------
// 時間計測とチェック
//------------------------------------------------------------------------
template <class T>
void
measure(const char *name, T f, std::vector<double> &data) {
  std::mt19937 mt(2);
  const int N = 4000000;
  const double ninv = 1.0 / static_cast<double>(N);
  const auto s = std::chrono::system_clock::now();
  double ave = 0.0;
  for (int i = 0; i < N; i++) {
    unsigned int s = f(mt);
    int n = __builtin_popcount(s);
    data[n] += ninv;
    ave += n * ninv;
  }
  const auto e = std::chrono::system_clock::now();
  const auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count();
  std::cerr << name << elapsed << "[ms]" << " average = " << ave << std::endl;
}
//------------------------------------------------------------------------
// 二項係数を求める
//------------------------------------------------------------------------
long int
binom(long int n, long int k) {
  k = std::min(k, n - k);
  return k == 0 ? 1 : binom(n - 1, k - 1) * n / k;
}
//------------------------------------------------------------------------
// 素直に32回乱数をふる
//------------------------------------------------------------------------
struct RandomBitsSimple {
  const double p;
  template <class RNG>
  unsigned int
  operator()(RNG &r) {
    std::uniform_real_distribution<double> ud(0.0, 1.0);
    unsigned int s = 0;
    for (int i = 0; i < 32; i++) {
      if (ud(r) < p) s |= (1 << i);
    }
    return s;
  }
  RandomBitsSimple(const double _p): p(_p) {}
};
//------------------------------------------------------------------------
std::vector<double>
generate_binomial(double p) {
  std::vector<double> v;
  for (int i = 0; i < 32; i++) {
    v.push_back(binom(32, i) * pow(p, i) * pow(1.0 - p, 32 - i));
  }
  return v;
}
//------------------------------------------------------------------------
// 二項分布で何ビット立てるか決めてから
// Robert Floyd' sampling algorithmでビットシャッフル
//------------------------------------------------------------------------
struct RandomBitsShuffle {
  WalkerAlias wa;
  std::uniform_real_distribution<double> ud;

  // 32ビット中mビットをランダムに立てる
  template <class RNG>
  unsigned int
  get_shuffled_bits(int m, RNG &r) {
    const int n = 32;
    unsigned int s = 0;
    for (int j = n - m + 1; j <= n; j++) {
      int k = r() % j;
      unsigned int t = 1 << k;
#ifdef ELIMINATE_BRANCH
      s = s | t | ((t & s) << (j - k - 1) & (1 << (j - 1)));
#else
      if (t & s) {
        s |= 1 << (j - 1);
      } else {
        s |= t;
      }
#endif
    }
    return s;
  }

  template <class RNG>
  unsigned int
  operator()(RNG &r) {
    // Walker's Alias Method で何ビット立てるか決める
    int m = wa.get_index(r);
    if (m > 16) {
      return ~get_shuffled_bits(32 - m, r);
    } else {
      return get_shuffled_bits(m, r);
    }
  }
  RandomBitsShuffle(double p) : wa(generate_binomial(p)) {}
};
//------------------------------------------------------------------------
// 二項分布で何ビット立てるか決めてから
// Robert Floyd' sampling algorithmでビットシャッフル
// 0.625を作ってからの差分
//------------------------------------------------------------------------
std::mt19937 mt(1);
struct RandomBitsShuffleMorita {
  WalkerAlias wa;
  // 32ビット中mビットをランダムに立てる
  template <class RNG>
  unsigned int
  get_shuffled_bits(int m, RNG &r) {
    const int n = 32;
    unsigned int s = 0;
    for (int j = n - m + 1; j <= n; j++) {
      int k = r() % j;
      unsigned int t = 1 << k;
#ifdef ELIMINATE_BRANCH
      s = s | t | ((t & s) << (j - k - 1) & (1 << (j - 1)));
#else
      if (t & s) {
        s |= 1 << (j - 1);
      } else {
        s |= t;
      }
#endif
    }
    return s;
  }

  template <class RNG>
  unsigned int
  operator()(RNG &r) {
    // Walker's Alias Method で何ビット立てるか決める
    int m = wa.get_index(r);
    unsigned int s = r() | (r() & r());
    if (m > 16) {
      return s | ~get_shuffled_bits(32 - m, r);
    } else {
      return s | get_shuffled_bits(m, r);
    }
  }
  RandomBitsShuffleMorita(double p) : wa(generate_binomial((p - 0.625) / (1.0 - 0.625))) {}
};
//------------------------------------------------------------------------
std::vector<double>
generate_poisson(double la) {
  std::vector<double> v;
  double x = exp(-la);
  int k = 0;
  v.push_back(x);
  while (k < la  || x > 1e-14) {
    k++;
    x = x * la / k;
    v.push_back(x);
  }
  return v;
}
//------------------------------------------------------------------------
//ポアソン分布の数だけORを取る
//------------------------------------------------------------------------
struct RandomBitsOr {
  WalkerAlias wa;
  template <class RNG>
  unsigned int
  operator()(RNG &r) {
    unsigned int s = 0;
    unsigned int n = wa.get_index(r);
    for (unsigned int i = 0; i < n; i++) {
      s |= (1 << (r() & 31));
    }
    return s;
  }
  RandomBitsOr(double p) : wa(generate_poisson(-32.0 * log((1.0 - p)))) {}
};
//------------------------------------------------------------------------
//ポアソン分布の数だけORを取る (Morita Ver.)
//------------------------------------------------------------------------
struct RandomBitsOrMorita {
  WalkerAlias wa;

  template <class RNG>
  unsigned int
  operator()(RNG &r) {
    unsigned int s = (r() | (r() & r()));
    unsigned int n = wa.get_index(r);
    for (unsigned int i = 0; i < n; i++) {
      s |= (1 << (r() & 31));
    }
    return s;
  }

  RandomBitsOrMorita(double p) : wa(generate_poisson(- 32.0 * log(8.0 / 3.0 * (1.0 - p)))) {}
};
//------------------------------------------------------------------------
int
main(void) {
  const double p = 0.65;
  std::vector<double> data(33, 0);
  std::vector<double> data2(33, 0);
  std::vector<double> data3(33, 0);
  std::vector<double> data4(33, 0);
  std::vector<double> data5(33, 0);
  // 厳密解
  std::vector<double> probs;
  for (int i = 0; i < 32; i++) {
    probs.push_back(binom(32, i) * pow(p, i) * pow(1.0 - p, 32 - i));
  }
  measure("simple     ", RandomBitsSimple(p), data);
  measure("shuffle    ", RandomBitsShuffle(p), data2);
  measure("shuffle(M) ", RandomBitsShuffleMorita(p), data3);
  measure("or         ", RandomBitsOr(p), data4);
  measure("or (M)     ", RandomBitsOrMorita(p), data5);

  for (int i = 0; i < 32; i++) {
    std::cout << i << " ";
    std::cout << probs[i] << " "; // Exact
    std::cout << data[i] << " ";// Simple
    std::cout << data2[i] << " "; // Shuffle
    std::cout << data3[i] << " "; // Shuffle (M)
    std::cout << data4[i] << " "; // Or
    std::cout << data5[i] << " "; // Or (M)
    std::cout << std::endl;
  }
}
//------------------------------------------------------------------------

walker.hpp

#pragma once

#include <vector>
#include <numeric>
#include <algorithm>
#include <random>
#include <iostream>

class WalkerAlias {
  private:
    unsigned int n;
    unsigned int nb;
    std::vector<int> index;
    std::vector<double> prob;
    std::vector<unsigned int> probi;
    std::mt19937 mt;
  public:
    WalkerAlias(std::vector<double> a) {
      nb = 1;
      while ((unsigned int)(1 << nb) <= a.size()) {
        nb++;
      }
      n = (1 << nb);
      index.resize(n, 0.0);
      prob.resize(n, 0.0);
      probi.resize(n, 0);
      for (unsigned int i = 0; i < n; i++) {
        prob[i] = a[i];
      }
      double ave = std::accumulate(prob.begin(), prob.end(), 0.0) / prob.size();
      for (auto &v : prob) {
        v /= ave;
      }
      std::vector<int> small, large;
      for (unsigned int i = 0; i < n; i++) {
        if (prob[i] < 1.0) {
          small.push_back(i);
        } else {
          large.push_back(i);
        }
        index[i] = i;
      }
      while (small.size() && large.size()) {
        const int j = small.back();
        small.pop_back();
        const int k = large.back();
        index[j] = k;
        prob[k] = prob[k] - 1.0 + prob[j];
        if (prob[k] < 1.0) {
          small.push_back(k);
          large.pop_back();
        }
      }
      unsigned int int_max = ~static_cast<unsigned int>(0);
      for (unsigned int i = 0; i < n; i++) {
        probi[i] = static_cast<unsigned int>(prob[i] * int_max);
      }
    }
    template<class RNG>
    int
    get_index(RNG &r) {
      int k = r() & (n - 1);
      if (r() > probi[k]) {
        return index[k];
      } else {
        return k;
      }
    }
};

$ g++ -O3 randbit.cpp
$ ./a.out  > test.dat
simple		1617[ms] average = 20.801
shuffle		509[ms] average = 20.8027
shuffle (M) 	175[ms] average = 20.7994
or 		408[ms] average = 20.7993
or (M)		148[ms] average = 20.7973