NOTE: The following code is by George Russell
C A permutation on N objects {0,...,N-1} is represented in two ways.
C 1) As an INTEGER array (say PI) on N elements,PI(0) to PI(N-1). PI(I)
C is a
C number between 0 and N-1 and is the element obtained by applying
C the permutation to I.
C 2) As an INTEGER (say PIN). We shall not be able to deal with
C permutations on
C more than about 11 elements in this form.
C PIN is between 0 and N!-1. The conversion is not rapid
C (but it only has to be done a few times for each permutation, so
C this doesn't matter). PIN is the sum for I from 1 to N-1 of I!*
C T(I) where T(I) is a number between 0 and I. The T(I) are such
C that the permutation is the product of the transpositions
C (T(N-1) N-1)(T(N-2) N-2) . . . (T(1) 1) where
C for T(K)=K the resulting transposition is just the identity.
C
C The first thing we need are conversion routines both ways.
C
C ITOARR
C
C converts a permutation in form 2 to form 1. N,PIN and PI are as
C above. The permutation is read from PIN; converted; and written to
C PI.
SUBROUTINE ITOARR(PIN,PI)
C Include the common block (copied from BLOCK DATA subprogram).
***************** Include here
C (end of copied section).
INTEGER PIN,PI(0:N-1)
C
C Method: Starting with the identity permutation, postmultiply by
C (T(N-1) N-1), (T(N-2) N-2) in that order.
C Remainder, Quotient, Loop variable, Temporary variable
INTEGER REM,QUOT,I,TEMP
C
C Initialise to identity
DO 10 I=0,N-1
PI(I)=I
10 CONTINUE
C Get T(I) by repeated dividing current remainder by I!.
REM=PIN
DO 20 I=N-1,1,-1
QUOT=REM/FACT(I)
REM=MOD(REM,FACT(I))
C QUOT is T(I). Apply swap.
TEMP=PI(I)
PI(I)=PI(QUOT)
PI(QUOT)=TEMP
20 CONTINUE
END
C
C
C ARRTOI
C
C converts a permutation in form 1 to form 2. N,PIN and PI are as
C above. The permutation is read from PI; converted; and written to
C PIN.
SUBROUTINE ARRTOI(PIN,PI)
C Include the common block (copied from BLOCK DATA subprogram).
***************** Include here
C (end of copied section).
INTEGER PIN,PI(0:N-1)
C
C Method: (1) Copy PI to PI2; let INVPI2 be the inverse of PI2
C (2) The inverse
C is (T(1) 1)(T(2) 2) . . . (T(N-1) N-1). For I from N-1 to 1,
C we change INVPI2 and PI2 gradually to the identity, keeping the
C two inverse to each other and using each other to undo the
C transpositions., undoing each transposition
C in turn.
C
INTEGER INVPI2(0:MAXN-1),PI2(0:MAXN-1),I,TEMP
C
C TI is T(I); INVI is what maps to I under PI.
INTEGER TI,INVI
C
C Initialise INVPI2 and PI2.
DO 10 I=0,N-1
PI2(I)=PI(I)
INVPI2(PI(I))=I
10 CONTINUE
C
C Accumulate sum of T(I)*I! in PIN.
PIN=0
DO 20 I=N-1,1,-1
TI=PI2(I)
INVI=INVPI2(I)
PIN=PIN+TI*FACT(I)
C Update INVPI2 and PI2.
TEMP=INVPI2(I)
INVPI2(I)=INVPI2(TI)
INVPI2(TI)=TEMP
TEMP=PI2(I)
PI2(I)=PI2(INVI)
PI2(INVI)=TEMP
20 CONTINUE
END