double devroye_tail_gamma_ran (gsl_rng *rng, double tail, double a)
{
  double U, E, X, Es; 

  if (a <= 0) {
    
    X = 0;
    
  }
  else if (a < 1) {

    while (1) {
      U = gsl_ran_flat (rng, 0., 1.);
      E = gsl_ran_exponential (rng, 1.);
      X = tail + E;
      if (X * pow (U, 1/(1-a)) <= a) break;
    }

  }
  else if (a == 1) {

    E = gsl_ran_exponential (rng, 1.);
    X = tail + E;

  }
  else {

    if (tail > a-1) { // mode left 

      while (1) {
	E = gsl_ran_exponential (rng, 1.);
	Es = gsl_ran_exponential (rng, 1.);
	X = tail + E / (1 - (a-1)/tail);
	if (X/tail - 1 + log (tail/X) <= Es/(a-1)) break;
      }

    }
    else {

      while (1) { // mode right
	X = gsl_ran_gamma (rng, a, 1);
	if (X >= tail) break;

      }
    }
  }
  return X;
}
