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