question: public static int binomialRand ( int n , double pp ) { int j , k ; double am , em , g , p , sq , t , y ; double pc , plog , pclog , en ; p = ( pp <= 0.5 ) ? pp : 1.0 - pp ; am = n * p ; if ( p == 0.0 ) { k = 0 ; } else if ( p == 1.0 ) { k = n ; } else if ( n < 50 ) { k = 0 ; for ( j = 0 ; j < n ; j ++ ) { if ( uniformRand ( ) < p ) { k ++ ; } } } else if ( am < 1.0 ) { g = Math . exp ( - am ) ; t = 1.0 ; k = - 1 ; do { k ++ ; t *= uniformRand ( ) ; } while ( t > g ) ; if ( k > n ) { k = n ; } } else { en = n ; g = lngamma ( en + 1.0 ) ; pc = 1.0 - p ; plog = Math . log ( p ) ; pclog = Math . log ( pc ) ; sq = Math . sqrt ( 2.0 * am * pc ) ; do { do { y = Math . tan ( Math . PI * uniformRand ( ) ) ; em = sq * y + am ; } while ( ( em < 0.0 ) || ( em >= en + 1.0 ) ) ; em = Math . floor ( em ) ; t = 1.2 * sq * ( 1.0 + y * y ) * Math . exp ( g - lngamma ( em + 1.0 ) - lngamma ( en - em + 1.0 ) + em * plog + ( en - em ) * pclog ) ; } while ( uniformRand ( ) > t ) ; k = ( int ) em ; } if ( p != pp ) { k = n - k ; } return ( k ) ; }
target:  Binomial random generator from Numerical Recipes