This isn't very efficient when a is large (there are better methods), but it's easy to implement.
Generate a random variates uniformly distributed over [0,1). Call them u1, u2, ... ua.
Compute Li = -log(1-Ui) for i = 1 to a. Each of the Li has the Exp(1) distribution.
Let S be the sum of all the Li.
Return S/b as the gamma(a,b) random variate.