random_neg_binomial Function

FUNCTION random_neg_binomial(sk, p) RESULT(ival)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
  ! INVERSION AND/OR THE REPRODUCTIVE PROPERTY.
  !    SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
  !       = the `power' parameter of the negative binomial
  !           (0 < REAL)
  !    P = BERNOULLI SUCCESS PROBABILITY
  !           (0 < REAL < 1)
  ! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
  ! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
  ! THE REPRODUCTIVE PROPERTY IS USED.
  REAL, INTENT(IN)   :: sk, p
  INTEGER            :: ival
  !     Local variables
  ! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER).
  REAL, PARAMETER    :: h = 0.7
  REAL               :: q, x, st, uln, v, r, s, y, g
  INTEGER            :: k, i, n
  IF (sk <= zero .OR. p <= zero .OR. p >= one) THEN
     WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
     STOP
  END IF
  q = one - p
  x = zero
  st = sk
  IF (p > h) THEN
     v = one/LOG(p)
     k = st
     DO i = 1,k
        DO
           CALL RANDOM_NUMBER(r)
           IF (r > zero) EXIT
        END DO
        n = v*LOG(r)
        x = x + n
     END DO
     st = st - k
  END IF
  s = zero
  uln = -LOG(vsmall)
  IF (st > -uln/LOG(q)) THEN
     WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK'
     STOP
  END IF
  y = q**st
  g = st
  CALL RANDOM_NUMBER(r)
  DO
     IF (y > r) EXIT
     r = r - y
     s = s + one
     y = y*p*g/s
     g = g + one
  END DO
  ival = x + s + half
  RETURN
END FUNCTION random_neg_binomial