MODULE hypergeometric

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60)
INTEGER, SAVE      :: ier
PUBLIC             :: hyperg, lnfact, lngamma, ier

CONTAINS

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)

!     Programmer: Alan Miller
!     http://users.bigpond.net.au
!     e-mail: amiller @ bigpond.net.au
!     Latest revision - 22 May 1997

INTEGER, INTENT(IN)  :: n(:,:)
REAL (dp)     :: prob

!     Local variables

INTEGER          :: ncell(1:4), ntot(1:4), nall
LOGICAL          :: ended
REAL (dp) :: zero = 0.d0, one = 1.d0

!     Check that all elements are positive or zero.

ier = 0
IF (n(1,1) < 0 .OR. n(1,2) < 0 .OR. n(2,1) < 0 .OR. n(2,2) < 0) THEN
  ier = 1
  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
  ier = 2
  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.

IF (ncell(2) < 10) THEN                ! Use divisions & multiplications
  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
ELSE                                   ! Use logs of factorials for speed
  prob = EXP( lnfact(ntot(1)) + lnfact(ntot(2)) + lnfact(ntot(3)) +   &
              lnfact(ntot(4)) - lnfact(nall) - lnfact(ncell(1)) -     &
              lnfact(ncell(2)) - lnfact(ncell(3)) - lnfact(ncell(4)) )
END IF

RETURN
END FUNCTION hyperg



SUBROUTINE order4(n)

!     Order 4 integers; largest first.

INTEGER, INTENT(IN OUT) :: n(:)

!     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



FUNCTION lnfact(n) RESULT(fn_val)
! Calculate the logarithm of n!

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

! Local variable
REAL (dp) :: zero = 0.0D0

SELECT CASE (n)
  CASE (0:1)
    fn_val = zero
  CASE DEFAULT
    fn_val = lngamma(DBLE(n+1))
END SELECT

RETURN
END FUNCTION lnfact



FUNCTION lngamma(z) RESULT(lanczos)

!  Uses Lanczos-type approximation to ln(gamma) for z > 0.
!  Reference:
!       Lanczos, C. 'A precision approximation of the gamma function',
!               J. SIAM Numer. Anal., B, 1, 86-96, 1964.
!  Accuracy: About 14 significant digits except for small regions
!            in the vicinity of 1 and 2.

!  Programmer: Alan Miller
!              1 Creswick Street, Brighton, Vic. 3187, Australia
!  Latest revision - 14 October 1996

REAL(dp), INTENT(IN) :: z
REAL(dp)             :: lanczos

! Local variables

REAL(dp)  :: a(9) = (/ 0.9999999999995183D0, 676.5203681218835D0, &
		       -1259.139216722289D0, 771.3234287757674D0, &
		       -176.6150291498386D0, 12.50734324009056D0, &
		       -0.1385710331296526D0, 0.9934937113930748D-05, &
			0.1659470187408462D-06 /), zero = 0.D0,   &
			one = 1.d0, lnsqrt2pi = 0.9189385332046727D0, &
			half = 0.5d0, sixpt5 = 6.5d0, seven = 7.d0, tmp
INTEGER   :: j

IF (z <= zero) THEN
  WRITE(*, *)'Error: zero or -ve argument for lngamma'
  RETURN
END IF

lanczos = zero
tmp = z + seven
DO j = 9, 2, -1
  lanczos = lanczos + a(j)/tmp
  tmp = tmp - one
END DO
lanczos = lanczos + a(1)
lanczos = LOG(lanczos) + lnsqrt2pi - (z + sixpt5) + (z - half)*LOG(z + sixpt5)
RETURN

END FUNCTION lngamma

END MODULE hypergeometric



PROGRAM t_hyperg

!     Test function hyperg for calculating hypergeometric probabilities

USE hypergeometric
IMPLICIT NONE
INTEGER  :: n(2,2)

DO
  WRITE(*, *)'Enter 4 integers: '
  READ(*, *) n
  IF (ier == 0) THEN
    WRITE(*, 900) hyperg(n)
    900 FORMAT(' Prob. = ', g14.6/)
  ELSE
    WRITE(*, *) 'Error in input values'
  END IF
END DO

STOP
END PROGRAM t_hyperg