I'm looking for C/C++ implementations of functions that return random variates multinomial and Dirichlet distributions. This is in the context of a calculation for posterior predictive p-values, part of which includes a MCMC step. I've been using Python, but am rewriting parts of my code in C++ for speed reasons. I've been using the implementations from numpy.random thus far.
My current options include
Pull the implementations from numpy.random, which are most likely written in C.
Use the R versions of these functions, which are conveniently available for Debian in the r-mathlib library. This library has been my default choice for such things for a long time, because the R people know their stuff when it comes to probabilistic computing.
However, I'd be open to other versions of these functions. Suggestions?
EDIT: r-mathlib doesn't appear to have a function for sampling from the Dirichet distribution, though according to Wikipedia's entry on the Dirichlet distribution I can sum Gamma variates. I wonder if that is a good way of doing this.
EDIT2: If possible, please comment on why you think the implementation you suggest is a good choice.