This should be faster than the gamma variate transformation we are currently
using, and it seems to work better for parameters smaller than one.

The algorithm is also used by the R language, however I did not consult their
implementation in order to avoid licensing problems.

Reference:

R. C. H. Cheng (1978).
Generating beta variates with nonintegral shape parameters.
Communications of the ACM 21, 317-322.
https://doi.org/10.1145/359460.359482