C Test LPG (Lexicographic permutation generation)
C including operation counts (how many times are the steps of the
C permutation routine executed to create all permutations of
C n distinct elements)
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Version history:
C 14.01.2003 Asymptotic values of counts added in comments
C 12.01.2003 Counters for element comparisons and index tests
C 09.01.2003 Counters yielding obvious results removed
C 08.01.2003 Initial version
C
C Diagnostic counts
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      INTEGER
     & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
     & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
      COMMON /DIAGNS/
     & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
     & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
C LCOMA counts all comparisons between array elements of A
C LCOMJ counts tests of array indices
      INTEGER LCOMA, LCOMJ, LCOMAO, LCOMJO
      COMMON /DIAGNC/ LCOMA, LCOMJ
C
      INTEGER A(20)
C
      OPEN ( UNIT=10, FILE='LPGTES.LIS', FORM='FORMATTED',
     &       STATUS='UNKNOWN', RECL=160 )
C Header line
      WRITE (10,1001)
C
C Loop over all cases (limitation to N=13 due to 32 bit counts)
      DO 10 N = 4, 13
C      L2EAS=0
C      L21TES=0
C Preset diagnostic counts
      L22XGE=0
      L22TES=0
      L22YGE=0
      L22YGT=0
      L2XLTZ=0
      L2XGEZ=0
      L3TES=0
      L3EASI=0
      L31TES=0
      L31LOP=0
      L31YGE=0
      L41REV=0
      L4LOP=0
      L4KLTL=0
C
      LCOMA = 0
      LCOMAO = 0
      LCOMJ = 0
      LCOMJO = 0
C Preset array to be permuted
      DO 20 I = 1, N
      A(I) = I
20    CONTINUE
C
C      IFAC = 0
30    CONTINUE
C      IFAC = IFAC + 1
C LPG returns NEXT=0 after N! calls (for distinct A(i))
C (IFAC exceeds range 2*31-1 for N=13)
      CALL LPG ( N, A, NEXT)
CDIAG
      IF ( LCOMA .GT. 1000000000 ) THEN
        LCOMAO = LCOMAO + 1
        LCOMA = LCOMA - 1000000000
      ENDIF
      IF ( LCOMJ .GT. 1000000000 ) THEN
        LCOMJO = LCOMJO + 1
        LCOMJ = LCOMJ - 1000000000
      ENDIF
CDIAG
C Terminate when reverse order is reached
      IF ( NEXT .EQ. 0 ) GOTO 40
      GOTO 30
40    CONTINUE
C Output of diagnostic counts 
      WRITE (10,1000) N,
     & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
     & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
1000  FORMAT ( I3, 14I11 )
1001  FORMAT ( '  N',5X,'L22XGE', 5X,'L22TES', 5X,'L22YGE',
     &               5X,'L22YGT', 5X,'L2XLTZ', 5X,'L2XGEZ',
     &               5X,' L3TES', 5X,'L3EASI', 5X,'L41TES',
     &               5X,'L31LOP', 5X,'L31YGE', 5X,'L41REV',
     &               5X,' L4LOP', 5X,'L4KLTL')
      WRITE (*,1002) N, LCOMAO, LCOMA, LCOMJO, LCOMJ
