! This file now contains TOMS algorithm 668 as a separate module.
! Latest revision - 24 July 2000
! amiller@bigpond.net.au
! http://users.bigpond.net.au/amiller/


MODULE random_hypergeometric
! Generate random deviates from the hypergeometric distribution.
! The notation used for the 2 x 2 table is:

!   x  r-x |  r
!          |
!  g-x     | N-r
! _________|_____
!          |
!   g  N-g |  N

! where N, r, g are fixed, and x is the random value to be generated.

! Two methods are included.

! r_hyperg is a simple method which is suitable when MIN(r,g) is small,
!          say <= 10.

! h_alias is much faster for larger values of MIN(r,g) but requires a table
!         of hypergeometric probabilities to be set up on the first call.
!         It is not suitable if the marginal totals (N, r, g) change.

! Latest revision - 24 July 2000

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), ALLOCATABLE, SAVE  :: p(:), f(:)
INTEGER, ALLOCATABLE, SAVE    :: l(:)
INTEGER, SAVE                 :: Nlast = 0, rlast = 0, glast = 0, xmax


CONTAINS


FUNCTION r_hyperg(N, r, g) RESULT(x)

! This is algorithm 3.1 from page 52 of:
!   Dagpunar, John (1988) `Principles of Random Variate Generation'
!   published by Clarendon Press, Oxford.   ISBN 0-19-852202-9

INTEGER, INTENT(IN)  :: N, r, g
INTEGER              :: x

! Local variables

INTEGER  :: gg, i, NN
REAL     :: p, u

gg = g
NN = N
i = 1
x = 0
DO
  IF (gg <= 0 .OR. i > r) EXIT
  p = REAL(gg) / NN
  CALL RANDOM_NUMBER(u)
  IF (u < p) THEN
    x = x + 1
    gg = gg - 1
  END IF
  i = i + 1
  NN = NN - 1
END DO

RETURN
END FUNCTION r_hyperg



FUNCTION h_alias(N, r, g) RESULT(x)

! This is algorithm HALIAS from:

! Kachitvichyanukul, V. & Schmeiser, B. (1985)
! `Computer Generation of Hypergeometric Random Variates',
! J. of Statistical Computation & Simulation, vol.22, 127-145.

INTEGER, INTENT(IN)  :: N, r, g
INTEGER              :: x

! Local variables
INTEGER    :: i, i0, nn(2,2)
REAL (dp)  :: u, p_last

! If this is the first call, or N, r or g has changed, call alias
! to set up new tables.

IF (N /= Nlast .OR. r /= rlast .OR. g /= glast) THEN
  Nlast = N
  rlast = r
  glast = g
  xmax = MIN(r,g)
  IF (ALLOCATED( p )) DEALLOCATE( p, f, l )
  ALLOCATE( p(0:xmax), f(0:xmax), l(0:xmax) )

! Calculate the set of hypergeometric probabilities
! First calculate the maximum probability which occurs near
! i = (r+1)(g+1)/(N+2)

  i0 = NINT( REAL((r+1)*(g+1)) / REAL(N+2) )
  i0 = MIN(i0, xmax)
  nn(1,1) = i0
  nn(1,2) = r - i0
  nn(2,1) = g - i0
  nn(2,2) = N - r - g + i0
  p(i0) = hyperg(nn)

  p_last = p(i0)
  DO i = i0-1, 0, -1
    p(i) = p_last * REAL( (i+1)*(N-r-g+i+1), KIND=dp ) /  &
                    REAL( (r-i)*(g-i), KIND=dp )
    p_last = p(i)
  END DO

  p_last = p(i0)
  DO i = i0+1, xmax
    p(i) = p_last * REAL( (r-i+1)*(g-i+1), KIND=dp ) /  &
                    REAL( i*(N-r-g+i), KIND=dp )
    p_last = p(i)
  END DO

! Call ALIAS to set up its tables of F and L.

  CALL alias(p, xmax+1, f, l)
END IF

! Generate an integer in the range 0 to XMAX
CALL RANDOM_NUMBER(u)
x = u * (xmax + 1)

