public double nextGamma(double shape, double scale) {
		if (shape < 1) {
			// [1]: p. 228, Algorithm GS

			while (true) {
				// Step 1:
				final double u = this.nextUniform(0, 1);
				final double bGS = 1 + shape / FastMath.E;
				final double p = bGS * u;

				if (p <= 1) {
					// Step 2:

					final double x = FastMath.pow(p, 1 / shape);
					final double u2 = this.nextUniform(0.0, 1);

					if (u2 > FastMath.exp(-x)) {
						// Reject
						continue;
					} else {
						return scale * x;
					}
				} else {
					// Step 3:

					final double x = -1 * FastMath.log((bGS - p) / shape);
					final double u2 = this.nextUniform(0, 1);

					if (u2 > FastMath.pow(x, shape - 1)) {
						// Reject
						continue;
					} else {
						return scale * x;
					}
				}
			}
		}

		// Now shape >= 1

		final RandomGenerator generator = this.getRan();
		final double d = shape - 0.333333333333333333;
		final double c = 1.0 / (3 * FastMath.sqrt(d));

		while (true) {
			final double x = generator.nextGaussian();
			final double v = (1 + c * x) * (1 + c * x) * (1 + c * x);

			if (v <= 0) {
				continue;
			}

			final double xx = x * x;
			final double u = this.nextUniform(0, 1);

			// Squeeze
			if (u < 1 - 0.0331 * xx * xx) {
				return scale * d * v;
			}

			if (FastMath.log(u) < 0.5 * xx + d * (1 - v + FastMath.log(v))) {
				return scale * d * v;
			}
		}
	}