public double nextGamma(double a) throws Exception {
  if (a <= 0.0)   throw new Exception("Shape parameter of Gamma distribution" + "must be greater than 0!");
 else   if (a == 1.0)   return nextExponential();
 else   if (a < 1.0) {
    double b=1.0 + Math.exp(-1.0) * a, p, x, condition;
    do {
      p=b * super.nextDouble();
      if (p < 1.0) {
        x=Math.exp(Math.log(p) / a);
        condition=x;
      }
 else {
        x=-Math.log((b - p) / a);
        condition=(1.0 - a) * Math.log(x);
      }
    }
 while (nextExponential() < condition);
    return x;
  }
 else {
    double b=a - 1.0, D=Math.sqrt(b), D1, x1, x2, xl, f1, f2, x4, x5, xr, f4, f5, p1, p2, p3, p4;
    if (a <= 2.0) {
      D1=b / 2.0;
      x1=0.0;
      x2=D1;
      xl=-1.0;
      f1=0.0;
    }
 else {
      D1=D - 0.5;
      x2=b - D1;
      x1=x2 - D1;
      xl=1.0 - b / x1;
      f1=Math.exp(b * Math.log(x1 / b) + 2.0 * D1);
    }
    f2=Math.exp(b * Math.log(x2 / b) + D1);
    x4=b + D;
    x5=x4 + D;
    xr=1.0 - b / x5;
    f4=Math.exp(b * Math.log(x4 / b) - D);
    f5=Math.exp(b * Math.log(x5 / b) - 2.0 * D);
    p1=2.0 * f4 * D;
    p2=2.0 * f2 * D1 + p1;
    p3=f5 / xr + p2;
    p4=-f1 / xl + p3;
    double u, w=Double.MAX_VALUE, x=b, v, xp;
    while (Math.log(w) > (b * Math.log(x / b) + b - x) || x < 0.0) {
      u=super.nextDouble() * p4;
      if (u <= p1) {
        w=u / D - f4;
        if (w <= 0.0)         return (b + u / f4);
        if (w <= f5)         return (x4 + (w * D) / f5);
        v=super.nextDouble();
        x=x4 + v * D;
        xp=2.0 * x4 - x;
        if (w >= f4 + (f4 - 1) * (x - x4) / (x4 - b))         return xp;
        if (w <= f4 + (b / x4 - 1) * f4 * (x - x4))         return x;
        if ((w < 2.0 * f4 - 1.0) || (w < 2.0 * f4 - Math.exp(b * Math.log(xp / b) + b - xp)))         continue;
        return xp;
      }
 else       if (u <= p2) {
        w=(u - p1) / D1 - f2;
        if (w <= 0.0)         return (b - (u - p1) / f2);
        if (w <= f1)         return (x1 + w * D1 / f1);
        v=super.nextDouble();
        x=x1 + v * D1;
        xp=2.0 * x2 - x;
        if (w >= f2 + (f2 - 1) * (x - x2) / (x2 - b))         return xp;
        if (w <= f2 * (x - x1) / D1)         return x;
        if ((w < 2.0 * f2 - 1.0) || (w < 2.0 * f2 - Math.exp(b * Math.log(xp / b) + b - xp)))         continue;
        return xp;
      }
 else       if (u < p3) {
        w=super.nextDouble();
        u=(p3 - u) / (p3 - p2);
        x=x5 - Math.log(u) / xr;
        if (w <= (xr * (x5 - x) + 1.0) / u)         return x;
        w=w * f5 * u;
      }
 else {
        w=super.nextDouble();
        u=(p4 - u) / (p4 - p3);
        x=x1 - Math.log(u) / xl;
        if (x < 0.0)         continue;
        if (w <= (xl * (x1 - x) + 1.0) / u)         return x;
        w=w * f1 * u;
      }
    }
    return x;
  }
}

