chiark / gitweb /
progs/perftest.c: Use from Glibc syscall numbers.
[catacomb] / utils / permute.lisp
1 ;;; -*-lisp-*-
2
3 ;;; This file isn't a program as such: rather, it's a collection of handy
4 ;;; functions which can be used in an interactive session.
5
6 ;;;--------------------------------------------------------------------------
7 ;;; General permutation utilities.
8
9 (defun shuffle (v)
10   "Randomly permute the elements of the vector V.  Return V."
11   (let ((n (length v)))
12     (do ((k n (1- k)))
13         ((<= k 1) v)
14       (let ((i (random k)))
15         (unless (= i (1- k))
16           (rotatef (aref v i) (aref v (1- k))))))))
17
18 (defun identity-permutation (n)
19   "Return the do-nothing permutation on N elements."
20   (let ((v (make-array n :element-type 'fixnum)))
21     (dotimes (i n v) (setf (aref v i) i))))
22
23 (defun invert-permutation (p)
24   "Given a permutation P, return its inverse."
25   (let* ((n (length p)) (p-inv (make-array n :element-type 'fixnum)))
26     (dotimes (i n) (setf (aref p-inv (aref p i)) i))
27     p-inv))
28
29 (defun next-permutation (v)
30   "Adjust V so that it reflects the next permutation in ascending order.
31
32    V should be a vector of real numbers.  Returns V if successful, or nil if
33    there are no more permutations."
34
35   ;; The tail of the vector consists of a sequence ... A, Z, Y, X, ..., where
36   ;; Z > Y > X ... is in reverse order, and A < Z.  The next permutation is
37   ;; then the smallest out of Z, Y, X, ... which is larger than A, followed
38   ;; by the remaining elements in ascending order.
39   ;;
40   ;; Equivalently, reverse the tail Z, Y, X, ... so we have A, ... X, Y, Z,
41   ;; and swap A with the next larger element.
42
43   (let ((n (length v)))
44     (cond ((< n 2) nil)
45           (t (let* ((k (1- n))
46                     (x (aref v k)))
47                (loop (when (zerop k) (return-from next-permutation nil))
48                      (decf k)
49                      (let ((y (aref v k)))
50                        (when (prog1 (< y x)
51                                (setf x y))
52                          (return))))
53                (do ((i (1+ k) (1+ i))
54                     (j (1- n) (1- j)))
55                    ((> i j))
56                  (rotatef (aref v i) (aref v j)))
57                (do ((i (- n 2) (1- i)))
58                    ((or (<= i k) (< (aref v i) x))
59                     (rotatef (aref v k) (aref v (1+ i)))))
60                v)))))
61
62 (defun make-index-mask (w mask-expr)
63   "Construct a bitmask based on bitwise properties of the bit indices.
64
65    The function returns a W-bit mask in which each bit is set if MASK-EXPR
66    of true of the bit's index.  MASK-EXPR may be one of the following:
67
68      * I -- an integer I is true if bit I of the bit index is set;
69      * (not EXPR) -- is true if EXPR is false;
70      * (and EXPR EXPR ...) -- is true if all of the EXPRs are true; and
71      * (or EXPR EXPR ...) -- is true if any of the EXPRs is true."
72
73   (let ((max-bit (1- (integer-length (1- w))))
74         (mask 0))
75     (dotimes (i w mask)
76       (labels ((interpret (expr)
77                  (cond ((and (integerp expr) (<= 0 expr max-bit))
78                         (logbitp expr i))
79                        ((and (consp expr) (eq (car expr) 'not)
80                              (null (cddr expr)))
81                         (not (interpret (cadr expr))))
82                        ((and (consp expr) (eq (car expr) 'and))
83                         (every #'interpret (cdr expr)))
84                        ((and (consp expr) (eq (car expr) 'or))
85                         (some #'interpret (cdr expr)))
86                        (t
87                         (error "unknown mask expression ~S" expr)))))
88         (when (interpret mask-expr)
89           (setf (ldb (byte 1 i) mask) 1))))))
90
91 (defun make-permutation-network (w steps)
92   "Construct a permutation network.
93
94    The integer W gives the number of bits to be acted upon.  The STEPS are a
95    list of instructions of the following forms:
96
97      * (SHIFT . MASK) -- a pair of integers is treated literally;
98
99      * (SHIFT MASK-EXPR) -- the SHIFT is literal, but the MASK-EXPR is
100        processed by `make-index-mask' to calculate the mask;
101
102      * (:invert I) -- make an instruction which inverts the sense of the
103        index bit I;
104
105      * (:exchange I J) -- make an instruction which exchanges index bits I
106        and J; or
107
108      * (:exchange-invert I J) -- make an instruction which exchanges and
109        inverts index bits I and J.
110
111    The output is a list of primitive (SHIFT . MASK) steps, indicating that
112    the bits of the input selected by MASK are to be swapped with the bits
113    selected by (ash MASK SHIFT)."
114
115   (let ((max-mask (1- (ash 1 w)))
116         (max-shift (1- w))
117         (max-bit (1- (integer-length (1- w))))
118         (list nil))
119     (dolist (step steps)
120       (cond ((and (consp step)
121                   (integerp (car step)) (<= 0 (car step) max-shift)
122                   (integerp (cdr step)) (<= 0 (cdr step) max-mask))
123              (push step list))
124             ((and (consp step)
125                   (integerp (car step)) (<= 0 (car step) max-shift)
126                   (null (cddr step)))
127              (push (cons (car step) (make-index-mask w (cadr step))) list))
128             ((and (consp step)
129                   (eq (car step) :invert)
130                   (integerp (cadr step)) (<= 0 (cadr step) max-bit)
131                   (null (cddr step)))
132              (let ((i (cadr step)))
133                (push (cons (ash 1 i) (make-index-mask w `(not ,i))) list)))
134             ((and (consp step)
135                   (eq (car step) :exchange)
136                   (integerp (cadr step)) (integerp (caddr step))
137                   (<= 0 (cadr step) (caddr step) max-bit)
138                   (null (cdddr step)))
139              (let ((i (cadr step)) (j (caddr step)))
140                (push (cons (- (ash 1 j) (ash 1 i))
141                            (make-index-mask w `(and ,i (not ,j))))
142                      list)))
143             ((and (consp step)
144                   (eq (car step) :exchange-invert)
145                   (integerp (cadr step)) (integerp (caddr step))
146                   (<= 0 (cadr step) (caddr step) max-bit)
147                   (null (cdddr step)))
148              (let ((i (cadr step)) (j (caddr step)))
149                (push (cons (+ (ash 1 i) (ash 1 j))
150                            (make-index-mask w `(and (not ,i) (not ,j))))
151                      list)))
152             (t
153              (error "unknown permutation step ~S" step))))
154     (nreverse list)))
155
156 ;;;--------------------------------------------------------------------------
157 ;;; Permutation network diagnostics.
158
159 (defun print-permutation-network (steps &optional (stream *standard-output*))
160   "Print a description of the permutation network STEPS to STREAM.
161
162    A permutation network consists of a list of pairs
163
164         (SHIFT . MASK)
165
166    indicating that the bits selected by MASK, and those SHIFT bits to the
167    left, should be exchanged.
168
169    The output is intended to be human-readable and is subject to change."
170
171   (let ((shiftwd 1) (maskwd 2))
172
173     ;; Determine suitable print widths for shifts and masks.
174     (dolist (step steps)
175       (let ((shift (car step)) (mask (cdr step)))
176         (let ((swd (1+ (floor (log shift 10))))
177               (mwd (ash 1 (- (integer-length (1- (integer-length mask)))
178                              2))))
179           (when (> swd shiftwd) (setf shiftwd swd))
180           (when (> mwd maskwd) (setf maskwd mwd)))))
181
182     ;; Print the display.
183     (pprint-logical-block (stream steps :prefix "(" :suffix ")")
184       (let ((first t))
185         (dolist (step steps)
186           (let ((shift (car step)) (mask (cdr step)))
187
188             ;; Separate entries with newlines.
189             (cond (first (setf first nil))
190                   (t (pprint-newline :mandatory stream)))
191
192             (let ((swaps nil))
193
194               ;; Determine the list of exchanges implied by the mask.
195               (dotimes (i (integer-length mask))
196                 (when (logbitp i mask)
197                   (push (cons i (+ i shift)) swaps)))
198               (setf swaps (nreverse swaps))
199
200               ;; Print the entry.
201               (format stream "~@<(~;~vD #x~(~v,'0X~) ~8I~:@_~W~;)~:>"
202                       shiftwd shift maskwd mask swaps))))))
203
204     ;; Print a final newline following the close parenthesis.
205     (terpri stream)))
206
207 (defun demonstrate-permutation-network
208     (n steps
209      &key reference
210           offsets
211           (stream *standard-output*))
212   "Print, on STREAM, a demonstration of the permutation STEPS.
213
214    The output is intended to be useful to human readers and is subject to
215    change.  Currently, it prints a sequence of diagrams on STREAM.  The left
216    hand side of each row shows a map of which bits are affected: `-' means
217    that the bit remains in the same position, `*' means that it moves
218    forward, and `#' means that it moves back; each `*' pairs with the
219    earliest unpaired `#' marker.  The right hand side shows the arrangement
220    of the original input bits.
221
222    If OFFSETS is non-nil, then also print a table of offsets showing how far
223    each bit has yet to move.
224
225    If REFERENCE is not nil, then print a final pair of diagrams.  This time,
226    the map shows `-' for correct bits and `x' for incorrect ones, with the
227    right hand side showing the expected arrangement of input bits.
228
229    The function returns non-nil if the STEPS resulted in the REFERENCE
230    permutation, and nil if either the STEPS are incorrect or no REFERENCE was
231    provided."
232
233   (flet ((apply-step (shift mask v)
234            (dotimes (k n)
235              (when (logbitp k mask)
236                (rotatef (aref v k) (aref v (+ k shift)))))))
237
238     (let* ((v (identity-permutation n))
239            (end (or reference
240                     (let ((e (identity-permutation n)))
241                       (dolist (step steps e)
242                         (let ((shift (car step)) (mask (cdr step)))
243                           (apply-step shift mask e))))))
244            (end-inv (invert-permutation end))
245            (mapwd (ceiling (sqrt n)))
246            (ixwd (length (format nil "~D" (1- n)))))
247
248       (flet ((show-stage (shift mask v)
249                (do ((i 0 (+ i mapwd)))
250                    ((>= i n))
251                  (write-string ";;      " stream)
252                  (dotimes (j mapwd)
253                    (let ((k (+ i j)))
254                      (when (plusp j) (write-char #\space stream))
255                      (write-char (cond ((>= k n)
256                                         #\space)
257                                        ((logbitp k mask)
258                                         #\*)
259                                        ((and (>= k shift)
260                                              (logbitp (- k shift) mask))
261                                         #\#)
262                                        (t
263                                         #\-))
264                                  stream)))
265                  (write-string " | " stream)
266
267                  (when offsets
268                    (dotimes (j mapwd)
269                      (let ((k (+ i j)))
270                        (when (plusp j) (write-char #\space stream))
271                        (cond ((>= k n)
272                               (format stream "~v@T" (1+ ixwd)))
273                              (t
274                               (format stream "~*~[~2:*~vD~:;~2:*~v@D~]"
275                                       (1+ ixwd)
276                                       (- (aref end-inv (aref v k)) k))))))
277                    (write-string " | " stream))
278
279                  (dotimes (j (min mapwd (- n i)))
280                    (let ((k (+ i j)))
281                      (when (plusp j) (write-char #\space stream))
282                      (format stream "~vD" ixwd (aref v k))))
283                  (terpri))))
284
285         (fresh-line)
286         (show-stage 0 0 v)
287
288         (dolist (step steps)
289           (let ((shift (car step)) (mask (cdr step)))
290             (apply-step shift mask v)
291             (format stream ";;~%")
292             (show-stage shift mask v)))
293
294         (let ((ok (not (null reference))))
295           (when reference
296             (format stream ";;~%")
297             (do ((i 0 (+ i mapwd)))
298                 ((>= i n))
299               (write-string ";; " stream)
300               (dotimes (j mapwd)
301                 (let ((k (+ i j)))
302                   (when (plusp j) (write-char #\space stream))
303                   (write-char (cond ((>= k n) #\space)
304                                     ((/= (aref v k) (aref reference k)) #\x)
305                                     (t #\-))
306                               stream)))
307               (write-string " | " stream)
308               (dotimes (j (min mapwd (- n i)))
309                 (let ((k (+ i j)))
310                   (when (plusp j) (write-char #\space stream))
311                   (format stream "~vD" ixwd (aref reference k))))
312               (terpri))
313             (unless (every #'= v reference)
314               (setf ok nil))
315             (format stream "~:[FAIL~;pass~]~%" ok))
316           ok)))))
317
318 ;;;--------------------------------------------------------------------------
319 ;;; Beneš networks.
320
321 (defun compute-benes-step (n p p-inv bit clear-input)
322   "Compute a single layer of a Beneš network.
323
324    N is a fixnum.  P is a vector of fixnums defining a permutation: for each
325    output bit position i (numbering the least significant bit 0), element i
326    gives the number of the input which should end up in that position; P-INV
327    gives the inverse permutation in the same form.  BIT is a power of 2 which
328    gives the distance between bits we should consider.  CLEAR-INPUT is
329    a (generalized) boolean: if true, we attempt to do no work on the input
330    side; if false, we try to do no work on the output side.  The length of P
331    must be at least (logior N BIT).
332
333    The output consists of a pair of masks M0 and M1, to be used on the input
334    and output sides respectively.  The permutation stage, applied to an input
335    X, is as follows:
336
337         (let ((tmp (logand (logxor x (ash x (- bit))) mask)))
338           (logxor x tmp (ash tmp bit)))
339
340    The critical property of the masks is that it's possible to compute an
341    inner permutation, mapping the output of the M0 stage to the input of the
342    M1 stage, such that (a) the overall composition of the three permutations
343    is the given permutation P, and (b) the distances that the bits need to
344    be moved by the inner permutation all have BIT clear.
345
346    The resulting permutation will attempt to avoid touching elements at
347    indices greater than N.  This attempt will succeed if all such elements
348    contain values equal to their indices.
349
350    By appropriately composing layers computed by this function, then, it's
351    possible to perform an arbitrary permutation of 2^n bits in 2 n - 1 simple
352    steps like the one above."
353
354   ;; Consider the following problem.  You're given two equally-sized
355   ;; collections of things, called `left' and `right'.  Each left thing is
356   ;; attached to exactly one right thing with a string, and /vice versa/.
357   ;; Furthermore, the left things, and the right things, are each linked
358   ;; together in pairs, so each pair has two strings coming out of it.  Our
359   ;; job is to paint the strings so that each linked pair of things has one
360   ;; red string and one blue string.
361   ;;
362   ;; This is quite straightforward.  Pick a pair whose strings aren't yet
363   ;; coloured, and colour one of its strings, chosen arbitrarily, red.  Find
364   ;; the pair at the other end of this red string.  If the two other things
365   ;; in these two pairs are connected, then paint that string blue and move
366   ;; on.  Otherwise, both things have an uncoloured string, so paint both of
367   ;; them blue and trace along these now blue strings to find two more thing
368   ;; pairs.  Again, the other thing in each pair has an uncoloured string.
369   ;; If these things share the /same/ string, paint it red and move on.
370   ;; Otherwise, paint both strings red, trace, and repeat.  Eventually, we'll
371   ;; find two things joined by the same string, each paired with another
372   ;; thing whose strings we've just painted the same colour.  Once we're
373   ;; done, we'll have painted a cycle's worth of strings, and each pair of
374   ;; things will have either both of its strings painted different colours,
375   ;; or neither of them.
376   ;;
377   ;; The right things are the integers 0, 1, ..., n - 1, where n is some
378   ;; power of 2, paired according to whether they differ only in BIT.  The
379   ;; left things are the same integers, connected to the right things
380   ;; according to the permutation P: the right thing labelled i is connected
381   ;; to the left thing labelled P(i).  Similarly, two left things are paired
382   ;; if their labels P(i) and P(j) differ only in BIT.  We're going to paint
383   ;; a string red if we're going to arrange to clear BIT in the labels at
384   ;; both ends, possibly by swapping the two labels, and paint it red if
385   ;; we're going to arrange to set BIT.  Once we've done this, later stages
386   ;; of the filter will permute the red- and blue-painted things
387   ;; independently.
388
389   (let ((m0 0) (m1 0) (done 0))
390
391     ;; Now work through the permutation cycles.
392     (do ((i (1- n) (1- i)))
393         ((minusp i))
394
395       ;; Skip if we've done this one already.
396       (unless (or (plusp (logand i bit))
397                   (logbitp i done))
398
399         ;; Find the other associated values.
400         (let* ((i0 i) (i1 (aref p-inv i))
401                (sense (cond ((>= (logior i0 bit) n) 0)
402                             (clear-input (logand i0 bit))
403                             (t (logand i1 bit)))))
404
405           #+noise
406           (format t ";; new cycle: i0 = ~D, j0 = ~D; i1 = ~D, j1 = ~D~%"
407                   i0 (logxor i0 bit)
408                   i1 (logxor i1 bit))
409
410           ;; Mark this index as done.
411           (setf (ldb (byte 1 i0) done) 1)
412           #+noise (format t ";;   done = #x~2,'0X~%" done)
413
414           ;; Now trace round the cycle.
415           (loop
416
417             ;; Mark this index as done.
418             (setf (ldb (byte 1 (logandc2 i0 bit)) done) 1)
419             #+noise (format t ";;   done = #x~2,'0X~%" done)
420
421             ;; Swap the input and output pairs if necessary.
422             (unless (= (logand i0 bit) sense)
423               #+noise
424               (format t ";;   swap input:  ~D <-> ~D~%"
425                       (logandc2 i0 bit) (logior i0 bit))
426               (setf (ldb (byte 1 (logandc2 i0 bit)) m0) 1))
427             (unless (= (logand i1 bit) sense)
428               #+noise
429               (format t ";;   swap output: ~D <-> ~D~%"
430                       (logandc2 i1 bit) (logior i1 bit))
431               (setf (ldb (byte 1 (logandc2 i1 bit)) m1) 1))
432
433             ;; Advance around the cycle.
434             (let* ((j0 (logxor i0 bit))
435                    (j1 (logxor i1 bit))
436                    (next-i1 (aref p-inv j0))
437                    (next-i0 (aref p j1)))
438               (when (= next-i0 j0) (return))
439               (setf i0 next-i0
440                     i1 next-i1
441                     sense (logxor sense bit)))
442
443             #+noise
444             (format t ";; advance: i0 = ~D, j0 = ~D; i1 = ~D, j1 = ~D~%"
445                     i0 (logxor i0 bit)
446                     i1 (logxor i1 bit))))))
447
448     (values m0 m1)))
449
450 (defun compute-final-benes-step (n p p-inv bit)
451   "Determine the innermost stage of a Beneš network.
452
453    N is a fixnum.  P is a vector of fixnums defining a permutation: for each
454    output bit position i (numbering the least significant bit 0), element i
455    gives the number of the input which should end up in that position; P-INV
456    gives the inverse permutation in the same form.  BIT is a power of 2 which
457    gives the distance between bits we should consider.  The length of P must
458    be at least (logior N BIT).
459
460    Furthermore, the ith element of P must be equal either to i or to
461    (logxor i BIT); and therefore P-INV must be equal to P.
462
463    Return the mask such that
464
465         (let ((tmp (logand (logxor x (ash x (- bit))) mask)))
466           (logxor x tmp (ash tmp bit)))
467
468    applies the permutation P to the bits of x."
469
470   (declare (ignorable p-inv))
471
472   (let ((m 0))
473     (dotimes (i n)
474       (unless (plusp (logand i bit))
475         (let ((x (aref p i)))
476           #+paranoid
477           (assert (= (logandc2 x bit) i))
478           #+paranoid
479           (assert (= x (aref p-inv i)))
480
481           (unless (= x i)
482             (setf (ldb (byte 1 i) m) 1)))))
483     m))
484
485 (defun apply-benes-step (p p-inv bit m0 m1)
486   "Apply input and output steps for a Beneš network to a permutation.
487
488    Given the permutation P and its inverse, and the distance BIT, as passed
489    to `compute-benes-step', and the masks M0 and M1 returned, determine and
490    return the necessary `inner' permutation to be applied between these
491    steps, and its inverse.
492
493    A permutation-network step, and, in particular, a Beneš step, is an
494    involution, so the change to the vectors P and P-INV can be undone by
495    calling the function again with the same arguments."
496
497   (flet ((swaps (p p-inv mask)
498            (dotimes (i0 (length p))
499              (when (logbitp i0 mask)
500                (let* ((j0 (logior i0 bit))
501                       (i1 (aref p-inv i0))
502                       (j1 (aref p-inv j0)))
503                  (setf (aref p i1) j0
504                        (aref p j1) i0)
505                  (rotatef (aref p-inv i0) (aref p-inv j0)))))))
506     (swaps p p-inv m0)
507     (swaps p-inv p m1)
508
509     #+paranoid
510     (let* ((n (length p)))
511       (dotimes (i n)
512         (assert (= (aref p (aref p-inv i)) i))
513         (assert (= (aref p-inv (aref p i)) i))))))
514
515 (defun benes-search (p)
516   "Given a bit permutation P, describe a Beneš network implementing P.
517
518    P is a sequence of fixnums defining a permutation: for each output bit
519    position i (numbering the least significant bit 0), element i gives the
520    number of the input which should end up in that position.
521
522    The return value is a list of steps of the form
523
524         (BIT MASK (X . Y) (X' . Y') ...)
525
526    To implement this permutation step:
527
528      * given an input X, compute
529
530         (let ((tmp (logand (logxor x (ash x (- bit))) mask)))
531           (logxor x tmp (ash tmp bit)))
532
533        or, equivalently,
534
535      * exchange the bits in the positions given in each of the pairs X, Y,
536        ..., where each Y = X + BIT."
537
538   (let* ((n (length p))
539          (w (ash 1 (integer-length (1- n))))
540          (p (let ((new (make-array w :element-type 'fixnum)))
541               (replace new p)
542               (do ((i n (1+ i)))
543                   ((>= i w))
544                 (setf (aref new i) i))
545               new))
546          (p-inv (invert-permutation p)))
547
548     (labels ((recurse (todo)
549                ;; Main recursive search.  DONE is a mask of the bits which
550                ;; have been searched.  Return the number of skipped stages
551                ;; and a list of steps (BIT M0 M1), indicating that (BIT M0)
552                ;; should be performed before the following stages, and
553                ;; (BIT M1) should be performed afterwards.
554                ;;
555                ;; The permutation `p' and its inverse `p-inv' will be
556                ;; modified and restored.
557
558                (cond ((zerop (logand todo (1- todo)))
559                       ;; Only one more bit left.  Use the more efficient
560                       ;; final-step computation.
561
562                       (let ((m (compute-final-benes-step n p p-inv todo)))
563                         (values (if m 0 1) (list (list todo m 0)))))
564
565                      (t
566                       ;; More searching to go.  We'll keep the result which
567                       ;; maximizes the number of skipped stages.
568                       (let ((best-list nil)
569                             (best-skips -1))
570
571                         (flet ((try (bit clear-input)
572                                  ;; Try a permutation with the given BIT and
573                                  ;; CLEAR-INPUT arguments to
574                                  ;; `compute-benes-step'.
575
576                                  ;; Compute the next step.
577                                  (multiple-value-bind (m0 m1)
578                                      (compute-benes-step n p p-inv
579                                                          bit clear-input)
580
581                                    ;; Apply the step and recursively
582                                    ;; determine the inner permutation.
583                                    (apply-benes-step p p-inv bit m0 m1)
584                                    (multiple-value-bind (nskip tail)
585                                        (recurse (logandc2 todo bit))
586                                      (apply-benes-step p p-inv bit m0 m1)
587
588                                      ;; Work out how good this network is.
589                                      ;; Keep it if it improves over the
590                                      ;; previous attempt.
591                                      (when (zerop m0) (incf nskip))
592                                      (when (zerop m1) (incf nskip))
593                                      (when (> nskip best-skips)
594                                        (setf best-list
595                                                (cons (list bit m0 m1)
596                                                      tail)
597                                              best-skips
598                                                nskip))))))
599
600                           ;; Work through each bit that we haven't done
601                           ;; already, and try skipping both the start and end
602                           ;; steps.
603                           (do ((bit 1 (ash bit 1)))
604                               ((>= bit w))
605                             (when (plusp (logand bit todo))
606                               (try bit nil)
607                               (try bit t))))
608                         (values best-skips best-list))))))
609
610       ;; Find the best permutation network.
611       (multiple-value-bind (nskip list) (recurse (1- w))
612         (declare (ignore nskip))
613
614         ;; Turn the list returned by `recurse' into a list of (SHIFT MASK)
615         ;; entries as expected by everything else.
616         (let ((head nil) (tail nil))
617           (dolist (step list (nconc (nreverse head) tail))
618             (destructuring-bind (bit m0 m1) step
619               (when (plusp m0) (push (cons bit m0) head))
620               (when (plusp m1) (push (cons bit m1) tail)))))))))
621
622 ;;;--------------------------------------------------------------------------
623 ;;; Special functions for DES permutations.
624
625 (defun benes-search-des (p &optional attempts)
626   "Search for a Beneš network for a DES 64-bit permutation.
627
628    P must be a sequence of 64 fixnums, each of which is between 0 and 64
629    inclusive.  In the DES convention, bits are numbered with the most-
630    significant bit being bit 1, and increasing towards the least-significant
631    bit, which is bit 64.  Each nonzero number must appear at most once, and
632    specifies which input bit must appear in that output position.  There may
633    also be any number of zero entries, which mean `don't care'.
634
635    This function searches for and returns a Beneš network which implements a
636    satisfactory permutation.  If ATTEMPTS is nil or omitted, then search
637    exhaustively, returning the shortest network.  Otherwise, return the
638    shortest network found after considering ATTEMPTS randomly chosen
639    matching permutations."
640
641   (let* ((n (length p))
642          (p (map '(vector fixnum)
643                  (lambda (x)
644                    (if (zerop x) -1
645                        (- 64 x)))
646                  (reverse p)))
647          (seen (make-hash-table))
648          (nmissing 0) (missing nil) (indices nil))
649
650     ;; Find all of the `don't care' slots, and keep track of the bits which
651     ;; have homes to go to.
652     (dotimes (i n)
653       (let ((x (aref p i)))
654         (cond ((minusp x)
655                (push i indices)
656                (incf nmissing))
657               (t (setf (gethash x seen) t)))))
658
659     ;; Fill in numbers of the input bits which don't have fixed places to go.
660     (setf missing (make-array nmissing :element-type 'fixnum))
661     (let ((j 0))
662       (dotimes (i n)
663         (unless (gethash i seen)
664           (setf (aref missing j) i)
665           (incf j)))
666       (assert (= j nmissing)))
667
668     ;; Run the search, printing successes as we find them to keep the user
669     ;; amused.
670     (let ((best nil) (best-length nil))
671       (loop
672         (cond ((eql attempts 0) (return best))
673               (attempts (shuffle missing) (decf attempts))
674               ((null (next-permutation missing)) (return best)))
675         (do ((idx indices (cdr idx))
676              (i 0 (1+ i)))
677             ((endp idx))
678           (setf (aref p (car idx)) (aref missing i)))
679         (let* ((benes (benes-search p)) (len (length benes)))
680           (when (or (null best-length)
681                     (< len best-length))
682             (setf best-length len
683                   best benes)
684             (print-permutation-network benes)
685             (force-output)))))))
686
687 ;;;--------------------------------------------------------------------------
688 ;;; Examples and useful runes.
689
690 #+example
691 (let* ((ip #(58 50 42 34 26 18 10  2
692              60 52 44 36 28 20 12  4
693              62 54 46 38 30 22 14  6
694              64 56 48 40 32 24 16  8
695              57 49 41 33 25 17  9  1
696              59 51 43 35 27 19 11  3
697              61 53 45 37 29 21 13  5
698              63 55 47 39 31 23 15  7))
699        (fixed-ip (map '(vector fixnum)
700                       (lambda (x) (- 64 x))
701                       (reverse ip)))
702
703        ;; The traditional network.
704        (trad-network
705         (make-permutation-network
706          64                             ;  5  4  3  2  1  0
707          '((:exchange-invert 2 5)       ; ~2  4  3 ~5  1  0
708            (:exchange-invert 1 4)       ; ~2 ~1  3 ~5 ~4  0
709            (:exchange-invert 0 3)       ; ~2 ~1 ~0 ~5 ~4 ~3
710            (:exchange-invert 3 4)       ; ~2  0  1 ~5 ~4 ~3
711            (:exchange-invert 4 5))))    ; ~0  2  1 ~5 ~4 ~3
712
713        ;; The new twizzle-optimized network.
714        (new-network
715         (make-permutation-network
716          64                             ;  5  4  3  2  1  0
717          '((:exchange-invert 2 5)       ; ~2  4  3 ~5  1  0
718            (:exchange-invert 4 5)       ; ~4  2  3 ~5  1  0
719            (:exchange        1 5)       ;  1  2  3 ~5 ~4  0
720            (:exchange        3 5)       ;  3  2  1 ~5 ~4  0
721            (:exchange-invert 0 5)))))   ; ~0  2  1 ~5 ~4 ~3
722
723   (fresh-line)
724
725   (let ((benes-network (benes-search fixed-ip)))
726     (print-permutation-network benes-network)
727     (demonstrate-permutation-network 64 benes-network :reference fixed-ip))
728   (terpri)
729   (print-permutation-network trad-network)
730   (demonstrate-permutation-network 64 trad-network :reference fixed-ip)
731   (terpri)
732   (print-permutation-network new-network)
733   (demonstrate-permutation-network 64 new-network :reference fixed-ip))
734
735 #+example
736 (benes-search-des #( 0  0  0  0
737                     57 49 41 33 25 17  9  1
738                     58 50 42 34 26 18 10  2
739                     59 51 43 35 27 19 11  3
740                     60 52 44 36
741                      0  0  0  0
742                     63 55 47 39 31 23 15  7
743                     62 54 46 38 30 22 14  6
744                     61 53 45 37 29 21 13  5
745                                 28 20 12  4))
746
747 #+example
748 (let ((pc2 (make-array '(8 6)
749                        :element-type 'fixnum
750                        :initial-contents '((14 17 11 24  1  5)
751                                            ( 3 28 15  6 21 10)
752                                            (23 19 12  4 26  8)
753                                            (16  7 27 20 13  2)
754                                            (41 52 31 37 47 55)
755                                            (30 40 51 45 33 48)
756                                            (44 49 39 56 34 53)
757                                            (46 42 50 36 29 32)))))
758   (benes-search-des
759    (make-array 64
760                :element-type 'fixnum
761                :initial-contents
762                  (loop for i in '(2 4 6 8 1 3 5 7)
763                        nconc (list 0 0)
764                        nconc (loop for j below 6
765                                    for x = (aref pc2 (1- i) j)
766                                    collect (if (<= x 32) (+ x 4) (+ x 8)))))
767    1000))