I've looked at your code but I fail to create a test case. I would like to compare the result of your code with other tools like R, octave. For instance using this https://svn.r-project.org/R/trunk/src/nmath/rgamma.c adapted accordingly in order to remove R dependencies.
The code is a bit cryptic, but is it possible that the 2 different algorithms for a<1 and a>1 are for efficiency reasons? This is probably obvious to you but I'll say it anway : I wouldn't expect, nor is it required, that the generated values agree but their sample distribution, however, should be the same.
int main() { typedef double value_type; typedef boost::mt19937 urng_type;
value_type shape = 1.0; value_type scale = 3.0; unsigned int seed = 5489u;
std::cout << "Using Boost.Random ..." << ::std::endl; urng_type urng(seed); std::cout << rgamma_rnd(shape, scale, urng) << std::endl; std::cout << rgamma_rnd(shape, scale, urng) << std::endl; std::cout << rgamma_rnd(shape, scale, urng) << std::endl;
std::cout << "Using Boost.Random + Boost.Statistic ..." << ::std::endl; urng = urng_type(seed); std::cout << rgamma_math(shape, scale, urng) << std::endl; std::cout << rgamma_math(shape, scale, urng) << std::endl; std::cout << rgamma_math(shape, scale, urng) << std::endl; }
I ran these and they agree (PS: I think you meant to pass the generator by reference, not value).