Function random_beta

numpy/random/src/distributions/distributions.c:411–457

* Note: random_beta assumes that a != 0 and b != 0. */

 *  Note: random_beta assumes that a != 0 and b != 0.
 */
double random_beta(bitgen_t *bitgen_state, double a, double b) {
  double Ga, Gb;

  if ((a <= 1.0) && (b <= 1.0)) {
    double U, V, X, Y, XpY;

    if (a < BETA_TINY_THRESHOLD && b < BETA_TINY_THRESHOLD) {
      /*
       * When a and b are this small, the probability that the
       * sample would be a double precision float that is not
       * 0 or 1 is less than approx. 1e-100.  So we use the
       * proportion a/(a + b) and a single uniform sample to
       * generate the result.
       */
      U = next_double(bitgen_state);
      return (a + b)*U < a;
    }

    /* Use Johnk's algorithm */

    while (1) {
      U = next_double(bitgen_state);
      V = next_double(bitgen_state);
      X = pow(U, 1.0 / a);
      Y = pow(V, 1.0 / b);
      XpY = X + Y;
      /* Reject if both U and V are 0.0, which is approx 1 in 10^106 */
      if ((XpY <= 1.0) && (U + V > 0.0)) {
        if (XpY > 0) {
          return X / XpY;
        } else {
          double logX = log(U) / a;
          double logY = log(V) / b;
          double logM = logX > logY ? logX : logY;
          logX -= logM;
          logY -= logM;

          return exp(logX - log(exp(logX) + exp(logY)));
        }
      }
    }
  } else {
    Ga = random_standard_gamma(bitgen_state, a);
    Gb = random_standard_gamma(bitgen_state, b);
    return Ga / (Ga + Gb);
  }
}

double random_chisquare(bitgen_t *bitgen_state, double df) {
  return 2.0 * random_standard_gamma(bitgen_state, df / 2.0);