Draw a number [0...n) from a loaded dice (Vose's alias method)

loadeddice.cc

//
// Draw a number [0...n) from a loaded dice (Vose's alias method)
//
#include <iostream>
#include <vector>

using namespace std;

unsigned int xor32() {
  static unsigned long y = 2463534242UL;
  y^=(y<<13); y^=(y>>17);
  return (y^=(y<<15));
}
double urand() {
  return xor32() / (1.0 + UINT_MAX);
}

struct dice {
  vector<double> p;
  vector<int> a;
  dice(const vector<double> &p_) : p(p_) { 
    int n = p.size(); a.resize(n);
    vector<int> S, L;
    for (int i = 0; i < n; ++i) {
      p[i] *= n;
      if (p[i] < 1) S.push_back(i); 
      else          L.push_back(i);
    }
    while (!S.empty() && !L.empty()) {
      int s = S.back(); S.pop_back();
      int l = L.back(); L.pop_back();
      p[l] += p[s] - 1; a[s] = l;
      if (p[l] < 1) S.push_back(l); 
      else          L.push_back(l);
    }
    for (auto s: S) p[s] = 1.0;
    for (auto l: L) p[l] = 1.0;
  }
  int draw() {
    int k = p.size() * urand();
    return urand() < p[k] ? k : a[k];
  }
};


int main() {
  dice d({1.0/10.0, 2.0/10.0, 3.0/10.0, 4.0/10.0});

  int m = 1000000;
  int count[4] = {0};
  for (int i = 0; i < m; ++i) {
    count[d.draw()] += 1;
  }
  for (int i = 0; i < 4; ++i) {
    cout << 1.0 * count[i] / m << endl;
  }
}