docs: @Reference(authors = "J. H. Ahrens, U. Dieter", //
        title = "Computer methods for sampling from gamma, beta, Poisson and binomial distributions", //
        booktitle = "Computing 12", //
        url = "https://doi.org/10.1007/BF02293108", //
        bibkey = "DBLP:journals/computing/AhrensD74")
    @Reference(authors = "J. H. Ahrens, U. Dieter", //
        title = "Generating gamma variates by a modified rejection technique", //
        booktitle = "Communications of the ACM 25", //
        url = "https://doi.org/10.1145/358315.358390", //
        bibkey = "DBLP:journals/cacm/AhrensD82")
    public static double nextRandom(double k, double theta, Random random) {
      /* Constants */
      final double q1 = 0.0416666664, q2 = 0.0208333723, q3 = 0.0079849875;
      final double q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332;
      final double q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320;
      final double a1 = 0.333333333, a2 = -0.249999949, a3 = 0.199999867;
      final double a4 = -0.166677482, a5 = 0.142873973, a6 = -0.124385581;
      final double a7 = 0.110368310, a8 = -0.112750886, a9 = 0.104089866;
      final double e1 = 1.000000000, e2 = 0.499999994, e3 = 0.166666848;
      final double e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826;
      final double e7 = 0.000247453;
  
      if(k < 1.0) { // Base case, for small k
        final double b = 1.0 + 0.36788794412 * k; // Step 1
        while(true) {
          final double p = b * random.nextDouble();
          if(p <= 1.0) { // when gds <= 1
            final double gds = FastMath.exp(FastMath.log(p) / k);
            if(FastMath.log(random.nextDouble()) <= -gds) {
              return (gds / theta);
            }
          }
          else { // when gds > 1
            final double gds = -FastMath.log((b - p) / k);
            if(FastMath.log(random.nextDouble()) <= ((k - 1.0) * FastMath.log(gds))) {
              return (gds / theta);
            }
          }
        }
      }
      else {
        // Step 1. Preparations
        final double ss, s, d;
        if(k != -1.0) {
          ss = k - 0.5;
          s = FastMath.sqrt(ss);
          d = 5.656854249 - 12.0 * s;
        }
        else {
          // For k == -1.0:
          ss = 0.0;
          s = 0.0;
          d = 0.0;
        }
        // Random vector of maximum length 1
        final double v1, /* v2, */v12;
        { // Temporary values - candidate
          double tv1, tv2, tv12;
          do {
            tv1 = 2.0 * random.nextDouble() - 1.0;
            tv2 = 2.0 * random.nextDouble() - 1.0;
            tv12 = tv1 * tv1 + tv2 * tv2;
          }
          while(tv12 > 1.0);
          v1 = tv1;
          /* v2 = tv2; */
          v12 = tv12;
        }
  
        // double b = 0.0, c = 0.0;
        // double si = 0.0, q0 = 0.0;
        final double b, c, si, q0;
  
        // Simpler accept cases & parameter computation
        {
          final double t = v1 * FastMath.sqrt(-2.0 * FastMath.log(v12) / v12);
          final double x = s + 0.5 * t;
          final double gds = x * x;
          if(t >= 0.0) {
            return (gds / theta); // Immediate acceptance
          }
  
          // Random uniform
          final double un = random.nextDouble();
          // Squeeze acceptance
          if(d * un <= t * t * t) {
            return (gds / theta);
          }
  
          if(k != -1.0) { // Step 4. Set-up for hat case
            final double r = 1.0 / k;
            q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r;
            if(k > 3.686) {
              if(k > 13.022) {
                b = 1.77;
                si = 0.75;
                c = 0.1515 / s;
              }
              else {
                b = 1.654 + 0.0076 * ss;
                si = 1.68 / s + 0.275;
                c = 0.062 / s + 0.024;
              }
            }
            else {
              b = 0.463 + s - 0.178 * ss;
              si = 1.235;
              c = 0.195 / s - 0.079 + 0.016 * s;
            }
          }
          else {
            // For k == -1.0:
            b = 0.0;
            c = 0.0;
            si = 0.0;
            q0 = 0.0;
          }
          // Compute v and q
          if(x > 0.0) {
            final double v = t / (s + s);
            final double q;
            if(Math.abs(v) > 0.25) {
              q = q0 - s * t + 0.25 * t * t + (ss + ss) * FastMath.log(1.0 + v);
            }
            else {
              q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
            }
            // Quotient acceptance:
            if(FastMath.log(1.0 - un) <= q) {
              return (gds / theta);
            }
          }
        }
  
        // Double exponential deviate t
        while(true) {
          double e, u, sign_u, t;
          // Retry until t is sufficiently large
          do {
            e = -FastMath.log(random.nextDouble());
            u = random.nextDouble();
            u = u + u - 1.0;
            sign_u = (u > 0) ? 1.0 : -1.0;
            t = b + (e * si) * sign_u;
          }
          while(t <= -0.71874483771719);
  
          // New v(t) and q(t)
          final double v = t / (s + s);
          final double q;
          if(Math.abs(v) > 0.25) {
            q = q0 - s * t + 0.25 * t * t + (ss + ss) * FastMath.log(1.0 + v);
          }
          else {
            q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
          }
          if(q <= 0.0) {
            continue; // retry
          }
          // Compute w(t)
          final double w;
          if(q > 0.5) {
            w = FastMath.exp(q) - 1.0;
          }
          else {
            w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * q + e1) * q;
          }
          // Hat acceptance
          if(c * u * sign_u <= w * FastMath.exp(e - 0.5 * t * t)) {
            final double x = s + 0.5 * t;
            return (x * x / theta);
          }
        }
      }
    }