! Accept it if it is <= f(x), otherwise use its alias.
CALL RANDOM_NUMBER(u)
IF (u > f(x)) x = l(x) - 1

RETURN
END FUNCTION h_alias



FUNCTION hyperg(n) RESULT(prob)

!  Calculate hypergeometric probability for a 2 x 2 table:
!      n(1,1)  n(1,2)
!      n(2,1)  n(2,2)

!  For large numbers, it would be faster to use Stirling's
!  approximation or similar for the factorials.

!  Programmer: Alan Miller
!  Latest revision - 9 January 1997

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n(2,2)
REAL (dp)            :: prob

!     Local variables

INTEGER    :: ncell(1:4), ntot(1:4), nall
LOGICAL    :: ended
REAL (dp)  :: zero = 0.0_dp, one = 1.0_dp

!     Check that all elements are positive or zero.

IF (n(1,1) < 0 .OR. n(1,2) < 0 .OR. n(2,1) < 0 .OR. n(2,2) < 0) THEN
  WRITE(*, *) '** Negative cell frequency in argument to HYPERG **'
  prob = zero
  RETURN
END IF

!     Copy the cell frequencies into ncell and order them.
!     Calculate the row & column totals and order them.

ncell(1) = n(1,1)
ncell(2) = n(1,2)
ncell(3) = n(2,1)
ncell(4) = n(2,2)
CALL order4(ncell)

ntot(1) = n(1,1) + n(1,2)
ntot(2) = n(2,1) + n(2,2)
ntot(3) = n(1,1) + n(2,1)
ntot(4) = n(1,2) + n(2,2)
CALL order4(ntot)

nall = ncell(1) + ncell(2) + ncell(3) + ncell(4)
IF (nall == 0) THEN
  WRITE(*, *) '** Call to HYPERG with total cell frequency = 0 **'
  prob = zero
  RETURN
END IF

!     Calculate hyperg as:

!       ntot(1)!   ntot(2)!    ntot(3)!    ntot(4)!        1
!       -------- . --------- . --------- . --------- . ---------
!        nall!     ncell(1)!   ncell(2)!   ncell(3)!   ncell(4)!

!     The order of multiplications and divisions minimizes the chance
!     of overflows or underflows without using logarithms.

prob = one
DO
  ended = .true.
  IF (nall > ntot(1)) THEN
    prob = prob / nall
    nall = nall - 1
    ended = .false.
  END IF

  IF (ntot(2) > ncell(1)) THEN
    prob = prob * ntot(2)
    ntot(2) = ntot(2) - 1
    ended = .false.
  END IF

  IF (ntot(3) > ncell(2)) THEN
    prob = prob * ntot(3)
    ntot(3) = ntot(3) - 1
    ended = .false.
  END IF

  IF (ntot(4) > ncell(3)) THEN
    prob = prob * ntot(4)
    ntot(4) = ntot(4) - 1
    ended = .false.
  END IF

  IF (ncell(4) > 1) THEN
    prob = prob / ncell(4)
    ncell(4) = ncell(4) - 1
    ended = .false.
  END IF
  IF (ended) EXIT
END DO

RETURN
END FUNCTION hyperg



SUBROUTINE order4(n)

!     Order 4 integers; largest first.

IMPLICIT NONE
INTEGER, INTENT(IN OUT) :: n(4)

!     Local variable
INTEGER :: nn

IF (n(1) < n(2)) THEN
  nn = n(1)
  n(1) = n(2)
  n(2) = nn
END IF

IF (n(3) < n(4)) THEN
  nn = n(3)
  n(3) = n(4)
  n(4) = nn
END IF

IF (n(2) < n(3)) THEN
  nn = n(2)
  n(2) = n(3)
  n(3) = nn
  IF (n(3) < n(4)) THEN
    nn = n(3)
    n(3) = n(4)
    n(4) = nn
  END IF
  IF (n(1) < n(2)) THEN
    nn = n(1)
    n(1) = n(2)
    n(2) = nn
    IF (n(2) < n(3)) THEN
      nn = n(2)
      n(2) = n(3)
      n(3) = nn
    END IF
  END IF
