C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/combenum.for
C
C===Enumerate Combinations of an N-Vector         HDK 3/19/93
C To run in WATFOR77 issue: WATFOR77 COMBENUM
C   or in VS Fortran issue: FVCG COMBENUM
C
C   In this example V is a 3-vector of binary values.
      PARAMETER (N=3)
      INTEGER V(N), IC(N),BINOM,NCOMB,I,J,K
      DATA V/N*0/
C
      WRITE(6,*) 'V=',(V(K),K=1,N)
C---There are clearly  2**N - 1  possible non-zero binary values.
      DO 4 J=1,N
        NCOMB=BINOM(N,J)
        DO 3 I=1,NCOMB
          CALL COMB(N,J,I,IC)
          DO 1 K=1,N
1         V(K)=0
          DO 2 K=1,J
2         V(IC(K))=1
3       WRITE(6,*) 'V=',(V(K),K=1,N)
4     CONTINUE
      STOP
      END
      SUBROUTINE COMB(N, P, L, C)
C ACM Algorithm 515, by B.P. Buckles and M. Lybanon
C "Generation of a Vector from the Lexicographical Index"
C T.O.M.S. Vol. 3, No. 2, pp 180-182, June 1977.
C
C THIS SUBROUTINE FINDS THE COMBINATION SET OF N THINGS
C TAKEN P AT A TIME FOR A GIVEN LEXICOGRAPHICAL INDEX.
C N - NUMBER OF THINGS IN THE SET
C P - NUMBER OF THINGS IN EACH COMBINATION
C L - LEXICOGRAPHICAL INDEX OF COMBINATION SOUGHT
C C - OUTPUT ARRAY CONTAINING THE COMBINATION SET
C THE FOLLOWING RELATIONSHIPS MUST EXIST AMONG THE INPUT
C VARIABLES.  L MUST BE GREATER THAN OR EQUAL TO 1 AND LESS
C THAN OR EQUAL TO THE MAXIMUM LEXICOGRAPHICAL INDEX.
C P MUST BE LESS THAN OR EQUAL TO N AND GREATER THAN ZERO.
      INTEGER N, P, L, C(P), K, R, P1, BINOM
C INITIALIZE LOWER BOUND INDEX AT ZERO
      K = 0
C LOOP TO SELECT ELEMENTS IN ASCENDING ORDER
      P1 = P - 1
      C(1)=L
      IF(P1.LE.0) RETURN
      DO 20 I=1,P1
C SET LOWER BOUND ELEMENT NUMBER FOR NEXT ELEMENT VALUE
        C(I) = 0
        IF (I.NE.1) C(I) = C(I-1)
C LOOP TO CHECK VALIDITY OF EACH ELEMENT VALUE
   10   C(I) = C(I) + 1
        R = BINOM(N-C(I),P-I)
        K = K + R
        IF (K.LT.L) GO TO 10
        K = K - R
   20 CONTINUE
      C(P) = C(P1) + L - K
      RETURN
      END
      INTEGER FUNCTION BINOM(M, N)
C ACM ALGORITHM 160 TRANSLATED TO FORTRAN.  CALCULATES THE
C NUMBER OF COMBINATIONS OF M THINGS TAKEN N AT A TIME.
      INTEGER M, N, P, I, N1, R
      N1 = N
      P = M - N1
      IF (N1.GE.P) GO TO 10
      P = N1
      N1 = M - P
   10 R = N1 + 1
      IF (P.EQ.0) R = 1
      IF (P.LT.2) GO TO 30
      DO 20 I=2,P
        R = (R*(N1+I))/I
   20 CONTINUE
   30 BINOM = R
      RETURN
      END