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..
C INTEGER
C & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
C & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
C COMMON /DIAGNS/
C & L22XGE, L22TES, L22YGE, L22YGT, L2XLTZ, L2XGEZ, L3TES,
C & L3EASI, L31TES, L31LOP, L31YGE, L41REV, L4LOP, L4KLTL
C INTEGER LCOMA, LCOMJ
C 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)=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
C L22TES = L22TES + 1
CDIAG L22TES_sum = 0, 4, 25, 156, 1099, 8800, 79209, 792100, 8713111,
CDIAG 104557344, 1359245485
CDIAG for n = 3, 4, .... , 13 OEIS A079750
CDIAG a079750(3)=0, a079750(n) = n * a079750 + n for n>=4
C LCOMA = LCOMA + 1
CDIAG
IF ( Y .GE. X ) THEN
CDIAG
C L22YGE = L22YGE + 1
CDIAG L22YGE_sum = 0, 1, 6, 37, 260, 2081, 18730, 187301, 2060312,
CDIAG 24723745, 321408686
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
C 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
C L22YGT = L22YGT + 1
CDIAG L22YGT_sum = n!/3
C LCOMA = LCOMA + 1
CDIAG
IF ( X .LT. Z ) THEN
CDIAG
C L2XLTZ = L2XLTZ + 1
CDIAG L22XLTZ_sum = n!/6
CDIAG
A(N-2) = Z
A(N-1) = X
A(N) = Y
ELSE
CDIAG
C 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
C L3' Easy increase?
CDIAG
C L3TES = L3TES + 1
CDIAG L3TES_sum = n!/6 - 1
C LCOMA = LCOMA + 1
CDIAG
IF ( Y .LT. Z ) THEN
CDIAG
C 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
C 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
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
C L31LOP = L31LOP + 1
CDIAG L31LOP_sum = 0, 3, 21, 136, 967, 7757, 69841, 698446, 7682951,
CDIAG 921195467, 1198541137
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)
C LCOMA = LCOMA + 1
CDIAG
IF ( Y .GE. A(L) ) THEN
CDIAG
C L31YGE = L31YGE + 1
CDIAG L32YGE_sum = 0, 1, 8, 54, 388, 3119, 28092, 280948, 3090464,
CDIAG 37085613, 482113024
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
C L41REV = L41REV + 1
CDIAG L41REV_sum = n!/6 - 1
K = J + 2
L = N - 1
30 CONTINUE
CDIAG
C L4LOP = L4LOP + 1
CDIAG L4LOP_sum = 0, 3, 23, 148, 1054, 8453, 76109, 761126, 8372436,
CFIAG 100469287, 1306100803
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
C LCOMJ = LCOMJ + 1
CDIAG
IF ( K .LT. L ) THEN
CDIAG
C L4KLTL = L4KLTL + 1
CDIAG L4KLTL_sum = 0, 0, 4, 29, 215, 1734, 15360, 156327, 1719637,
CDIAG 20635688, 268264004
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