1 ################################################################################
   2 
   3 # Generate bag-permutations of m elements, ordered by nearly adjacent
   4 # transposition; generator value = m+1 on termination;
   5 # assumes m > 0; 4 sentinels; no pre-initialisation.
   6 # Maple program, Fred Lunnon, Maynooth 30/10/08
   7 # Version with type-checking time 7x slower than without under Maple 9
   8 # --- commented out again!
   9 # A 'run' is defined to be a sequence of equally ranked elements in
  10 # consecutive locations (and with consecutive identifiers), such that all
  11 # but the righthand end (highest location) are drifting right (d = +1).
  12 # The drift of the righthand element at j determines that of the entire run.
  13 
  14 # A run is 'blocked' when the element adjacent in its direction of drift at i
  15 # has equal or lower rank. The element adjacent to a run at its lefthand end
  16 # is of different rank; at its righthand end may be another of equal rank.
  17 # Generator employs Steinhaus-Johnson-Trotter adjacent transposition strategy
  18 # to exchange the earliest (by identifier) unblocked run at j with the single
  19 # element at i. Individual elements within a run behave as though under the
  20 # action of a combination generator by nearly-adjacent transpositions;
  21 # all permutations of the user bag of ranks are transpositions, most of
  22 # which are actually adjacent.
  23 # Once an unblocked element is found, the attached run is shifted one place
  24 # along direction d, the end elements transposed, and all but the rightmost
  25 # drift set leftwards; the drift of all earlier elements is reversed.
  26 # In order to avoid unnecesary and nested loops, reversal takes place during
  27 # search, and the left end of a run is retained in l for future reference.
  28 # Spatial blocking sentinels ranked lowest are installed at either end, and
  29 # temporal termination sentinels ranked highest placed off the righthand end.
  30 
  31 # Case m = 0 is excluded, since iterator signals immediate termination;
  32 # could be fixed, at the cost of extra testing at its start.
  33 # The generator is stable, in the sense that elements of the same rank remain
  34 # always in the same order. For this reason it might be regarded as a
  35 # permutation generator with restrictions on the ordering of specified subsets.
  36 # Next move takes place in constant amortised time; always by transposition,
  37 # which will be adjacent when possible (for large m, at least 82.7% of moves).
  38 
  39 # Each element has unique identifier p, lower p drifting more frequently;
  40 # R_[k] = rank (positive integer user symbol) of element at k;
  41 # P_[k] = identifier p of element at k, 0 <= p < m+3;
  42 # Q_[p] = location k of element p, 0 <= k < m+3;
  43 # D_[k] = current drift direction d of element at k, -1 <= d <= +1.
  44 # i,j,k,l locate elements in current permutation;
  45 # d,e = direction of element at j,i; l holds left end of current run.
  46 # Arrays are extended by sentinels indexed 0 at left, m+1,m+2,m+3 at right;
  47 # list [...] subscripts run from 1 to op(list).
  48 
  49 # Initialiser: Input list of m > 0 rank values (non-decreasing order);
  50 # adjacency improves for increasing partition;
  51 # output R_, P_, Q_, D_ 1-dim arrays of integer.
  52 #initbagperm := proc (rank :: list(integer), R_, P_, Q_, D_) :: NULL;
  53 # local m :: integer, j :: integer, r :: array(integer),
  54 # r_1 :: integer, r_m :: integer;
  55 initbagperm := proc (rank, R_, P_, Q_, D_)
  56 local m, j, r, r_1, r_m;
  57 m := nops(rank); r := sort(rank);
  58 if m <= 0 then print(`bagperm: length <= 0`)
  59 else r_1, r_m := r[1], r[m];
  60 R_ := array(0..m+3, [r_1-1, seq(r[j], j = 1..m), r_1-1, r_m+1, r_m+2]);
  61 P_ := array(0..m+3, [0, seq(j, j = 1..m), 0, m+1, m+2]); # permed identities
  62 Q_ := array(0..m+3, [0, seq(j, j = 1..m), m+2, m+3, 0]); # inverse locations
  63 D_ := array(0..m+3, [0, seq(+1, j = 1..m), 0, +1, +1]); # permed drifts
  64 fi; NULL end;
  65 
  66 # Iterator: Output permuted ranks in R_[1..m];
  67 # result = location j of run drifted; on termination j = m+1.
  68 #nextbagperm := proc (R_ :: array(integer), P_ :: array(integer),
  69 # Q_ :: array(integer), D_ :: array(integer)) :: integer;
  70 # local i :: integer, j :: integer, k :: integer, l :: integer,
  71 # d :: integer, e :: integer, p :: integer;
  72 nextbagperm := proc (R_, P_, Q_, D_)
  73 local i, j, k, l, d, e, p;
  74 # locate earliest unblocked element at j, starting at blocked element 0
  75 j, i, d := 0, 0, 0;
  76 while R_[j] >= R_[i] do
  77 D_[j] := -d; # blocked at j; reverse drift d pre-emptively
  78 j := Q_[P_[j]+1]; d := D_[j]; i := j+d; # next element at j, neighbour at i
  79 if R_[j-1] <> R_[j] then l := j # save left end of run in l
  80 elif d < 0 then i := l-1 fi od; # restore left end at head of run
  81 # shift run of equal rank from i-d,i-2d,...,l to i,i-d,...,l+d
  82 k := i; if d < 0 then l := j fi;
  83 e := D_[i]; p := P_[i]; # save neighbour drift e and identifier p
  84 while k <> l do
  85 P_[k] := P_[k-d]; Q_[P_[k]] := k;
  86 D_[k] := -1; k := k-d od; # reset drifts of run tail elements
  87 R_[l], R_[i] := R_[i], R_[l]; # transpose user ranks
  88 D_[l], D_[i] := e, d; # restore drifts of head and neighbour
  89 P_[l], Q_[p] := p, l; # wrap neighbour around to other end
  90 j end;
  91 
  92 # Test bag-perms --- listing 
  93 parts := [2, 2, 2]; # chosen partition of set 
  94 n := nops(parts): m := add(parts[i], i = 1..n): # length and weight 
  95 symb := [seq(seq(i, j = 1..parts[i]), i = 1..n)]: 
  96 cardi := combinat[multinomial](m, op(parts)): # number of bag-perms 
  97 initbagperm(symb, R_,P_,Q_,D_): 
  98 n := 1: print(n, [seq(R_[i], i = 1..m)]); 
  99 while nextbagperm(R_,P_,Q_,D_) <= m and n < cardi + 2 do 
 100   n := n+1; print(n, [seq(R_[i], i = 1..m)]) od: 
 101 print(n, cardi); # equal 
 102 
 103 ################################################################################