1002  FORMAT (I3, I4, I9.9, I4, I9.9)
10    CONTINUE
      END
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      SUBROUTINE LPG ( N, A, NEXT )
C
C Generate next permutation of N elements A(i) in
C lexicographic order
C
C The method used is
C "Algorithm L (Lexicographic permutation generation)"
C described in Chapter 7.2.1.2
C of Donald E. Knuth's The Art of Computer Programming, Volume 4,
C Combinatorial Algorithms, C Volume 4A, Enumeration and Backtracking.
C
C with the modification proposed in the answer to
C exercise 1, which is based on "Algorithm 28,
C PERMUTATIONS OF THE ELEMENTS OF A VECTOR IN LEXICOGRAPHIC ORDER"
C by J. P. N. Phillips. The Computer Journal, Volume 10,
C Issue 3: November 1967. (Algorithms supplement),
C also available online:
C http://www3.oup.co.uk/computer_journal/hdb/Volume_10/Issue_03/default.html
C
C The printed version of TAOCP Vol.4 is expected to be available
C in the year 2007. D. Knuth has put some chapters online, including
C Pre-fascicle 2b (generating all permutations), available at
C http://www-cs-faculty.stanford.edu/~knuth/fasc2b.ps.gz
C Most of the comments are copied verbally from Knuth's
C description
C
C This subroutine is a special version including counters for
C all paths of the algorithm. The main result found is:
C
C To create all permutations of n distinct elements the
C number of comparisons between the array elements becomes
C 2.410756*n! for large n (e.g. n>8)
C The exact number of comparisons c(n) is given by the recursion
C c(3)=11, c(n)=n*c(n-1)+n*(n+1)/2 for n>=4
C OEIS: A079884(n) = 2*n! - 1 + A079750(n) + A079753(n)
C c(3)...c(14): 11,54,285,1731,12145,97196,874809,8748145,
C               96229661,1154756010,15011828221,210165595199
C
C The required number of index tests (test for termination and
C test in the final reversion loop) becomes
C 0.2613625*n! for large n, if the test for n=3 is excluded.
C The exact number of index comparisons i(n) is given by the recursion
C i(3)=0, i(n)=n*i(n-1)+1+(n-1)*floor((n-1)/2) for n>=4
C OEIS: A079885(n) = A079751(n) + A079755(n)
C i(3)...i(14): 0,4,29,185,1314,10534,94839,948427,10432748,
C               125193032,1627509489,22785132925
C If n=3 is included the additionally required termination test
C adds n!/6 index comparisons, increasing the number of index
C comparisons to 0.428029*n! (63.8% more index comparisons).
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Change history:
C 10.01.2003 Results from diagnostic counts included as comments
C 08.01.2003 Diagnostic counts
C 06.01.2003 Optional code to treat N<4
C 05.01.2003 Initial version
C
C Call parameters:
      INTEGER N, A(N), NEXT
C N: number of elements to permute. N>=4 must hold (not checked)
C A: Integer array holding the elements.
C    to be preset externally (e.g. by 1,2,3,...,n-1,n)
C    A is overwritten by the lexicographically next permutation
C NEXT: =1 if A has been successfully updated
C       =0 if A was already in reverse lexicographic order
C          i.e. if a(n)>=a(n-1)>=....>=a(2)>=a(1)
C          In this case no change is made to A
C
C Local variables
      INTEGER J, K, L, X, Y, Z, H
C
C Diagnostic counts 
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      INTEGER
     & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
     & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
      COMMON /DIAGNS/
     & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
     & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
      INTEGER LCOMA, LCOMJ
      COMMON /DIAGNC/ LCOMA, LCOMJ

C
C Preset return status
      NEXT = 1
C
C In a generally applicable version for all N all lines starting
C with comment CS have to be activated
CSC Start part necessary for N<3
CS      IF ( N .LE. 1 ) THEN
CS        NEXT = 0
CS        RETURN
CS      ENDIF
CS      IF ( N .EQ. 2 ) THEN
CS        IF ( A(1) .LT. A(2) ) THEN
CS          H = A(1)
CS          A(1) = A(2)
CS          A(2) = H
CS        ELSE
CS          NEXT = 0
CS        ENDIF
CS        RETURN
CS      ENDIF
CSC End part necessary for N<3
C
C Find the largest j such that a(j) can be increased
C
C Step L2 finds j=n-1 half of the time when the elements are distinct,
C because exactly n!/2 of the n! permutations have a(n-1)<a(n)
C
C L2' Easiest case?
      Y = A(N-1)
      Z = A(N)
CDIAG
      LCOMA = LCOMA + 1
CDIAG
      IF ( Y .LT. Z ) THEN
