'Gamma distribution.  Not the gamma function!
Conditions on the parameters are alpha > 0 and beta > 0.
The probability distribution function is:
x ** (alpha - 1) * math.exp(-x / beta)
pdf(x) =  --------------------------------------
math.gamma(alpha) * beta ** alpha'

def gammavariate(self, alpha, beta):

    if ((alpha <= 0.0) or (beta <= 0.0)):
        raise ValueError, 'gammavariate:    alpha    and    beta    must    be    >    0.0'
    random = self.random
    if (alpha > 1.0):
        ainv = _sqrt(((2.0 * alpha) - 1.0))
        bbb = (alpha - LOG4)
        ccc = (alpha + ainv)
        while 1:
            u1 = random()
            if (not (1e-07 < u1 < 0.9999999)):
                continue
            u2 = (1.0 - random())
            v = (_log((u1 / (1.0 - u1))) / ainv)
            x = (alpha * _exp(v))
            z = ((u1 * u1) * u2)
            r = ((bbb + (ccc * v)) - x)
            if ((((r + SG_MAGICCONST) - (4.5 * z)) >= 0.0) or (r >= _log(z))):
                return (x * beta)
    else:
        if (alpha == 1.0):
            u = random()
            while (u <= 1e-07):
                u = random()
            return ((- _log(u)) * beta)
        while 1:
            u = random()
            b = ((_e + alpha) / _e)
            p = (b * u)
            if (p <= 1.0):
                x = (p ** (1.0 / alpha))
            else:
                x = (- _log(((b - p) / alpha)))
            u1 = random()
            if (p > 1.0):
                if (u1 <= (x ** (alpha - 1.0))):
                    break
            elif (u1 <= _exp((- x))):
                break
        return (x * beta)