END IF

RETURN
END SUBROUTINE order4



SUBROUTINE alias(p, n, f, l)

! Set up the aliases for the method described in:

! Kronmal, R.A. & Peterson, A.V. (1979)
! `On the Alias Method for Generating Random Variables from a Discrete
! Distribution', The American Statistician, vol.33(4), 214-218.

! The algorithm improves slightly on the basic algorithm given in:

! Walker, A.J. (1977)
! `An Efficient Method for Generating Discrete Random Variables with General
! Distributions', ACM Transactions on Mathematical Software, vol.3, 253-256.

! Arguments:
! p(1:n)  INPUT   Contains the discrete distribution
!                 N.B. It must have a finite range
! n       INPUT
! f(1:n)  OUTPUT  The cut-off values for each value of the random variable
! l(1:n)  OUTPUT  The aliases (at least one will be undefined, in which case
!                 it is not needed as the corresponding f(i) >= 1.

REAL (dp), INTENT(IN)   :: p(:)
INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(OUT)  :: f(:)
INTEGER, INTENT(OUT)    :: l(:)

! Local variables

INTEGER    :: i, ig, is, gptr(n), j, k, sptr(n-1)
REAL (dp), PARAMETER  :: nearly_one = 1.0_dp - 100.0_dp*EPSILON(1.0_dp)

! Step 1
f(1:n) = n * p(1:n)

! Step 2
! G = set of i such that f(i) >= 1
! S = the others
ig = 0
is = 0
DO i = 1, n
  IF (f(i) < 1.0) THEN
    is = is + 1
    sptr(is) = i
  ELSE
    ig = ig + 1
    gptr(ig) = i
  END IF
END DO

! Step 3
! Main loop
DO
  IF (is == 0) EXIT

! Step 4
! Choose element k in G, and element j in S.
! We choose the last ones on each list.
  k = gptr(ig)
  j = sptr(is)
! This finalizes operations on element j, so we can remove it from S.
  is = is - 1

! Step 5
  l(j) = k

! Step 6
  f(k) = f(k) - (1.0 - f(j))

! Step 7
! f(k) may now be < 1.
! If so, remove it from G and add it to S.
  IF (f(k) < nearly_one) THEN
    ig = ig - 1
    is = is + 1
    sptr(is) = k
  END IF
END DO

RETURN
END SUBROUTINE alias

END MODULE random_hypergeometric



MODULE toms668
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

INTEGER, SAVE       :: iseed = 123456789

CONTAINS


SUBROUTINE h2pec(kk, nn1, nn2, iseed1, jx)

!  HYPERGEOMETRIC RANDOM VARIATE GENERATOR

!  METHOD
!     IF (MODE - MAX(0,KK-NN2) .LT. 10), USE THE INVERSE CDF.
!        OTHERWISE, USE ALGORITHM H2PE: ACCEPTANCE-REJECTION VIA
!        THREE REGION COMPOSITION.  THE THREE REGIONS ARE A
!        RECTANGLE, AND EXPONENTIAL LEFT AND RIGHT TAILS.
!     H2PE  REFERS TO HYPERGEOMETRIC-2 POINTS-EXPONENTIAL TAILS.
!     H2PEC REFERS TO H2PE AND "COMBINED."  THUS H2PE IS THE RESEARCH RESULT
!        AND H2PEC IS THE IMPLEMENTATION OF A COMPLETE USABLE ALGORITHM.

!  REFERENCE
!     VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,

!     "COMPUTER GENERATION OF HYPERGEOMETRIC RANDOM VARIATES,"
!     JOURNAL OF STATISTICAL COMPUTATION AND SIMULATION,
!     22(1985), 2, 1985, 127-145.

!  REQUIRED SUBPROGRAMS
!     AFC() : A DOUBLE-PRECISION FUNCTION TO EVALUATE
!                THE LOGARITHM OF THE FACTORIAL.
!     RAND(): A UNIFORM (0,1) RANDOM NUMBER GENERATOR.

!  ARGUMENTS
!     NN1   : NUMBER OF WHITE BALLS          (INPUT)
!     NN2   : NUMBER OF BLACK BALLS          (INPUT)
!     KK    : NUMBER OF BALLS TO BE DRAWN    (INPUT)
!     ISEED : RANDOM NUMBER SEED  (INPUT AND OUTPUT)
!     JX    : NUMBER OF WHITE BALLS DRAWN   (OUTPUT)

!  STRUCTURAL VARIABLES
!     REJECT: LOGICAL FLAG TO REJECT THE VARIATE GENERATE BY H2PE.
!     SETUP1: LOGICAL FLAG TO SETUP FOR NEW VALUES OF NN1 OR NN2.
!     SETUP2: LOGICAL FLAG TO SETUP FOR NEW VALUES OF KK.
!     IX    : INTEGER CANDIDATE VALUE.
!     M     : DISTRIBUTION MODE.
!     MINJX : DISTRIBUTION LOWER BOUND.
!     MAXJX : DISTRIBUTION UPPER BOUND.
!     KS    : SAVED VALUE OF KK FROM THE LAST CALL TO H2PEC.
!     N1S   : SAVED VALUE OF NN1 FROM THE LAST CALL TO H2PEC.
!     N2S   : SAVED VALUE OF NN2 FROM THE LAST CALL TO H2PEC.
!     K,N1,N2: ALTERNATE VARIABLES FOR KK, NN1, AND NN2
!                (ALWAYS (N1 .LE. N2) AND (K .LE. (N1+N2)/2)).
!     TN    : TOTAL NUMBER OF WHITE AND BLACK BALLS

!  INVERSE-TRANSFORMATION VARIABLES
!     CON   : NATURAL LOGARITHM  OF SCALE.
!     P     : CURRENT SCALED PROBABILITY FOR THE INVERSE CDF.
!     SCALE : A BIG CONSTANT (1.E25) USED TO SCALE THE
!                PROBABILITY TO AVOID NUMERICAL UNDERFLOW
!     U     : THE UNIFORM VARIATE BETWEEN (0, 1.E25).
!     W     : SCALED HYPERGEOMETRIC PROBABILITY OF MINJX.

!  H2PE VARIABLES
!     S     : DISTRIBUTION STANDARD DEVIATION.
!     D     : HALF THE AREA OF THE RECTANGLE.
!     XL    : LEFT END OF THE RECTANGLE.
!     XR    : RIGHT END OF THE RECTANGLE.
!     A     : A SCALING CONSTANT.
!     KL    : HIGHEST POINT OF THE LEFT-TAIL REGION.
!     KR    : HIGHEST POINT OF THE RIGHT-TAIL REGION.
!     LAMDL : RATE FOR THE LEFT EXPONENTIAL TAIL.
!     LAMDR : RATE FOR THE RIGHT EXPONENTIAL TAIL.
!     P1    : AREA OF THE RECTANGLE.
!     P2    : AREA OF THE LEFT EXPONENTIAL TAIL PLUS P1.
!     P3    : AREA OF THE RIGHT EXPONENTIAL TAIL PLUS P2.
!     U     : A UNIFORM (0,P3) RANDOM VARIATE USED FIRST TO SELECT
!                ONE OF THE THREE REGIONS AND THEN CONDITIONALLY TO
!                GENERATE A VALUE FROM THE REGION.
!     V     : U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM
!                VALUE OR TO ACCEPT OR REJECT THE CANDIDATE VALUE.
!     F     : THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
!                ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL.
!     I     : INDEX FOR EXPLICIT CALCULATION OF F FOR H2PE.

! THE FOLLOWING VARIABLES ARE TEMPORARY VARIABLES USED IN COMPUTING THE UPPER
! AND LOWER BOUNDS OF THE NATURAL LOGARITHM OF THE SCALED DENSITY.
! THE DETAILED DESCRIPTION IS GIVEN IN PROPOSITIONS 2 AND 3 OF THE APPENDIX IN
! THE REFERENCE.
!     Y, Y1, YM, YN, YK, NK, R, S, T, E, G, DG, GU, GL, XM, XN, XK, NM

!     Y     : PRELIMINARY CONTINUOUS CANDIDATE VALUE, FLOAT(IX)
!     UB    : UPPER BOUND FOR THE NATURAL LOGARITHM OF THE SCALED DENSITY.
!     ALV   : NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V.
!     DR, DS, DT, DE: ONE OF MANY TERMS SUBTRACTED FROM THE UPPER
!                BOUND TO OBTAIN THE LOWER BOUND ON THE NATURAL
!                LOGARITHM OF THE SCALED DENSITY.
!     DELTAU: A CONSTANT, THE VALUE 0.0034 IS OBTAINED BY SETTING
!                N1 = N2 = 200, K = 199, M = 100, AND Y = 50 IN
!                THE FUNCTION DELTA_U IN LEMMA 1 AND ROUNDING THE
!                VALUE TO FOUR DECIMAL PLACES.
!     DELTAL: A CONSTANT, THE VALUE 0.0078 IS OBTAINED BY SETTING
!                N1 = N2 = 200, K = 199, M = 100, AND Y = 50 IN
!                THE FUNCTION DELTA_L IN LEMMA 1 AND ROUNDING THE
!                VALUE TO FOUR DECIMAL PLACES.

INTEGER, INTENT(IN)      :: kk
INTEGER, INTENT(IN)      :: nn1
INTEGER, INTENT(IN)      :: nn2
INTEGER, INTENT(IN OUT)  :: iseed1
INTEGER, INTENT(OUT)     :: jx

REAL (dp)      :: p, u, w, a, xl, xr
REAL           :: alv, d, dg, de, dr, ds, dt, e, f, g, gl, gu, kl, kr,  &
                  lamdl, lamdr, nk, nm, p1, p2, p3, r, s, t, ub, v,  &
                  xk, xm, xn, y, y1, yk, ym, yn
REAL, SAVE     :: tn
INTEGER        :: i, ix, k, n1, n2
LOGICAL, SAVE  :: reject, setup1, setup2
INTEGER, SAVE  :: ks = -1, n1s = -1, n2s = -1, m, minjx, maxjx
REAL (dp), PARAMETER  :: con = 57.56462733_dp, deltal = 0.0078_dp,  &
                         deltau = 0.0034_dp, scale = 1.0e+25_dp

!*****CHECK PARAMETER VALIDITY

IF ( nn1 < 0  .OR. nn2 < 0 .OR. kk  < 0 .OR. kk > nn1 + nn2 ) THEN
  jx     = -1
  RETURN
END IF
iseed = iseed1

!*****IF NEW PARAMETER VALUES, INITIALIZE

reject = .true.
setup1 = .false.
setup2 = .false.
IF (nn1 /= n1s .OR. nn2 /= n2s)  THEN
  setup1 = .true.
  setup2 = .true.
ELSE IF (kk /= ks)  THEN
  setup2 = .true.
END IF

IF (setup1)  THEN
  n1s   = nn1
  n2s   = nn2
  tn    = nn1 + nn2
  IF (nn1 <= nn2)  THEN
    n1 = nn1
    n2 = nn2
  ELSE
    n1 = nn2
    n2 = nn1
  END IF
END IF

IF (setup2)  THEN
  ks    = kk
  IF (kk+kk >= tn)  THEN
    k  = tn - kk
  ELSE
    k  = kk
  END IF
END IF

IF (setup1 .OR. setup2)  THEN
  m     = (k+1.) * (n1+1.) / (tn+2.)
  minjx = MAX (0, k-n2)
  maxjx = MIN (n1, k)
END IF

!*****GENERATE RANDOM VARIATE

IF (minjx == maxjx)  THEN
  
!        ...DEGENERATE DISTRIBUTION...
  
  ix      = maxjx
  RETURN
ELSE IF (m-minjx < 10)  THEN
  
!        ...INVERSE TRANSFORMATION...
  
  IF (setup1 .OR. setup2)  THEN
    IF (k < n2) THEN
      w = EXP (con + afc(n2) + afc(n1+n2-k) - afc(n2-k) - afc(n1+n2))
    ELSE
      w = EXP (con + afc(n1) + afc(k) - afc(k-n2) - afc(n1+n2))
    END IF
  END IF
  
  10 p  = w
  ix = minjx
  u  = rand () * scale
  20 IF (u > p)  THEN
    u  = u - p
    p  = p * (n1-ix)*(k-ix)
    ix = ix + 1
    p  = p / ix / (n2-k+ix)
    IF (ix > maxjx)  GO TO 10
    GO TO 20
  END IF
ELSE
  
!        ...H2PE...
  
  IF (setup1 .OR. setup2)  THEN
    s     = SQRT ((tn-k) * k * n1 * n2 / (tn-1) / tn /tn)
    
!           ...REMARK:  D IS DEFINED IN REFERENCE WITHOUT INT.
!           THE TRUNCATION CENTERS THE CELL BOUNDARIES AT 0.5
    
    d     = INT (1.5*s) + .5
    xl    = m - d + .5
    xr    = m + d + .5
    a     = afc(m) + afc(n1-m) + afc(k-m) + afc(n2-k+m)
    kl    = EXP (a - afc(INT(xl)) - afc(INT(n1-xl))  &
            - afc(INT(k-xl)) - afc(INT(n2-k+xl)))
    kr    = EXP (a - afc(INT(xr-1)) - afc(INT(n1-xr+1))  &
            - afc(INT(k-xr+1)) - afc(INT(n2-k+xr-1)))
    lamdl = -LOG (xl * (n2-k+xl) / (n1-xl+1) / (k-xl+1))
    lamdr = -LOG ((n1-xr+1) * (k-xr+1) / xr / (n2-k+xr))
    p1    = d + d
    p2    = p1 + kl / lamdl
    p3    = p2 + kr / lamdr
  END IF
  
  30 u    = rand () * p3
  v       = rand ()
  IF (u < p1)  THEN
    
!           ...RECTANGULAR REGION...
    
    ix    = xl + u
  ELSE IF (u <= p2)  THEN
    
!           ...LEFT TAIL...
    
    ix    = xl + LOG(v)/lamdl
    IF (ix < minjx)  GO TO 30
    v     = v * (u-p1) * lamdl
  ELSE
    
!           ...RIGHT TAIL...
    
    ix    = xr - LOG(v)/lamdr
    IF (ix > maxjx)  GO TO 30
    v     = v * (u-p2) * lamdr
  END IF
  
!        ...ACCEPTANCE/REJECTION TEST...
  
  IF (m < 100 .OR. ix <= 50)  THEN
    
!           ...EXPLICIT EVALUATION...
    
    f     = 1.0
    IF (m < ix)  THEN
      DO  i = m+1,ix
        f   = f * (n1-i+1) * (k-i+1) / (n2-k+i) / i
      END DO
    ELSE IF (m > ix)  THEN
      DO  i = ix+1,m
        f   = f * i * (n2-k+i) / (n1-i) / (k-i)
      END DO
    END IF
    IF (v <= f)  THEN
      reject = .false.
    END IF
  ELSE
    
!        ...SQUEEZE USING UPPER AND LOWER BOUNDS...
    
    y   = ix
    y1  = y + 1.
    ym  = y - m
    yn  = n1 - y + 1.
    yk  = k - y + 1.
    nk  = n2 - k + y1
    r   = -ym / y1
    s   = ym / yn
    t   = ym / yk
    e   = -ym / nk
    g   = yn * yk / (y1*nk) - 1.
    dg  = 1.
    IF (g < 0.) dg = 1.+g
    gu  = g * (1.+g*(-.5+g/3.))
    gl  = gu - .25 * (g*g)**2 / dg
    xm  = m + .5
    xn  = n1 - m + .5
    xk  = k - m + .5
    nm  = n2 - k + xm
    ub  = y * gu - m * gl + deltau + xm * r * (1.+r*(-.5+r/3.))  &
          + xn * s * (1.+s*(-.5+s/3.)) + xk * t * (1.+t*(-.5+t/3.))  &
          + nm * e * (1.+e*(-.5+e/3.))
    
!           ...TEST AGAINST UPPER BOUND...
    
    alv = LOG(v)
    IF (alv > ub)  THEN
      reject = .true.
    ELSE
      
!              ...TEST AGAINST LOWER BOUND...
      
      dr = xm * (r*r)**2
      IF (r < 0.)  dr = dr / (1.+r)
      ds = xn * (s*s)**2
      IF (s < 0.)  ds = ds / (1.+s)
      dt = xk * (t*t)**2
      IF (t < 0.)  dt = dt / (1.+t)
      de = nm * (e*e)**2
      IF (e < 0.)  de = de / (1.+e)
      IF (alv < ub-.25*(dr+ds+dt+de) +(y+m)*(gl-gu)-deltal) THEN
        reject = .false.
      ELSE
        
!                 ...STIRLING'S FORMULA TO MACHINE ACCURACY...
        
        IF (alv <= a - afc(ix) - afc(n1-ix) - afc(k-ix) - afc(n2-k+ix) ) THEN
          reject = .false.
        ELSE
          reject = .true.
        END IF
      END IF
    END IF
  END IF
  IF (reject)  GO TO 30
END IF

!*****RETURN APPROPRIATE VARIATE

IF (kk + kk >= tn)  THEN
  IF (nn1 > nn2)  THEN
    ix = kk - nn2 + ix
  ELSE
    ix =  nn1 - ix
  END IF
ELSE
  IF (nn1 > nn2)  ix = kk - ix
END IF
jx = ix

RETURN
END SUBROUTINE h2pec



FUNCTION afc(i) RESULT(fn_val)

!     FUNCTION TO EVALUATE LOGARITHM OF THE FACTORIAL I
!        IF (I .GT. 7), USE STIRLING'S APPROXIMATION
!        (It actually uses the less accurate De Moivre approximation)
!           OTHERWISE,  USE TABLE LOOKUP

INTEGER, INTENT(IN)  :: i
REAL (dp)            :: fn_val

REAL (dp) :: di

REAL (dp), PARAMETER  :: al(8) = (/ 0._dp, 0._dp, 0.6931471806_dp,  &
                         1.791759469_dp, 3.178053830_dp, 4.787491743_dp, &
                         6.579251212_dp, 8.525161361_dp /)

IF (i <= 7)  THEN
  fn_val = al(i+1)
ELSE
  di  = i
  fn_val = (di+0.5_dp) * LOG(di) - di + 0.08333333333333_dp/di  &
           - 0.00277777777777_dp/di/di/di  + 0.9189385332_dp
END IF

RETURN
END FUNCTION afc



FUNCTION rand() RESULT(fn_val)

!  UNIFORM RANDOM NUMBER GENERATOR
!  REFERENCE:  L. SCHRAGE,
!     "A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR,"
!     ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, 5(1979), 132-138.

REAL                :: fn_val

INTEGER             :: xhi, xalo, leftlo, fhi, k
INTEGER, PARAMETER  :: a = 16807, b15 = 32768, b16 = 65536, p = 2147483647

xhi    = iseed / b16
xalo   = (iseed-xhi * b16) * a
leftlo = xalo/b16
fhi    = xhi * a + leftlo
k      = fhi / b15
iseed  = (((xalo - leftlo * b16) - p) + (fhi - k * b15) * b16) + k
IF (iseed < 0) iseed = iseed + p
fn_val = iseed * 4.656612875E-10

RETURN
END FUNCTION rand

END MODULE toms668




PROGRAM test_random_hyperg
! TEST 1
! ------
! Test the random generation of hypergeometric variables by generating
! the number of spades in a hand of 13 cards out of 52.
! 1 million replicates used.

USE random_hypergeometric
USE toms668
IMPLICIT NONE

REAL     :: t1, t2
INTEGER  :: freq(0:13), c, d, h, i, j, s

WRITE(*, *) '  TEST 1'
WRITE(*, *)

! First using r_hyperg
freq = 0
CALL CPU_TIME(t1)
DO i = 1, 1000000
  j = r_hyperg(52, 13, 13)
  freq(j) = freq(j) + 1
END DO
CALL CPU_TIME(t2)
WRITE(*, '(a, f8.2, a)') ' Time taken by r_hyperg = ', t2-t1, 'secs.'
WRITE(*, '(a/9i7/5i7)') ' Frequencies', freq

! Now using h_alias
freq = 0
CALL CPU_TIME(t1)
DO i = 1, 1000000
  j = h_alias(52, 13, 13)
  freq(j) = freq(j) + 1
END DO
CALL CPU_TIME(t2)
WRITE(*, '(a, f8.2, a)') ' Time taken by h_alias = ', t2-t1, 'secs.'
WRITE(*, '(a/9i7/5i7)') ' Frequencies', freq

! Now using TOMS algorithm 668

freq = 0
CALL CPU_TIME(t1)
DO i = 1, 1000000
  CALL h2pec(13, 13, 39, iseed, j)
  freq(j) = freq(j) + 1
END DO
CALL CPU_TIME(t2)
WRITE(*, '(a, f8.2, a)') ' Time taken by TOMS668 = ', t2-t1, 'secs.'
WRITE(*, '(a/9i7/5i7)') ' Frequencies', freq

! TEST 2
! ------
! This tests the set-up time as the parameters change on every call.
! The number of clubs, diamonds, hearts and spades are found in a hand
! of 13 cards out of 52 by calling the random number generator 3 times
! with different parameters.

WRITE(*, *)
WRITE(*, *) '  TEST 2'
WRITE(*, *)

! First using r_hyperg
freq = 0
CALL CPU_TIME(t1)
DO i = 1, 100000
  c = r_hyperg(52, 13, 13)
  d = r_hyperg(39, 13, 13-c)
  h = r_hyperg(26, 13, 13-c-d)
  s = 13 - c - d - h
  freq(c) = freq(c) + 1
  freq(d) = freq(d) + 1
  freq(h) = freq(h) + 1
  freq(s) = freq(s) + 1
END DO
CALL CPU_TIME(t2)
WRITE(*, '(a, f8.2, a)') ' Time taken by r_hyperg = ', t2-t1, 'secs.'
WRITE(*, '(a/9i7/5i7)') ' Frequencies', freq

! Now using h_alias

freq = 0
CALL CPU_TIME(t1)
DO i = 1, 100000
  c = h_alias(52, 13, 13)
  d = h_alias(39, 13, 13-c)
  h = h_alias(26, 13, 13-c-d)
  s = 13 - c - d - h
  freq(c) = freq(c) + 1
  freq(d) = freq(d) + 1
  freq(h) = freq(h) + 1
  freq(s) = freq(s) + 1
END DO
CALL CPU_TIME(t2)
WRITE(*, '(a, f8.2, a)') ' Time taken by h_alias = ', t2-t1, 'secs.'
WRITE(*, '(a/9i7/5i7)') ' Frequencies', freq

! Now using TOMS algorithm 668

freq = 0
CALL CPU_TIME(t1)
DO i = 1, 100000
  CALL h2pec(13, 13, 39, iseed, c)
  CALL h2pec(13-c, 13, 26, iseed, d)
  CALL h2pec(13-c-d, 13, 13, iseed, h)
  s = 13 - c - d - h
  freq(c) = freq(c) + 1
  freq(d) = freq(d) + 1
  freq(h) = freq(h) + 1
  freq(s) = freq(s) + 1
END DO
CALL CPU_TIME(t2)
WRITE(*, '(a, f8.2, a)') ' Time taken by TOMS668 = ', t2-t1, 'secs.'
WRITE(*, '(a/9i7/5i7)') ' Frequencies', freq

STOP
END PROGRAM test_random_hyperg