In principle this should be possible by means of the "scaling property": Gamma(shape,scale) ~ scale*Gamma(shape,1)
This would simply translate into boost::mt19937 rng; boost::gamma_distribution<> gamma(shape); // 1-parameter Gamma distribution boost::variate_generator
rvg(rng, gamma); double r = scale*rvg();
However, when shape == 1 the implementation uses the fact that Gamma(1) ~ Exp(1). This is rigth when scale==1 but not when scale != 1 since the exact relation is Gamma(1,scale) ~ Exp(1/scale)
An example of this would be gamma(1,0.5) or gamma(1,2.0)? FYI I ran these kolmogorov-smirnov statistics which seem to converge to zero: gamma(1,0.5) (10,0.0974325) (110,0.0464457) (1110,0.0173321) (11110,0.0108157) (111110,0.00259073) (1111110,0.000503333) gamma(1,2) (10,0.0974325) (110,0.0464457) (1110,0.0173321) (11110,0.0108157) (111110,0.00259073) (1111110,0.000503333) To rule out a fluke I've compared make_random_generator(urng,d1) to d where d and d1 only differ by 0.01 in their shape xor scale, and the statistic is larger by a factor of 10 and 100 respectively, at the last iteration... HTH