public double nextGamma ( double a ) throws Exception { if ( a <= 0.0 ) throw new Exception ( "Shape parameter of Gamma distribution" + "must be greater than 0!" ) ; else if ( a == 1.0 ) return nextExponential ( ) ; else if ( a < 1.0 ) { double b = 1.0 + Math . exp ( - 1.0 ) * a , p , x , condition ; do { p = b * super . nextDouble ( ) ; if ( p < 1.0 ) { x = Math . exp ( Math . log ( p ) / a ) ; condition = x ; } else { x = - Math . log ( ( b - p ) / a ) ; condition = ( 1.0 - a ) * Math . log ( x ) ; } } while ( nextExponential ( ) < condition ) ; return x ; } else { double b = a - 1.0 , D = Math . sqrt ( b ) , D1 , x1 , x2 , xl , f1 , f2 , x4 , x5 , xr , f4 , f5 , p1 , p2 , p3 , p4 ; if ( a <= 2.0 ) { D1 = b / 2.0 ; x1 = 0.0 ; x2 = D1 ; xl = - 1.0 ; f1 = 0.0 ; } else { D1 = D - 0.5 ; x2 = b - D1 ; x1 = x2 - D1 ; xl = 1.0 - b / x1 ; f1 = Math . exp ( b * Math . log ( x1 / b ) + 2.0 * D1 ) ; } f2 = Math . exp ( b * Math . log ( x2 / b ) + D1 ) ; x4 = b + D ; x5 = x4 + D ; xr = 1.0 - b / x5 ; f4 = Math . exp ( b * Math . log ( x4 / b ) - D ) ; f5 = Math . exp ( b * Math . log ( x5 / b ) - 2.0 * D ) ; p1 = 2.0 * f4 * D ; p2 = 2.0 * f2 * D1 + p1 ; p3 = f5 / xr + p2 ; p4 = - f1 / xl + p3 ; double u , w = Double . MAX_VALUE , x = b , v , xp ; while ( Math . log ( w ) > ( b * Math . log ( x / b ) + b - x ) || x < 0.0 ) { u = super . nextDouble ( ) * p4 ; if ( u <= p1 ) { w = u / D - f4 ; if ( w <= 0.0 ) return ( b + u / f4 ) ; if ( w <= f5 ) return ( x4 + ( w * D ) / f5 ) ; v = super . nextDouble ( ) ; x = x4 + v * D ; xp = 2.0 * x4 - x ; if ( w >= f4 + ( f4 - 1 ) * ( x - x4 ) / ( x4 - b ) ) return xp ; if ( w <= f4 + ( b / x4 - 1 ) * f4 * ( x - x4 ) ) return x ; if ( ( w < 2.0 * f4 - 1.0 ) || ( w < 2.0 * f4 - Math . exp ( b * Math . log ( xp / b ) + b - xp ) ) ) continue ; return xp ; } else if ( u <= p2 ) { w = ( u - p1 ) / D1 - f2 ; if ( w <= 0.0 ) return ( b - ( u - p1 ) / f2 ) ; if ( w <= f1 ) return ( x1 + w * D1 / f1 ) ; v = super . nextDouble ( ) ; x = x1 + v * D1 ; xp = 2.0 * x2 - x ; if ( w >= f2 + ( f2 - 1 ) * ( x - x2 ) / ( x2 - b ) ) return xp ; if ( w <= f2 * ( x - x1 ) / D1 ) return x ; if ( ( w < 2.0 * f2 - 1.0 ) || ( w < 2.0 * f2 - Math . exp ( b * Math . log ( xp / b ) + b - xp ) ) ) continue ; return xp ; } else if ( u < p3 ) { w = super . nextDouble ( ) ; u = ( p3 - u ) / ( p3 - p2 ) ; x = x5 - Math . log ( u ) / xr ; if ( w <= ( xr * ( x5 - x ) + 1.0 ) / u ) return x ; w = w * f5 * u ; } else { w = super . nextDouble ( ) ; u = ( p4 - u ) / ( p4 - p3 ) ; x = x1 - Math . log ( u ) / xl ; if ( x < 0.0 ) continue ; if ( w <= ( xl * ( x1 - x ) + 1.0 ) / u ) return x ; w = w * f1 * u ; } } return x ; } }
generate a value of a variate following standard gamma distribution with shape parameter a . < p > if a > 1 , it uses a rejection method developed by minh ( 1988 ) " generating gamma variates " , acm trans . on math . software , vol . 14 , no . 3 , pp261 - 266 . < br > if a < 1 , it uses the algorithm " gs " developed by ahrens and dieter ( 1974 ) " computer methods for sampling from gamma , beta , poisson and binomial distributions " , computing , 12 ( 1974 ) , pp223 - 246 , and further implemented in fortran by ahrens , kohrt and dieter ( 1983 ) " algorithm 599 : sampling from gamma and poisson distributions " , acm trans . on math . software , vol . 9 no . 2 , pp255 - 257 . < p > variates related to standard gamma can be generated using simple transformations .
