SUBROUTINE ncr(n, r, ncomb, ier)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2000-01-20  Time: 18:08:52

!     Calculate the number of different combinations of r objects out of n.

!     ier = 0 if no error is detected
!         = 1 if n < 1
!         = 2 if r < 0
!         = 3 if r > n
!         = 4 if nCr > 1.e+308, i.e. if it overflows.  In this case, the
!                natural log of nCr is returned.

!     Programmer: Alan.Miller @ cmis.csiro.au
!     Latest revision - 28 July 1988 (Fortran 77 version)

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

INTEGER, INTENT(IN)     :: n
INTEGER, INTENT(IN)     :: r
REAL (dp), INTENT(OUT)  :: ncomb
INTEGER, INTENT(OUT)    :: ier

INTEGER :: rr, i, nn

INTERFACE
  FUNCTION lngamma(x) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)  :: x
    REAL (dp)              :: fn_val
  END FUNCTION lngamma
END INTERFACE

IF (n < 1) THEN
  ier = 1
ELSE IF (r < 0) THEN
  ier = 2
ELSE IF (r > n) THEN
  ier = 3
ELSE
  ier = 0
END IF
IF (ier /= 0) RETURN

IF (r <= n-r) THEN
  rr = r
ELSE
  rr = n - r
END IF

IF (rr == 0) THEN
  ncomb = 1.0_dp
  RETURN
END IF

IF (rr > 25) THEN
  ncomb = lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1))
  IF (ncomb > 709._dp) THEN
    ier = 4
  ELSE
    ncomb = EXP(ncomb)
  END IF
  RETURN
END IF

ncomb = n
i = 1
nn = n
DO
  IF (i == rr) RETURN
  nn = nn - 1
  i = i + 1
  ncomb = (ncomb * nn) / REAL(i)
END DO

RETURN
END SUBROUTINE nCr



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

INTERFACE
  SUBROUTINE ncr(n, r, ncomb, ier)
    IMPLICIT NONE
    INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
    INTEGER, INTENT(IN)     :: n
    INTEGER, INTENT(IN)     :: r
    REAL (dp), INTENT(OUT)  :: ncomb
    INTEGER, INTENT(OUT)    :: ier
  END SUBROUTINE ncr
END INTERFACE

INTEGER    :: n, r, ier
REAL (dp)  :: result

DO
  WRITE(*, '(a)', ADVANCE='NO') ' Enter n, r : '
  READ(*, *) n, r
  CALL nCr(n, r, result, ier)
  IF (ier /= 0) THEN
    WRITE(*, *) ' Error, IER = ', ier
    IF (ier == 4) WRITE(*, '(a, f12.5)') ' Ln(nCr) = ', result
  ELSE
    WRITE(*, '(a, g16.8)') ' nCr = ', result
  END IF
END DO

STOP
END PROGRAM test_nCr



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

IMPLICIT NONE
INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)
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