question: 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 ; } }
target:  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.