define sample_gamma(alpha){
  /*
  A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001
  https://dl.acm.org/doi/pdf/10.1145/358407.358414
  see also the references/ folder
  Note that the Wikipedia page for the gamma distribution includes a scaling parameter
  k or beta
  https://en.wikipedia.org/wiki/Gamma_distribution
  such that gamma_k(alpha, k) = k * gamma(alpha)
  or gamma_beta(alpha, beta) = gamma(alpha) / beta
  So far I have not needed to use this, and thus the second parameter is by default 1.
  */

  if (alpha >= 1) {
    d = alpha - (1/3);
    c = 1.0 / sqrt(9.0 * d);
    while (1) {
      v=-1
      while(v<=0) {
        x = sample_unit_normal();
        v = 1 + c * x;
      }
      v = v * v * v;
      u = sample_unit_uniform();
      if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */
        /* 
        the 0.0331 doesn't inspire much confidence
        however, this isn't the whole story
        by knowing that Condition 1 implies condition 2
        we realize that this is just a way of making the algorithm faster
        i.e., of not using the logarithms
        */
        return d * v;
      }
      if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */
        return d * v;
      }
    }
  } else {
    return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha);
    /* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */
  }
}

define sample_beta(a, b)
{
    /* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */
    gamma_a = sample_gamma(a);
    gamma_b = sample_gamma(b);
    return gamma_a / (gamma_a + gamma_b);
}