I'm trying to write code (in C) to sample from a negative binomial distribution parameterized by $r,p$, where $r$ is not necessarily integer (also called Polya distribution).
I've found a number of explanations, and it looks like I could sample $y$ from a Gamma(r,1) distribution, transform $y$ somehow, and then use it to sample from a Poisson distribution, but that means I need to find a Gamma distribution with non-integer $r$.
All of this is a bit scary, and I was wondering if anyone could point me to some source code to do it. I don't much care what language. I have access to gamma functions, which I'm pretty sure will be needed.