multivar hypergeometric in rcpp

mvhyper.cpp

// destination: array to fill the drawn "balls"
// source: number of "balls" in each "urn"
// n: number of draws to take
// k: number of "urns"
void rmhyper(int* destination, int const* source, int n, int k){
  int sum, x, y;
  size_t i;
  if(n < 0 || k < 0){Rcpp::stop("Invalid parameters of distribution");}
  
  // total number of "balls"
  for(i = 0, sum = 0; i < k; i++){
    y = source[i];
    if(y < 0){Rcpp::stop("Cannot have a negative number of balls in an urn");}
    sum += y;
  }
  if(n > sum){Rcpp::stop("Distribution undefined for n > sum");}
  
  for(i=0; i<k-1; i++){
    // generate ouput by calling rhyper k-1 times
    y = source[i];
    x = (int)R::rhyper((double)y,(double)sum-y,(double)n);
    n -= x;
    sum -= y;
    destination[i] = x;
  }
  // get the last one
  destination[i] = n;
};

// [[Rcpp::export]]
std::vector<int> test_rmhyper(const std::vector<int> balls, const int n){
  int k = balls.size();
  std::vector<int> out(k,0);
  rmhyper(out.data(),balls.data(),n,k);
  return out;
};