protected int hmdu(int N,int M,int n,RandomEngine randomGenerator){
  int I, K;
  double p, nu, c, d, U;
  if (N != N_last || M != M_last || n != n_last) {
    N_last=N;
    M_last=M;
    n_last=n;
    Mp=(double)(M + 1);
    np=(double)(n + 1);
    N_Mn=N - M - n;
    p=Mp / (N + 2.0);
    nu=np * p;
    if ((m=(int)nu) == nu && p == 0.5) {
      mp=m--;
    }
 else {
      mp=m + 1;
    }
    fm=Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m) + Arithmetic.logFactorial(M) - Arithmetic.logFactorial(M - m) - Arithmetic.logFactorial(m) - Arithmetic.logFactorial(N) + Arithmetic.logFactorial(N - n) + Arithmetic.logFactorial(n));
    b=(int)(nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n / (double)N) + 1.0));
    if (b > n)     b=n;
  }
  for (; ; ) {
    if ((U=randomGenerator.raw() - fm) <= 0.0)     return (m);
    c=d=fm;
    for (I=1; I <= m; I++) {
      K=mp - I;
      c*=(double)K / (np - K) * ((double)(N_Mn + K) / (Mp - K));
      if ((U-=c) <= 0.0)       return (K - 1);
      K=m + I;
      d*=(np - K) / (double)K * ((Mp - K) / (double)(N_Mn + K));
      if ((U-=d) <= 0.0)       return (K);
    }
    for (K=mp + m; K <= b; K++) {
      d*=(np - K) / (double)K * ((Mp - K) / (double)(N_Mn + K));
      if ((U-=d) <= 0.0)       return (K);
    }
  }
}

protected int hmdu ( int N , int M , int n , RandomEngine randomGenerator ) { int I , K ; double p , nu , c , d , U ; if ( N != N_last || M != M_last || n != n_last ) { N_last = N ; M_last = M ; n_last = n ; Mp = ( double ) ( M + 1 ) ; np = ( double ) ( n + 1 ) ; N_Mn = N - M - n ; p = Mp / ( N + 2.0 ) ; nu = np * p ; if ( ( m = ( int ) nu ) == nu && p == 0.5 ) { mp = m -- ; } else { mp = m + 1 ; } fm = Math . exp ( Arithmetic . logFactorial ( N - M ) - Arithmetic . logFactorial ( N_Mn + m ) - Arithmetic . logFactorial ( n - m ) + Arithmetic . logFactorial ( M ) - Arithmetic . logFactorial ( M - m ) - Arithmetic . logFactorial ( m ) - Arithmetic . logFactorial ( N ) + Arithmetic . logFactorial ( N - n ) + Arithmetic . logFactorial ( n ) ) ; b = ( int ) ( nu + 11.0 * Math . sqrt ( nu * ( 1.0 - p ) * ( 1.0 - n / ( double ) N ) + 1.0 ) ) ; if ( b > n ) b = n ; } for ( ; ; ) { if ( ( U = randomGenerator . raw ( ) - fm ) <= 0.0 ) return ( m ) ; c = d = fm ; for ( I = 1 ; I <= m ; I ++ ) { K = mp - I ; c *= ( double ) K / ( np - K ) * ( ( double ) ( N_Mn + K ) / ( Mp - K ) ) ; if ( ( U -= c ) <= 0.0 ) return ( K - 1 ) ; K = m + I ; d *= ( np - K ) / ( double ) K * ( ( Mp - K ) / ( double ) ( N_Mn + K ) ) ; if ( ( U -= d ) <= 0.0 ) return ( K ) ; } for ( K = mp + m ; K <= b ; K ++ ) { d *= ( np - K ) / ( double ) K * ( ( Mp - K ) / ( double ) ( N_Mn + K ) ) ; if ( ( U -= d ) <= 0.0 ) return ( K ) ; } } }
returns a random number from the distribution .
