1 ;;;; COMBINATORICS package for Maxima
3 ;;;; Maxima functions to work with permutations.
4 ;;;; $cyclep : predicate function for cycles.
5 ;;;; $permp : predicate function for permutations.
6 ;;;; $permult : permutations product.
7 ;;;; $perm_inverse : inverse of a permutation.
8 ;;;; $permute : permutes a list according to a permutation.
9 ;;;; $apply_cycles : permutes a list according to a list of cycles.
10 ;;;; $perm_undecomp : converts a list of cycles into a permutation.
11 ;;;; $perm_cycles : decomposes a permutation into cycles.
12 ;;;; $perm_lex_rank : permutation's position in lexicographic order.
13 ;;;; $perm_lex_unrank : permutation's position in lexicographic order.
14 ;;;; $perm_length : minimum number of adjacent transpositions.
15 ;;;; $perm_decomp : minimum set of adjacent transpositions to get p.
16 ;;;; $perm_parity : parity of a permutation (0 or 1).
17 ;;;; $perms_lex : lexicographic list of all n-degree permutations
18 ;;;; or permutations within some ranks.
19 ;;;; $perm_lex_next : finds next permutation in lexicographic order.
20 ;;;; $perm_rank : permutation's position in Trotter-Johnson order.
21 ;;;; $perm_unrank : permutation's position in Trotter-Johnson order.
22 ;;;; $perms : n-degree permutations in mimimum-change order
23 ;;;; (Trotter-Johnson) or permutations within some ranks.
24 ;;;; $perm_next : finds next permutation in Trotter-Johnson order.
25 ;;;; $random_perm : generates a random permutation of degree n.
28 ;;;; pm : a permutation as a maxima list
29 ;;;; pa : a permutation as a lisp array
30 ;;;; lpm : a lisp list of pm's (in reverse order)
31 ;;;; lmpm : a maxima list of pm's
32 ;;;; cm : a cycle as a maxima list
33 ;;;; ca : a cycle as a lisp array
34 ;;;; lcm : a lisp list of cycles as maxima lists
35 ;;;; lmcm : a maxima list of cycles as maxima lists
38 ;;;; Turn a pa into a pm:
39 ;;;; (concatenate 'list `((mlist simp)) pa)
41 ;;;; Turn an lpm into a maxima list
42 ;;;; `((mlist simp) ,@(nreverse lpm))
45 ;;;; Apply transposition [i,j] to a pa (i and j from 1 to n)
46 ;;;; (array-transposition pa i j)
48 ;;;; Apply transposition [i,i+1] to a pa (i from 1 to n-1)
49 ;;;; (array-adjacent-transposition pa i)
51 ;;;; Turn a pa into a pm and push it to the top of an lpm
52 ;;;; (push-array-mlist pa lpm)
54 ;;;; Copyright (C) 2017 Jaime Villate
56 ;;;; This program is free software; you can redistribute it and/or
57 ;;;; modify it under the terms of the GNU General Public License
58 ;;;; as published by the Free Software Foundation; either version 2
59 ;;;; of the License, or (at your option) any later version.
61 ;;;; This program is distributed in the hope that it will be useful,
62 ;;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
63 ;;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
64 ;;;; GNU General Public License for more details.
66 ;;;; You should have received a copy of the GNU General Public License
67 ;;;; along with this program; if not, write to the Free Software
68 ;;;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
69 ;;;; MA 02110-1301, USA
72 (macsyma-module combinatorics
)
74 (defmacro push-array-mlist
(pa lpm
)
75 `(push (concatenate (quote list
) (quote ((mlist simp
))) ,pa
) ,lpm
))
77 (defmacro array-adjacent-transposition
(pa i
)
78 `(rotatef (aref ,pa
(1- ,i
)) (aref ,pa
,i
)))
80 (defmacro array-transposition
(pa i j
)
81 `(rotatef (aref ,pa
(1- ,i
)) (aref ,pa
(1- ,j
))))
83 ;;; $cyclep returns true if its argument is a maxima list of length n or
84 ;;; lower, whose elements are integers between 1 and n, without repetitions.
86 ;; a cycle must be a maxima list
87 (or (and (integerp n
) (> n
1)) (return-from $cyclep nil
))
88 (or ($listp cm
) (return-from $cyclep nil
))
90 ;; cycle elements must be positive integers, less or equal to n
91 (or (and (integerp i
) (plusp i
) (<= i n
)) (return-from $cyclep nil
))
92 ;; and cannot be repeated
93 (or (< (count i
(rest cm
)) 2) (return-from $cyclep nil
)))
96 ;;; $permp returns true if its argument is a maxima list of length n,
97 ;;; whose elements are the integers 1, 2, ...n, without repetitions.
99 (or ($listp pm
) (return-from $permp nil
))
100 (let ((n ($length pm
)))
101 ;; a permutation of degree n is also a valid cycle of degree n
104 (defun check-perm (pm caller
)
105 (let ((s "~M: argument should be a permutation; found: ~M"))
106 (or ($permp pm
) (merror (intl:gettext s
) caller pm
))))
108 (defun check-cycle (cm n caller
)
109 (let ((s "~M: argument should be a cycle of degree ~M; found: ~M"))
110 (or ($cyclep cm n
) (merror (intl:gettext s
) caller n cm
))))
112 (defun check-list-or-set (lsm caller
)
113 (let ((s "~M: argument should be a list or set; found: ~M"))
114 (or ($listp lsm
) ($setp lsm
) (merror (intl:gettext s
) caller lsm
))))
116 (defun check-list (lm caller
)
117 (let ((s "~M: argument should be a list; found: ~M"))
118 (or ($listp lm
) (merror (intl:gettext s
) caller lm
))))
120 (defun check-list-length (lm n caller
)
121 (let ((s "~M: argument should be a list of length ~M; found: ~M"))
122 (or (eql ($length lm
) n
) (merror (intl:gettext s
) caller n lm
))))
124 (defun check-pos-integer (n caller
)
125 (let ((s "~M: argument must be a positive integer; found: ~M"))
126 (or (and (integerp n
) (> n
0)) (merror (intl:gettext s
) caller n
))))
128 (defun check-integer-n1-n2 (r n1 n2 caller
)
129 (let ((s "~M: argument must be an integer between ~M and ~M; found: ~M"))
130 (or (and (integerp r
) (>= r n1
) (<= r n2
))
131 (merror (intl:gettext s
) caller n1 n2 r
))))
133 ;;; $permult multiplies permutation pm by the permutations in list lmpm
134 (defun $permult
(pm &rest lmpm
)
135 (check-perm pm
"permult")
136 (let ((n ($length pm
)) (rm (copy-list pm
)))
137 ;; multiplies the current product, rm, by each of the permutations in lmpm
139 (check-list-length qm n
"permult")
140 (check-perm qm
"permult")
141 ;; the new product, rm, is the current product, rm, times
142 ;; permutation qm: rm = [qm[rm[1]], qm[rm[2]], ..., qm[rm[n]]]
143 (loop for j from
1 to n do
144 (setf (nth j rm
) (nth (nth j rm
) qm
))))
147 ;;; $perm_inverse computes the inverse of permutation p
148 (defun $perm_inverse
(pm)
149 (check-perm pm
"perm_inverse")
150 (let ((inverse (copy-list pm
)))
151 ;; if the element pm[j] of the permutation is equal to i, then
152 ;; the element inverse[i] of the inverse should be equal to j.
153 (dotimes (j ($length pm
))
154 (setf (nth (nth (1+ j
) pm
) inverse
) (1+ j
)))
157 ;;; $permute permutes the elements of list (or set) lsm according
158 ;;; to the permutation pm
159 (defun $permute
(pm lsm
)
160 (check-perm pm
"permute")
161 ; (check-list-or-set lsm "permute")
162 (check-list-length lsm
($length pm
) "permute")
164 (dolist (j (reverse (rest pm
)))
165 (setq result
(cons (nth j lsm
) result
)))
166 (cons (car lsm
) result
)))
168 ;;; $apply_cycles permutes the elements of list (or set) lsm according to the
169 ;;; list of cycles lmcm, applying the cycles on the end of the list first
170 (defun $apply_cycles
(lmcm lsm
)
171 (check-list lmcm
"apply_cycles")
172 (check-list-or-set lsm
"apply_cycles")
173 (let ((n ($length lsm
)) (r (copy-list lsm
)) m tmp
)
174 ;; iterate the cycles, the last in the list first
175 (dolist (cm (reverse (rest lmcm
)))
176 (check-cycle cm n
"apply_cycles")
177 (setq m
($length cm
))
178 (setq tmp
(nth (nth 1 cm
) r
))
179 (loop for i from
1 to
(1- m
) do
180 (setf (nth (nth i cm
) r
) (nth (nth (1+ i
) cm
) r
)))
181 (setf (nth (nth m cm
) r
) tmp
))
184 ;;; $perm_undecomp converts a list of cycles into a permutation,
185 ;;; equal to their product, by applying $apply_cycles to [1, 2,..., n]
186 (defun $perm_undecomp
(lmcm n
)
187 (check-list lmcm
"perm_undecomp")
188 (check-pos-integer n
"perm_undecomp")
190 lmcm
`((mlist simp
) ,@(loop for i from
1 to n collecting i
))))
192 ;;; $perm_cycles decomposes permutation p into a product of canonical
193 ;;; cycles: with lower indices first.
194 (defun $perm_cycles
(pm)
195 (check-perm pm
"perm_cycles")
196 (let* (i j k cm
(n ($length pm
)) (result '((mlist simp
)))
197 (v (make-array n
:element-type
'bit
:initial-element
1)))
198 ;; the i'th bit in bit array v equal to 1 means index i+1 has not been
199 ;; added to any cycles. The next cycle will then start where the first
200 ;; value of 1 is found in v
201 (while (position 1 v
)
202 (setf i
(position 1 v
) (bit v i
) 0)
203 ;; i= v index where a 1 was found. j=i+1 first index in the current
204 ;; cycle. k=next index in the current cycle
205 (setq j
(1+ i
) k j cm
`((mlist simp
) ,j
))
206 ;; if p[k] is different from k, there's a non-trivial cycle
207 (while (/= (nth k pm
) j
)
208 (setf k
(nth k pm
) cm
(append cm
(list k
)) (bit v
(1- k
)) 0))
209 ;; a trivial cycle with just one index (length 2) will not be saved
210 (and (> (length cm
) 2) (setq result
(append result
(list cm
)))))
213 ;;; $perm_lex_rank finds the position of the given permutation in
214 ;;; the lexicographic ordering of permutations (from 1 to n!).
215 ;;; Algorithm 2.15: from Kreher & Stinson (1999). Combinatorial Algorithms.
216 (defun $perm_lex_rank
(pm)
217 (check-perm pm
"perm_lex_rank")
218 (let ((r 1) (n ($length pm
)) (qm (copy-list pm
)))
219 (loop for j from
1 to n do
220 (incf r
(* (1- (nth j qm
)) (factorial (- n j
))))
221 (loop for i from
(1+ j
) to n do
222 (when (> (nth i qm
) (nth j qm
))
223 (setf (nth i qm
) (1- (nth i qm
))))))
226 ;;; $perm_length finds the minimum number of adjacent transpositions
227 ;;; necessary to write permutation p as a product of adjacent transpositions.
228 (defun $perm_length
(pm)
229 (check-perm pm
"perm_length")
230 (let ((d 0) (n ($length pm
)) (qm (copy-list pm
)))
231 (loop for j from
1 to n do
232 (incf d
(* (1- (nth j qm
))))
233 (loop for i from
(1+ j
) to n do
234 (when (> (nth i qm
) (nth j qm
))
235 (setf (nth i qm
) (1- (nth i qm
))))))
238 ;;; $perm_decomp finds the minimum set of adjacent transpositions
239 ;;; whose product equals permutation pm.
240 (defun $perm_decomp
(pm)
241 (check-perm pm
"perm_decomp")
242 (let (lcm (n ($length pm
)) (qm (copy-list pm
)))
243 (loop for j from
1 to
(1- n
) do
244 (when (>= (1- (nth j qm
)) 1)
245 (loop for i from
(+ (- (nth j qm
) 2) j
) downto j do
246 (setq lcm
(cons `((mlist simp
) ,i
,(1+ i
)) lcm
))))
247 (loop for k from
(1+ j
) to n do
248 (when (> (nth k qm
) (nth j qm
))
249 (setf (nth k qm
) (1- (nth k qm
))))))
250 (cons '(mlist simp
) lcm
)))
252 ;;; $perm_parity finds the parity of permutation p (0 or 1).
253 (defun $perm_parity
(pm)
254 (check-perm pm
"perm_parity")
255 (perm-parity ($length pm
) (rest pm
)))
257 ;;; perm-parity finds the parity of n first elements of p (lisp list)
258 ;;; Algorithm 2.19 from Kreher & Stinson (1999). Combinatorial Algorithms.
259 (defun perm-parity (n p
)
260 (let (i (c 0) (a (make-array n
:element-type
'bit
:initial-element
0)))
262 (when (= (bit a j
) 0)
264 (setf (bit a j
) 1 i j
)
265 (while (/= (1- (nth i p
)) j
)
266 (setf i
(1- (nth i p
)) (bit a i
) 1))))
269 ;;; $perms_lex returns a list of permutations of degree n in lexicographic order
270 (defun $perms_lex
(n &optional r0 rf
)
271 (check-pos-integer n
"perms_lex")
272 (when r0
(check-integer-n1-n2 r0
1 (factorial n
) "perms_lex"))
273 (when rf
(check-integer-n1-n2 rf r0
(factorial n
) "perms_lex"))
274 (let ((pa (make-array n
:element-type
'fixnum
)) lpm
)
277 (perm-lex-unrank n r0 pa
)
278 (push-array-mlist pa lpm
)
280 (dotimes (j (- rf r0
))
282 (push-array-mlist pa lpm
))))
284 (dotimes (i n
) (setf (aref pa i
) (1+ i
)))
285 (push-array-mlist pa lpm
)
286 (while (perm-lex-next n pa
)
287 (push-array-mlist pa lpm
))))
288 `((mlist simp
) ,@(nreverse lpm
))))
290 ;;; $perm_lex_next finds the next permutation in lexicographic order
291 (defun $perm_lex_next
(pm)
292 (check-perm pm
"perm_lex_next")
293 (let* ((n ($length pm
)) (pa (make-array n
:element-type
'fixnum
)))
295 (setf (aref pa i
) (nth (1+ i
) pm
)))
296 (unless (perm-lex-next n pa
)
297 (return-from $perm_lex_next nil
))
298 (concatenate 'list
'((mlist simp
)) pa
)))
300 ;;; $perm_lex_unrank finds the n-degree permutation with lexicographic rank r
301 (defun $perm_lex_unrank
(n r
)
302 (check-pos-integer n
"perm_lex_unrank")
303 (check-integer-n1-n2 r
1 (factorial n
) "perm_lex_unrank")
304 (let ((pa (make-array n
:element-type
'fixnum
)))
305 (perm-lex-unrank n r pa
)
306 (concatenate 'list
`((mlist simp
)) pa
)))
308 ;;; perm-lex-next finds next permutation in the lexicographic order.
309 ;;; Based on Algorithm 2.14 from Kreher & Stinson (1999). Combinatorial
311 (defun perm-lex-next (n pa
)
312 (declare (type (simple-array fixnum
*) pa
))
313 (declare (type fixnum n
))
314 (let ((i (1- n
)) (j n
))
315 (while (and (> i
0) (< (aref pa i
) (aref pa
(1- i
))))
317 (when (= i
0) (return-from perm-lex-next nil
))
318 (while (< (aref pa
(1- j
)) (aref pa
(1- i
)))
320 (array-transposition pa i j
)
321 (dotimes (k (floor (/ (- n i
) 2)))
322 (array-transposition pa
(+ k i
1) (- n k
)))
325 ;;; perm-lex-unrank finds the n-degree permutation of in position
326 ;;; r (from 1 to n!) in the lexicographic ordering of permutations.
327 ;;; Algorithm 2.16: from Kreher & Stinson (1999). Combinatorial Algorithms.
328 (defun perm-lex-unrank (n r pa
)
329 (declare (type (simple-array fixnum
*) pa
))
330 (declare (type fixnum n
))
331 (declare (type fixnum r
))
333 (setf (aref pa
(1- n
)) 1)
335 (setq d
(/ (mod (1- r
) (factorial (+ j
2))) (factorial (1+ j
))))
336 (decf r
(* d
(factorial (1+ j
))))
337 (setf (aref pa
(- n j
2)) (1+ d
))
338 (loop for i from
(- n j
1) to
(1- n
) do
339 (when (> (aref pa i
) d
)
340 (setf (aref pa i
) (1+ (aref pa i
))))))
343 ;;; $perms returns a list of the permutations of order n
344 ;;; in Trotter-Johnson ordering, using a simpler algorithm
345 (defun $perms
(n &optional r0 rf
)
346 (check-pos-integer n
"perms")
347 (when r0
(check-integer-n1-n2 r0
1 (factorial n
) "perms"))
348 (when rf
(check-integer-n1-n2 rf r0
(factorial n
) "perms"))
349 (let ((pa (make-array n
:element-type
'fixnum
)) lpm
)
352 (perm-unrank n r0 pa
)
353 (push-array-mlist pa lpm
)
355 (loop for r from r0 to
(- rf
1) do
356 (array-adjacent-transposition pa
(transposition-next n r
))
357 (push-array-mlist pa lpm
))))
359 (dotimes (i n
) (setf (aref pa i
) (1+ i
)))
360 (push-array-mlist pa lpm
)
361 (loop for r from
1 to
(1- (factorial n
)) do
362 (array-adjacent-transposition pa
(transposition-next n r
))
363 (push-array-mlist pa lpm
))))
364 `((mlist simp
) ,@(nreverse lpm
))))
366 ;;; $perm_next finds the next permutation in Trotter-Johnson ordering
367 (defun $perm_next
(pm)
368 (check-perm pm
"perm_next")
369 (let ((n ($length pm
)) (r ($perm_rank pm
)) (qm (copy-list pm
)) i
)
370 (when (= r
(factorial n
)) (return-from $perm_next nil
))
371 (setq i
(transposition-next n r
))
372 (rotatef (nth i qm
) (nth (1+ i
) qm
))
375 ;;; $perm_unrank finds the n-degree permutation with Trotter-Johnson rank r
376 (defun $perm_unrank
(n r
)
377 (check-pos-integer n
"perm_unrank")
378 (check-integer-n1-n2 r
1 (factorial n
) "perm_unrank")
379 (let ((pa (make-array n
:element-type
'fixnum
)))
381 (concatenate 'list
`((mlist simp
)) pa
)))
383 ;;; $perm_rank finds the position of the given permutation in
384 ;;; the Trotter-Johnson ordering of permutations (from 1 to n!).
385 ;;; Algorithm 2.17: from Kreher & Stinson (1999). Combinatorial Algorithms.
386 (defun $perm_rank
(pm)
387 (check-perm pm
"perm_rank")
388 (let (i k
(r 1) (n ($length pm
)))
389 (loop for j from
2 to n do
391 (while (/= (nth i pm
) j
)
392 (and (< (nth i pm
) j
) (incf k
))
395 (setq r
(- (+ (* j r
) 1) k
))
396 (setq r
(+ (* j
(1- r
)) k
))))
399 ;;; perm-unrank returns the n-degree permutation in position r
400 ;;; (from 1 to n!) in the Trotter-Johnson ordering of permutations.
401 ;;; Algorithm 2.18: from Kreher & Stinson (1999). Combinatorial Algorithms.
402 (defun perm-unrank (n r pa
)
403 (declare (type (simple-array fixnum
*) pa
))
404 (declare (type fixnum n
))
405 (declare (type fixnum r
))
408 (loop for j from
2 to n do
409 (setq r1
(floor (/ (* (1- r
) (factorial j
)) (factorial n
))))
410 (setq k
(- r1
(* j r2
)))
413 (loop for i from
(1- j
) downto
(- j k
) do
414 (setf (aref pa i
) (aref pa
(1- i
))))
415 (setf (aref pa
(- j k
1)) j
))
417 (loop for i from
(1- j
) downto
(1+ k
) do
418 (setf (aref pa i
) (aref pa
(1- i
))))
419 (setf (aref pa k
) j
)))
422 ;;; transposition-next finds the first index of the adjacent transposition
423 ;;; to be applied to the n-degree permutation of Trotter-Johnson rank r, in
424 ;;; order to get the permutation of rank r+1 (r must be between 1 and n!-1)
425 (defun transposition-next (n r
)
427 (setf (values m k
) (floor r n
))
431 (when (= (mod m
2) 1) (incf d
))
432 (setf (values m k
) (floor r n
)))
437 ;;; $random_perm generates a random permutation of degree n, using algorithm
438 ;;; 5.4 from Reingold, Nievergelt and Deo. Combinatorial algorithms: theory and
440 (defun $random_perm
(n)
441 (check-pos-integer n
"random_perm")
442 (let ((pm (cons '(mlist) (loop for i from
1 to n collecting i
))))
443 (loop for i from n downto
2 do
444 (rotatef (nth i pm
) (nth (1+ ($random i
)) pm
)))