CDIAG (Checked to be n!/2, then removed)
C        L2EAS = L2EAS + 1
CDIAG
        A(N-1) = Z
        A(N) = Y
        RETURN
      ENDIF
C
C L2.1' Next easiest case?
      X = A(N-2)
CDIAG (Checked to be n!/2, then removed)
C      L21TES = L21TES + 1
      LCOMA = LCOMA + 1
CDIAG
      IF ( X .GE. Y ) THEN
C L2.2' Find j
        J = N - 3
CDIAG
        L22XGE = L22XGE + 1
CDIAG L22XGE_sum = n!/6
CDIAG
C Treat special case N=3
CSC Start part necessary for N=3
CN>=4  LCOMJ = LCOMJ + 1
CS       IF ( J .EQ. 0 ) THEN
CS         NEXT = 0
CS         RETURN
CS       ENDIF
CSC End part necessary for N=3
        Y = A(J)
10      CONTINUE
CDIAG
        L22TES = L22TES + 1
CDIAG L22TES_sum = 0, 4, 25, 156, 1099, 8800, 79209, 792100, 8713111,
CDIAG              104557344, 1359245485  -> 0.2182818*n! for large n
CDIAG for n  = 3, 4, .... , 13  OEIS A079750
CDIAG a079750(3)=0, a079750(n) = n * a079750 + n for n>=4
        LCOMA = LCOMA + 1
CDIAG
        IF ( Y .GE. X ) THEN
CDIAG
          L22YGE = L22YGE + 1
CDIAG L22YGE_sum = 0, 1, 6, 37, 260, 2081, 18730, 187301, 2060312,
CDIAG              24723745, 321408686  -> 0.0516151*n! for large n
CDIAG for n  = 3, 4, .... , 13  OEIS A079751
CDIAG a079751(3)=0, a079751(n) = n * a079751(n-1) + 1 for n>=4
CDIAG
          J = J - 1
C Reverse lexicographic order reached?
CDIAG
          LCOMJ = LCOMJ + 1
CDIAG
          IF ( J .EQ. 0 ) THEN
C Terminate
            NEXT = 0
            RETURN
          ENDIF
          X = Y
          Y = A(J)
          GOTO 10
        ENDIF
C
C At this point j is the largest subscript such that we have already
C visited all permutations beginning with a(1)...a(j). Therefore
C the lexicographically next permutation will increase the value
C of a(j).
C
      ELSE
C
C Test for next easiest case passed, permute last 3 elements
C
CDIAG
        L22YGT = L22YGT + 1
CDIAG L22YGT_sum = n!/3
        LCOMA = LCOMA + 1
CDIAG
        IF ( X .LT. Z ) THEN
CDIAG
          L2XLTZ = L2XLTZ + 1
CDIAG L22XLTZ_sum = n!/6
CDIAG
          A(N-2) = Z
          A(N-1) = X
          A(N)   = Y
        ELSE
CDIAG
          L2XGEZ = L2XGEZ + 1
CDIAG L22XGEZ_sum = n!/6
CDIAG
          A(N-2) = Y
          A(N-1) = Z
          A(N)   = X
        ENDIF
        RETURN
      ENDIF
C
C L3 Increase a(j) by the smallest feasible amount
C
C If a(j) >= a(l) decrease l repeatedly until a(j)<a(l).
C Then interchange a(j)<->a(l)
C
C L3' Easy increase?
CDIAG
      L3TES = L3TES + 1
CDIAG L3TES_sum = n!/6 - 1
      LCOMA = LCOMA + 1
CDIAG
      IF ( Y .LT. Z ) THEN
CDIAG
        L3EASI = L3EASI + 1
CDIAG L3EASI_sum = L22YGE_sum (OEIS A079751)
        A(J)   = Z
        A(J+1) = Y
        A(N)   = X
C Done, go to L4.1'
        GOTO 41
      ENDIF
C L3.1' Increase A(J)
      L = N - 1
