protected int generateBinomial(int n,double p){
  final double C1_3=0.33333333333333333;
  final double C5_8=0.62500000000000000;
  final double C1_6=0.16666666666666667;
  final int DMAX_KM=20;
  int bh, i, K, Km, nK;
  double f, rm, U, V, X, T, E;
  if (n != n_last || p != p_last) {
    n_last=n;
    p_last=p;
    par=Math.min(p,1.0 - p);
    q=1.0 - par;
    np=n * par;
    if (np <= 0.0)     return -1;
    rm=np + par;
    m=(int)rm;
    if (np < 10) {
      p0=Math.exp(n * Math.log(q));
      bh=(int)(np + 10.0 * Math.sqrt(np * q));
      b=Math.min(n,bh);
    }
 else {
      rc=(n + 1.0) * (pq=par / q);
      ss=np * q;
      i=(int)(2.195 * Math.sqrt(ss) - 4.6 * q);
      xm=m + 0.5;
      xl=(double)(m - i);
      xr=(double)(m + i + 1L);
      f=(rm - xl) / (rm - xl * par);
      ll=f * (1.0 + 0.5 * f);
      f=(xr - rm) / (xr * q);
      lr=f * (1.0 + 0.5 * f);
      c=0.134 + 20.5 / (15.3 + (double)m);
      p1=i + 0.5;
      p2=p1 * (1.0 + c + c);
      p3=p2 + c / ll;
      p4=p3 + c / lr;
    }
  }
  if (np < 10) {
    double pk;
    K=0;
    pk=p0;
    U=randomGenerator.raw();
    while (U > pk) {
      ++K;
      if (K > b) {
        U=randomGenerator.raw();
        K=0;
        pk=p0;
      }
 else {
        U-=pk;
        pk=(double)(((n - K + 1) * par * pk) / (K * q));
      }
    }
    return ((p > 0.5) ? (n - K) : K);
  }
  for (; ; ) {
    V=randomGenerator.raw();
    if ((U=randomGenerator.raw() * p4) <= p1) {
      K=(int)(xm - U + p1 * V);
      return (p > 0.5) ? (n - K) : K;
    }
    if (U <= p2) {
      X=xl + (U - p1) / c;
      if ((V=V * c + 1.0 - Math.abs(xm - X) / p1) >= 1.0)       continue;
      K=(int)X;
    }
 else     if (U <= p3) {
      if ((X=xl + Math.log(V) / ll) < 0.0)       continue;
      K=(int)X;
      V*=(U - p2) * ll;
    }
 else {
      if ((K=(int)(xr - Math.log(V) / lr)) > n)       continue;
      V*=(U - p3) * lr;
    }
    if ((Km=Math.abs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss) {
      f=1.0;
      if (m < K) {
        for (i=m; i < K; ) {
          if ((f*=(rc / ++i - pq)) < V)           break;
        }
      }
 else {
        for (i=K; i < m; ) {
          if ((V*=(rc / ++i - pq)) > f)           break;
        }
      }
      if (V <= f)       break;
    }
 else {
      V=Math.log(V);
      T=-Km * Km / (ss + ss);
      E=(Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
      if (V <= T - E)       break;
      if (V <= T + E) {
        if (n != n_prev || par != p_prev) {
          n_prev=n;
          p_prev=par;
          nm=n - m + 1;
          ch=xm * Math.log((m + 1.0) / (pq * nm)) + Arithmetic.stirlingCorrection(m + 1) + Arithmetic.stirlingCorrection(nm);
        }
        nK=n - K + 1;
        if (V <= ch + (n + 1.0) * Math.log((double)nm / (double)nK) + (K + 0.5) * Math.log(nK * pq / (K + 1.0)) - Arithmetic.stirlingCorrection(K + 1) - Arithmetic.stirlingCorrection(nK))         break;
      }
    }
  }
  return (p > 0.5) ? (n - K) : K;
}

protected int generateBinomial ( int n , double p ) { final double C1_3 = 0.33333333333333333 ; final double C5_8 = 0.62500000000000000 ; final double C1_6 = 0.16666666666666667 ; final int DMAX_KM = 20 ; int bh , i , K , Km , nK ; double f , rm , U , V , X , T , E ; if ( n != n_last || p != p_last ) { n_last = n ; p_last = p ; par = Math . min ( p , 1.0 - p ) ; q = 1.0 - par ; np = n * par ; if ( np <= 0.0 ) return - 1 ; rm = np + par ; m = ( int ) rm ; if ( np < 10 ) { p0 = Math . exp ( n * Math . log ( q ) ) ; bh = ( int ) ( np + 10.0 * Math . sqrt ( np * q ) ) ; b = Math . min ( n , bh ) ; } else { rc = ( n + 1.0 ) * ( pq = par / q ) ; ss = np * q ; i = ( int ) ( 2.195 * Math . sqrt ( ss ) - 4.6 * q ) ; xm = m + 0.5 ; xl = ( double ) ( m - i ) ; xr = ( double ) ( m + i + 1L ) ; f = ( rm - xl ) / ( rm - xl * par ) ; ll = f * ( 1.0 + 0.5 * f ) ; f = ( xr - rm ) / ( xr * q ) ; lr = f * ( 1.0 + 0.5 * f ) ; c = 0.134 + 20.5 / ( 15.3 + ( double ) m ) ; p1 = i + 0.5 ; p2 = p1 * ( 1.0 + c + c ) ; p3 = p2 + c / ll ; p4 = p3 + c / lr ; } } if ( np < 10 ) { double pk ; K = 0 ; pk = p0 ; U = randomGenerator . raw ( ) ; while ( U > pk ) { ++ K ; if ( K > b ) { U = randomGenerator . raw ( ) ; K = 0 ; pk = p0 ; } else { U -= pk ; pk = ( double ) ( ( ( n - K + 1 ) * par * pk ) / ( K * q ) ) ; } } return ( ( p > 0.5 ) ? ( n - K ) : K ) ; } for ( ; ; ) { V = randomGenerator . raw ( ) ; if ( ( U = randomGenerator . raw ( ) * p4 ) <= p1 ) { K = ( int ) ( xm - U + p1 * V ) ; return ( p > 0.5 ) ? ( n - K ) : K ; } if ( U <= p2 ) { X = xl + ( U - p1 ) / c ; if ( ( V = V * c + 1.0 - Math . abs ( xm - X ) / p1 ) >= 1.0 ) continue ; K = ( int ) X ; } else if ( U <= p3 ) { if ( ( X = xl + Math . log ( V ) / ll ) < 0.0 ) continue ; K = ( int ) X ; V *= ( U - p2 ) * ll ; } else { if ( ( K = ( int ) ( xr - Math . log ( V ) / lr ) ) > n ) continue ; V *= ( U - p3 ) * lr ; } if ( ( Km = Math . abs ( K - m ) ) <= DMAX_KM || Km + Km + 2L >= ss ) { f = 1.0 ; if ( m < K ) { for ( i = m ; i < K ; ) { if ( ( f *= ( rc / ++ i - pq ) ) < V ) break ; } } else { for ( i = K ; i < m ; ) { if ( ( V *= ( rc / ++ i - pq ) ) > f ) break ; } } if ( V <= f ) break ; } else { V = Math . log ( V ) ; T = - Km * Km / ( ss + ss ) ; E = ( Km / ss ) * ( ( Km * ( Km * C1_3 + C5_8 ) + C1_6 ) / ss + 0.5 ) ; if ( V <= T - E ) break ; if ( V <= T + E ) { if ( n != n_prev || par != p_prev ) { n_prev = n ; p_prev = par ; nm = n - m + 1 ; ch = xm * Math . log ( ( m + 1.0 ) / ( pq * nm ) ) + Arithmetic . stirlingCorrection ( m + 1 ) + Arithmetic . stirlingCorrection ( nm ) ; } nK = n - K + 1 ; if ( V <= ch + ( n + 1.0 ) * Math . log ( ( double ) nm / ( double ) nK ) + ( K + 0.5 ) * Math . log ( nK * pq / ( K + 1.0 ) ) - Arithmetic . stirlingCorrection ( K + 1 ) - Arithmetic . stirlingCorrection ( nK ) ) break ; } } } return ( p > 0.5 ) ? ( n - K ) : K ; }
* binomial - distribution - acceptance rejection / inversion * * acceptance rejection method combined with inversion for generating binomial random numbers with parameters n ( number of trials ) and p ( probability of success ) . for min ( n * p , n * ( 1 - p ) ) < 10 the inversion method is applied : the random numbers are generated via sequential search , starting at the lowest index k = 0 . the cumulative probabilities are avoided by using the technique of chop - down . for min ( n * p , n * ( 1 - p ) ) > = 10 acceptance rejection is used : the algorithm is based on a hat - function which is uniform in the centre region and exponential in the tails . a triangular immediate acceptance region in the centre speeds up the generation of binomial variates . if candidate k is near the mode , f ( k ) is computed recursively starting at the mode m . the acceptance test by stirling ' s formula is modified according to w . hoermann ( 1992 ) : the generation of binomial random variates , to appear in j . statist . comput . simul . if p < . 5 the algorithm is applied to parameters n , p . otherwise p is replaced by 1 - p , and k is replaced by n - k . * * function : - samples a random number from the binomial distribution with parameters n and p and is valid for n * min ( p , 1 - p ) > 0 . reference : - v . kachitvichyanukul , b . w . schmeiser ( 1988 ) : binomial random variate generation , communications of the acm 31 , 216 - 222 . subprograms : - stirlingcorrection ( ) . . . correction term of the stirling approximation for log ( k ! ) ( series in 1 / k or table values for small k ) with long int k - randomgenerator . . . ( 0 , 1 ) - uniform engine * *
