This is what I use at the moment.

It is the GSL function modified to work with CuRAND.

**device** double ran_gamma (curandState localState, const double a, const double b){

/* assume a > 0 */

if (a < 1){

double u = curand_uniform_double(&localState);

return ran_gamma (localState, 1.0 + a, b) * pow (u, 1.0 / a);

}

{

double x, v, u;

double d = a - 1.0 / 3.0;

double c = (1.0 / 3.0) / sqrt (d);

```
while (1){
do{
x = curand_normal_double(&localState);
v = 1.0 + c * x;
} while (v <= 0);
v = v * v * v;
u = curand_uniform_double(&localState);
if (u < 1 - 0.0331 * x * x * x * x)
break;
if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
break;
}
return b * d * v;
```

}

}