CDIAG
      L31TES = L31TES + 1
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
CDIAG L31TES_sum = 0, 2, 13, 82, 579, 4638, 41749, 417498, 4592487,
CDIAG              55109854, 716428113  -> 0.1150515*n! for large n
CDIAG for n  = 3, 4, .... , 13  OEIS A079752
CDIAG a079752(3)=0, a079752(n) = n * a079752(n-1) + n - 2 for n>=4
CDIAG
20    CONTINUE
CDIAG
      L31LOP = L31LOP + 1
CDIAG L31LOP_sum = 0, 3, 21, 136, 967, 7757, 69841, 698446, 7682951,
CDIAG              92195467, 1198541137   -> 0.1924742*n! for large n
CDIAG for n  = 3, 4, .... , 13  OEIS A079753
CDIAG a079753(3)=0, a079753(n)=n*a079753(n-1) + (n-1)*(n-2)/2 for n>=4
CDIAG a079753(n)=a079752(n)+a079754(n)
      LCOMA = LCOMA + 1
CDIAG
      IF ( Y .GE. A(L) ) THEN
CDIAG
      L31YGE = L31YGE + 1
CDIAG L32YGE_sum = 0, 1, 8, 54, 388, 3119, 28092, 280948, 3090464,
CDIAG              37085613, 482113024  -> 0.0774227*n! for large n
CDIAG for n  = 3, 4, .... , 13  OEIS A079754
CDIAG a079754(3)=0, a079754(n)=n*a079754(n-1) + (n-2)*(n-3)/2 for n>=4
CDIAG
        L = L - 1
        GOTO 20
      ENDIF
C
C Since a(j+1) >= ... >= a(n), element a(l) is the smallest element
C greater than a(j) that can legitimately follow a(1)...a(j-1) in a
C permutation. Before the interchange we have a(j+1)>=...>=a(l-1)>=a(l)
C >a(j)>=a(l+1)>=...>=a(n)
C
C Interchange
C
      A(J) = A(L)
      A(L) = Y
C
C After the interchange, we have
C a(j+1)>=...>=a(l-1)>=a(j)>a(l)>=a(l+1)>=...>=a(n)
C
C Find the lexicographically least way to extend the new a(1)...a(j)
C to a complete pattern:
C The first permutation beginning with the current prefix a(1)...a(j)
C is a(1)...a(j)a(n)...a(j+1), and step L4 produces it by doing
C floor (n-j)/2 interchanges
C
C L4 Reverse a(j+1)...a(n)
C L4' Begin to reverse
      A(N) = A(J+1)
      A(J+1) = Z
C L4.1' Reverse a(j+1)...a(n-1)
41    CONTINUE
CDIAG
      L41REV = L41REV + 1
CDIAG L41REV_sum = n!/6 - 1
      K = J + 2
      L = N - 1
30    CONTINUE
CDIAG
      L4LOP = L4LOP + 1
CDIAG L4LOP_sum = 0, 3, 23, 148, 1054, 8453, 76109, 761126, 8372436,
CFIAG             100469287, 1306100803  -> 0.2097473*n! for large n
CDIAG    for n  = 3, 4, .... , 13  OEIS A079755
CDIAG a079755(3) = 0
CDIAG a079755(n) = n * a079755(n-1) + (n-1)*floor((n-1)/2) for n>=4
      LCOMJ = LCOMJ + 1
CDIAG
      IF ( K .LT. L ) THEN
CDIAG
        L4KLTL = L4KLTL + 1
CDIAG L4KLTL_sum = 0, 0, 4, 29, 215, 1734, 15360, 156327, 1719637,
CDIAG              20635688, 268264004  -> 0.0430806*n! for large n
CDIAG     for n  = 3, 4, .... , 13  OEIS A079756
CDIAG a079756(4) = 0
CDIAG a079756(n) = n*a079756(n-1) + (n-1)*(floor((n-1)/2)-1) for n>=5
CDIAG
        H = A(K)
        A(K) = A(L)
        A(L) = H
        K = K + 1
        L = L - 1
        GOTO 30
      ENDIF
      RETURN
C End of SUBROUTINE LPG
      END