gamma distribution f=x^(alpha-1)*exp(-x/theta)/(gamma(alpha)*theta^alpha)
alpha has to be bigger than 1, for alpha<1 use gammaD(alpha)=gammaD(alpha+1)*pow(r.uniform!(T),1/alpha)
from Marsaglia and Tsang, ACM Transaction on Mathematical Software, Vol. 26, N. 3
2000, p 363-372
gamma distribution f=x^(alpha-1)*exp(-x/theta)/(gamma(alpha)*theta^alpha) alpha has to be bigger than 1, for alpha<1 use gammaD(alpha)=gammaD(alpha+1)*pow(r.uniform!(T),1/alpha) from Marsaglia and Tsang, ACM Transaction on Mathematical Software, Vol. 26, N. 3 2000, p 363-372