utils/permute.lisp (demonstrate-permutation-network): Optionally show offsets.

[catacomb] / utils / permute.lisp

;;; -*-lisp-*-

;;; This file isn't a program as such: rather, it's a collection of handy
;;; functions which can be used in an interactive session.

;;;--------------------------------------------------------------------------
;;; General permutation utilities.

(defun shuffle (v)
  "Randomly permute the elements of the vector V.  Return V."
  (let ((n (length v)))
    (do ((k n (1- k)))
	((<= k 1) v)
      (let ((i (random k)))
	(unless (= i (1- k))
	  (rotatef (aref v i) (aref v (1- k))))))))

(defun identity-permutation (n)
  "Return the do-nothing permutation on N elements."
  (let ((v (make-array n :element-type 'fixnum)))
    (dotimes (i n v) (setf (aref v i) i))))

(defun invert-permutation (p)
  "Given a permutation P, return its inverse."
  (let* ((n (length p)) (p-inv (make-array n :element-type 'fixnum)))
    (dotimes (i n) (setf (aref p-inv (aref p i)) i))
    p-inv))

(defun next-permutation (v)
  "Adjust V so that it reflects the next permutation in ascending order.

   V should be a vector of real numbers.  Returns V if successful, or nil if
   there are no more permutations."

  ;; The tail of the vector consists of a sequence ... A, Z, Y, X, ..., where
  ;; Z > Y > X ... is in reverse order, and A < Z.  The next permutation is
  ;; then the smallest out of Z, Y, X, ... which is larger than A, followed
  ;; by the remaining elements in ascending order.
  ;;
  ;; Equivalently, reverse the tail Z, Y, X, ... so we have A, ... X, Y, Z,
  ;; and swap A with the next larger element.

  (let ((n (length v)))
    (cond ((< n 2) nil)
	  (t (let* ((k (1- n))
		    (x (aref v k)))
	       (loop (when (zerop k) (return-from next-permutation nil))
		     (decf k)
		     (let ((y (aref v k)))
		       (when (prog1 (< y x)
			       (setf x y))
			 (return))))
	       (do ((i (1+ k) (1+ i))
		    (j (1- n) (1- j)))
		   ((> i j))
		 (rotatef (aref v i) (aref v j)))
	       (do ((i (- n 2) (1- i)))
		   ((or (<= i k) (< (aref v i) x))
		    (rotatef (aref v k) (aref v (1+ i)))))
	       v)))))

(defun make-index-mask (w mask-expr)
  "Construct a bitmask based on bitwise properties of the bit indices.

   The function returns a W-bit mask in which each bit is set if MASK-EXPR
   of true of the bit's index.  MASK-EXPR may be one of the following:

     * I -- an integer I is true if bit I of the bit index is set;
     * (not EXPR) -- is true if EXPR is false;
     * (and EXPR EXPR ...) -- is true if all of the EXPRs are true; and
     * (or EXPR EXPR ...) -- is true if any of the EXPRs is true."

  (let ((max-bit (1- (integer-length (1- w))))
	(mask 0))
    (dotimes (i w mask)
      (labels ((interpret (expr)
		 (cond ((and (integerp expr) (<= 0 expr max-bit))
			(logbitp expr i))
		       ((and (consp expr) (eq (car expr) 'not)
			     (null (cddr expr)))
			(not (interpret (cadr expr))))
		       ((and (consp expr) (eq (car expr) 'and))
			(every #'interpret (cdr expr)))
		       ((and (consp expr) (eq (car expr) 'or))
			(some #'interpret (cdr expr)))
		       (t
			(error "unknown mask expression ~S" expr)))))
	(when (interpret mask-expr)
	  (setf (ldb (byte 1 i) mask) 1))))))

(defun make-permutation-network (w steps)
  "Construct a permutation network.

   The integer W gives the number of bits to be acted upon.  The STEPS are a
   list of instructions of the following forms:

     * (SHIFT . MASK) -- a pair of integers is treated literally;

     * (SHIFT MASK-EXPR) -- the SHIFT is literal, but the MASK-EXPR is
       processed by `make-index-mask' to calculate the mask;

     * (:invert I) -- make an instruction which inverts the sense of the
       index bit I;

     * (:exchange I J) -- make an instruction which exchanges index bits I
       and J; or

     * (:exchange-invert I J) -- make an instruction which exchanges and
       inverts index bits I and J.

   The output is a list of primitive (SHIFT . MASK) steps, indicating that
   the bits of the input selected by MASK are to be swapped with the bits
   selected by (ash MASK SHIFT)."

  (let ((max-mask (1- (ash 1 w)))
	(max-shift (1- w))
	(max-bit (1- (integer-length (1- w))))
	(list nil))
    (dolist (step steps)
      (cond ((and (consp step)
		  (integerp (car step)) (<= 0 (car step) max-shift)
		  (integerp (cdr step)) (<= 0 (cdr step) max-mask))
	     (push step list))
	    ((and (consp step)
		  (integerp (car step)) (<= 0 (car step) max-shift)
		  (null (cddr step)))
	     (push (cons (car step) (make-index-mask w (cadr step))) list))
	    ((and (consp step)
		  (eq (car step) :invert)
		  (integerp (cadr step)) (<= 0 (cadr step) max-bit)
		  (null (cddr step)))
	     (let ((i (cadr step)))
	       (push (cons (ash 1 i) (make-index-mask w `(not ,i))) list)))
	    ((and (consp step)
		  (eq (car step) :exchange)
		  (integerp (cadr step)) (integerp (caddr step))
		  (<= 0 (cadr step) (caddr step) max-bit)
		  (null (cdddr step)))
	     (let ((i (cadr step)) (j (caddr step)))
	       (push (cons (- (ash 1 j) (ash 1 i))
			   (make-index-mask w `(and ,i (not ,j))))
		     list)))
	    ((and (consp step)
		  (eq (car step) :exchange-invert)
		  (integerp (cadr step)) (integerp (caddr step))
		  (<= 0 (cadr step) (caddr step) max-bit)
		  (null (cdddr step)))
	     (let ((i (cadr step)) (j (caddr step)))
	       (push (cons (+ (ash 1 i) (ash 1 j))
			   (make-index-mask w `(and (not ,i) (not ,j))))
		     list)))
	    (t
	     (error "unknown permutation step ~S" step))))
    (nreverse list)))

;;;--------------------------------------------------------------------------
;;; Permutation network diagnostics.

(defun print-permutation-network (steps &optional (stream *standard-output*))
  "Print a description of the permutation network STEPS to STREAM.

   A permutation network consists of a list of pairs

	(SHIFT . MASK)

   indicating that the bits selected by MASK, and those SHIFT bits to the
   left, should be exchanged.

   The output is intended to be human-readable and is subject to change."

  (let ((shiftwd 1) (maskwd 2))

    ;; Determine suitable print widths for shifts and masks.
    (dolist (step steps)
      (let ((shift (car step)) (mask (cdr step)))
	(let ((swd (1+ (floor (log shift 10))))
	      (mwd (ash 1 (- (integer-length (1- (integer-length mask)))
			     2))))
	  (when (> swd shiftwd) (setf shiftwd swd))
	  (when (> mwd maskwd) (setf maskwd mwd)))))

    ;; Print the display.
    (pprint-logical-block (stream steps :prefix "(" :suffix ")")
      (let ((first t))
	(dolist (step steps)
	  (let ((shift (car step)) (mask (cdr step)))

	    ;; Separate entries with newlines.
	    (cond (first (setf first nil))
		  (t (pprint-newline :mandatory stream)))

	    (let ((swaps nil))

	      ;; Determine the list of exchanges implied by the mask.
	      (dotimes (i (integer-length mask))
		(when (logbitp i mask)
		  (push (cons i (+ i shift)) swaps)))
	      (setf swaps (nreverse swaps))

	      ;; Print the entry.
	      (format stream "~@<(~;~vD #x~(~v,'0X~) ~8I~:@_~W~;)~:>"
		      shiftwd shift maskwd mask swaps))))))

    ;; Print a final newline following the close parenthesis.
    (terpri stream)))

(defun demonstrate-permutation-network
    (n steps
     &key reference
	  offsets
	  (stream *standard-output*))
  "Print, on STREAM, a demonstration of the permutation STEPS.

   The output is intended to be useful to human readers and is subject to
   change.  Currently, it prints a sequence of diagrams on STREAM.  The left
   hand side of each row shows a map of which bits are affected: `-' means
   that the bit remains in the same position, `*' means that it moves
   forward, and `#' means that it moves back; each `*' pairs with the
   earliest unpaired `#' marker.  The right hand side shows the arrangement
   of the original input bits.

   If OFFSETS is non-nil, then also print a table of offsets showing how far
   each bit has yet to move.

   If REFERENCE is not nil, then print a final pair of diagrams.  This time,
   the map shows `-' for correct bits and `x' for incorrect ones, with the
   right hand side showing the expected arrangement of input bits.

   The function returns non-nil if the STEPS resulted in the REFERENCE
   permutation, and nil if either the STEPS are incorrect or no REFERENCE was
   provided."

  (flet ((apply-step (shift mask v)
	   (dotimes (k n)
	     (when (logbitp k mask)
	       (rotatef (aref v k) (aref v (+ k shift)))))))

    (let* ((v (identity-permutation n))
	   (end (or reference
		    (let ((e (identity-permutation n)))
		      (dolist (step steps e)
			(let ((shift (car step)) (mask (cdr step)))
			  (apply-step shift mask e))))))
	   (end-inv (invert-permutation end))
	   (mapwd (ceiling (sqrt n)))
	   (ixwd (length (format nil "~D" (1- n)))))

      (flet ((show-stage (shift mask v)
	       (do ((i 0 (+ i mapwd)))
		   ((>= i n))
		 (write-string ";;	" stream)
		 (dotimes (j mapwd)
		   (let ((k (+ i j)))
		     (when (plusp j) (write-char #\space stream))
		     (write-char (cond ((>= k n)
					#\space)
				       ((logbitp k mask)
					#\*)
				       ((and (>= k shift)
					     (logbitp (- k shift) mask))
					#\#)
				       (t
					#\-))
				 stream)))
		 (write-string " | " stream)

		 (when offsets
		   (dotimes (j mapwd)
		     (let ((k (+ i j)))
		       (when (plusp j) (write-char #\space stream))
		       (cond ((>= k n)
			      (format stream "~v@T" (1+ ixwd)))
			     (t
			      (format stream "~*~[~2:*~vD~:;~2:*~v@D~]"
				      (1+ ixwd)
				      (- (aref end-inv (aref v k)) k))))))
		   (write-string " | " stream))

		 (dotimes (j (min mapwd (- n i)))
		   (let ((k (+ i j)))
		     (when (plusp j) (write-char #\space stream))
		     (format stream "~vD" ixwd (aref v k))))
		 (terpri))))

	(fresh-line)
	(show-stage 0 0 v)

	(dolist (step steps)
	  (let ((shift (car step)) (mask (cdr step)))
	    (apply-step shift mask v)
	    (format stream ";;~%")
	    (show-stage shift mask v)))

	(let ((ok (not (null reference))))
	  (when reference
	    (format stream ";;~%")
	    (do ((i 0 (+ i mapwd)))
		((>= i n))
	      (write-string ";;	" stream)
	      (dotimes (j mapwd)
		(let ((k (+ i j)))
		  (when (plusp j) (write-char #\space stream))
		  (write-char (cond ((>= k n) #\space)
				    ((/= (aref v k) (aref reference k)) #\x)
				    (t #\-))
			      stream)))
	      (write-string " | " stream)
	      (dotimes (j (min mapwd (- n i)))
		(let ((k (+ i j)))
		  (when (plusp j) (write-char #\space stream))
		  (format stream "~vD" ixwd (aref reference k))))
	      (terpri))
	    (unless (every #'= v reference)
	      (setf ok nil))
	    (format stream "~:[FAIL~;pass~]~%" ok))
	  ok)))))

;;;--------------------------------------------------------------------------
;;; Beneš networks.

(defun compute-benes-step (n p p-inv bit clear-input)
  "Compute a single layer of a Beneš network.

   N is a fixnum.  P is a vector of fixnums defining a permutation: for each
   output bit position i (numbering the least significant bit 0), element i
   gives the number of the input which should end up in that position; P-INV
   gives the inverse permutation in the same form.  BIT is a power of 2 which
   gives the distance between bits we should consider.  CLEAR-INPUT is
   a (generalized) boolean: if true, we attempt to do no work on the input
   side; if false, we try to do no work on the output side.  The length of P
   must be at least (logior N BIT).

   The output consists of a pair of masks M0 and M1, to be used on the input
   and output sides respectively.  The permutation stage, applied to an input
   X, is as follows:

	(let ((tmp (logand (logxor x (ash x (- bit))) mask)))
	  (logxor x tmp (ash tmp bit)))

   The critical property of the masks is that it's possible to compute an
   inner permutation, mapping the output of the M0 stage to the input of the
   M1 stage, such that (a) the overall composition of the three permutations
   is the given permutation P, and (b) the distances that the bits need to
   be moved by the inner permutation all have BIT clear.

   The resulting permutation will attempt to avoid touching elements at
   indices greater than N.  This attempt will succeed if all such elements
   contain values equal to their indices.

   By appropriately composing layers computed by this function, then, it's
   possible to perform an arbitrary permutation of 2^n bits in 2 n - 1 simple
   steps like the one above."

  ;; Consider the following problem.  You're given two equally-sized
  ;; collections of things, called `left' and `right'.  Each left thing is
  ;; attached to exactly one right thing with a string, and /vice versa/.
  ;; Furthermore, the left things, and the right things, are each linked
  ;; together in pairs, so each pair has two strings coming out of it.  Our
  ;; job is to paint the strings so that each linked pair of things has one
  ;; red string and one blue string.
  ;;
  ;; This is quite straightforward.  Pick a pair whose strings aren't yet
  ;; coloured, and colour one of its strings, chosen arbitrarily, red.  Find
  ;; the pair at the other end of this red string.  If the two other things
  ;; in these two pairs are connected, then paint that string blue and move
  ;; on.  Otherwise, both things have an uncoloured string, so paint both of
  ;; them blue and trace along these now blue strings to find two more thing
  ;; pairs.  Again, the other thing in each pair has an uncoloured string.
  ;; If these things share the /same/ string, paint it red and move on.
  ;; Otherwise, paint both strings red, trace, and repeat.  Eventually, we'll
  ;; find two things joined by the same string, each paired with another
  ;; thing whose strings we've just painted the same colour.  Once we're
  ;; done, we'll have painted a cycle's worth of strings, and each pair of
  ;; things will have either both of its strings painted different colours,
  ;; or neither of them.
  ;;
  ;; The right things are the integers 0, 1, ..., n - 1, where n is some
  ;; power of 2, paired according to whether they differ only in BIT.  The
  ;; left things are the same integers, connected to the right things
  ;; according to the permutation P: the right thing labelled i is connected
  ;; to the left thing labelled P(i).  Similarly, two left things are paired
  ;; if their labels P(i) and P(j) differ only in BIT.  We're going to paint
  ;; a string red if we're going to arrange to clear BIT in the labels at
  ;; both ends, possibly by swapping the two labels, and paint it red if
  ;; we're going to arrange to set BIT.  Once we've done this, later stages
  ;; of the filter will permute the red- and blue-painted things
  ;; independently.

  (let ((m0 0) (m1 0) (done 0))

    ;; Now work through the permutation cycles.
    (do ((i (1- n) (1- i)))
	((minusp i))

      ;; Skip if we've done this one already.
      (unless (or (plusp (logand i bit))
		  (logbitp i done))

	;; Find the other associated values.
	(let* ((i0 i) (i1 (aref p-inv i))
	       (sense (cond ((>= (logior i0 bit) n) 0)
			    (clear-input (logand i0 bit))
			    (t (logand i1 bit)))))

	  #+noise
	  (format t ";; new cycle: i0 = ~D, j0 = ~D; i1 = ~D, j1 = ~D~%"
		  i0 (logxor i0 bit)
		  i1 (logxor i1 bit))

	  ;; Mark this index as done.
	  (setf (ldb (byte 1 i0) done) 1)
	  #+noise (format t ";;   done = #x~2,'0X~%" done)

	  ;; Now trace round the cycle.
	  (loop

	    ;; Mark this index as done.
	    (setf (ldb (byte 1 (logandc2 i0 bit)) done) 1)
	    #+noise (format t ";;   done = #x~2,'0X~%" done)

	    ;; Swap the input and output pairs if necessary.
	    (unless (= (logand i0 bit) sense)
	      #+noise
	      (format t ";;   swap input:  ~D <-> ~D~%"
		      (logandc2 i0 bit) (logior i0 bit))
	      (setf (ldb (byte 1 (logandc2 i0 bit)) m0) 1))
	    (unless (= (logand i1 bit) sense)
	      #+noise
	      (format t ";;   swap output: ~D <-> ~D~%"
		      (logandc2 i1 bit) (logior i1 bit))
	      (setf (ldb (byte 1 (logandc2 i1 bit)) m1) 1))

	    ;; Advance around the cycle.
	    (let* ((j0 (logxor i0 bit))
		   (j1 (logxor i1 bit))
		   (next-i1 (aref p-inv j0))
		   (next-i0 (aref p j1)))
	      (when (= next-i0 j0) (return))
	      (setf i0 next-i0
		    i1 next-i1
		    sense (logxor sense bit)))

	    #+noise
	    (format t ";; advance: i0 = ~D, j0 = ~D; i1 = ~D, j1 = ~D~%"
		    i0 (logxor i0 bit)
		    i1 (logxor i1 bit))))))

    (values m0 m1)))

(defun compute-final-benes-step (n p p-inv bit)
  "Determine the innermost stage of a Beneš network.

   N is a fixnum.  P is a vector of fixnums defining a permutation: for each
   output bit position i (numbering the least significant bit 0), element i
   gives the number of the input which should end up in that position; P-INV
   gives the inverse permutation in the same form.  BIT is a power of 2 which
   gives the distance between bits we should consider.  The length of P must
   be at least (logior N BIT).

   Furthermore, the ith element of P must be equal either to i or to
   (logxor i BIT); and therefore P-INV must be equal to P.

   Return the mask such that

	(let ((tmp (logand (logxor x (ash x (- bit))) mask)))
	  (logxor x tmp (ash tmp bit)))

   applies the permutation P to the bits of x."

  (declare (ignorable p-inv))

  (let ((m 0))
    (dotimes (i n)
      (unless (plusp (logand i bit))
	(let ((x (aref p i)))
	  #+paranoid
	  (assert (= (logandc2 x bit) i))
	  #+paranoid
	  (assert (= x (aref p-inv i)))

	  (unless (= x i)
	    (setf (ldb (byte 1 i) m) 1)))))
    m))

(defun apply-benes-step (p p-inv bit m0 m1)
  "Apply input and output steps for a Beneš network to a permutation.

   Given the permutation P and its inverse, and the distance BIT, as passed
   to `compute-benes-step', and the masks M0 and M1 returned, determine and
   return the necessary `inner' permutation to be applied between these
   steps, and its inverse.

   A permutation-network step, and, in particular, a Beneš step, is an
   involution, so the change to the vectors P and P-INV can be undone by
   calling the function again with the same arguments."

  (flet ((swaps (p p-inv mask)
	   (dotimes (i0 (length p))
	     (when (logbitp i0 mask)
	       (let* ((j0 (logior i0 bit))
		      (i1 (aref p-inv i0))
		      (j1 (aref p-inv j0)))
		 (setf (aref p i1) j0
		       (aref p j1) i0)
		 (rotatef (aref p-inv i0) (aref p-inv j0)))))))
    (swaps p p-inv m0)
    (swaps p-inv p m1)

    #+paranoid
    (let* ((n (length p)))
      (dotimes (i n)
	(assert (= (aref p (aref p-inv i)) i))
	(assert (= (aref p-inv (aref p i)) i))))))

(defun benes-search (p)
  "Given a bit permutation P, describe a Beneš network implementing P.

   P is a sequence of fixnums defining a permutation: for each output bit
   position i (numbering the least significant bit 0), element i gives the
   number of the input which should end up in that position.

   The return value is a list of steps of the form

	(BIT MASK (X . Y) (X' . Y') ...)

   To implement this permutation step:

     * given an input X, compute

	(let ((tmp (logand (logxor x (ash x (- bit))) mask)))
	  (logxor x tmp (ash tmp bit)))

       or, equivalently,

     * exchange the bits in the positions given in each of the pairs X, Y,
       ..., where each Y = X + BIT."

  (let* ((n (length p))
	 (w (ash 1 (integer-length (1- n))))
	 (p (let ((new (make-array w :element-type 'fixnum)))
	      (replace new p)
	      (do ((i n (1+ i)))
		  ((>= i w))
		(setf (aref new i) i))
	      new))
	 (p-inv (invert-permutation p)))

    (labels ((recurse (todo)
	       ;; Main recursive search.  DONE is a mask of the bits which
	       ;; have been searched.  Return the number of skipped stages
	       ;; and a list of steps (BIT M0 M1), indicating that (BIT M0)
	       ;; should be performed before the following stages, and
	       ;; (BIT M1) should be performed afterwards.
	       ;;
	       ;; The permutation `p' and its inverse `p-inv' will be
	       ;; modified and restored.

	       (cond ((zerop (logand todo (1- todo)))
		      ;; Only one more bit left.  Use the more efficient
		      ;; final-step computation.

		      (let ((m (compute-final-benes-step n p p-inv todo)))
			(values (if m 0 1) (list (list todo m 0)))))

		     (t
		      ;; More searching to go.  We'll keep the result which
		      ;; maximizes the number of skipped stages.
		      (let ((best-list nil)
			    (best-skips -1))

			(flet ((try (bit clear-input)
				 ;; Try a permutation with the given BIT and
				 ;; CLEAR-INPUT arguments to
				 ;; `compute-benes-step'.

				 ;; Compute the next step.
				 (multiple-value-bind (m0 m1)
				     (compute-benes-step n p p-inv
							 bit clear-input)

				   ;; Apply the step and recursively
				   ;; determine the inner permutation.
				   (apply-benes-step p p-inv bit m0 m1)
				   (multiple-value-bind (nskip tail)
				       (recurse (logandc2 todo bit))
				     (apply-benes-step p p-inv bit m0 m1)

				     ;; Work out how good this network is.
				     ;; Keep it if it improves over the
				     ;; previous attempt.
				     (when (zerop m0) (incf nskip))
				     (when (zerop m1) (incf nskip))
				     (when (> nskip best-skips)
				       (setf best-list
					       (cons (list bit m0 m1)
						     tail)
					     best-skips
					       nskip))))))

			  ;; Work through each bit that we haven't done
			  ;; already, and try skipping both the start and end
			  ;; steps.
			  (do ((bit 1 (ash bit 1)))
			      ((>= bit w))
			    (when (plusp (logand bit todo))
			      (try bit nil)
			      (try bit t))))
			(values best-skips best-list))))))

      ;; Find the best permutation network.
      (multiple-value-bind (nskip list) (recurse (1- w))
	(declare (ignore nskip))

	;; Turn the list returned by `recurse' into a list of (SHIFT MASK)
	;; entries as expected by everything else.
	(let ((head nil) (tail nil))
	  (dolist (step list (nconc (nreverse head) tail))
	    (destructuring-bind (bit m0 m1) step
	      (when (plusp m0) (push (cons bit m0) head))
	      (when (plusp m1) (push (cons bit m1) tail)))))))))

;;;--------------------------------------------------------------------------
;;; Special functions for DES permutations.

(defun benes-search-des (p &optional attempts)
  "Search for a Beneš network for a DES 64-bit permutation.

   P must be a sequence of 64 fixnums, each of which is between 0 and 64
   inclusive.  In the DES convention, bits are numbered with the most-
   significant bit being bit 1, and increasing towards the least-significant
   bit, which is bit 64.  Each nonzero number must appear at most once, and
   specifies which input bit must appear in that output position.  There may
   also be any number of zero entries, which mean `don't care'.

   This function searches for and returns a Beneš network which implements a
   satisfactory permutation.  If ATTEMPTS is nil or omitted, then search
   exhaustively, returning the shortest network.  Otherwise, return the
   shortest network found after considering ATTEMPTS randomly chosen
   matching permutations."

  (let* ((n (length p))
	 (p (map '(vector fixnum)
		 (lambda (x)
		   (if (zerop x) -1
		       (- 64 x)))
		 (reverse p)))
	 (seen (make-hash-table))
	 (nmissing 0) (missing nil) (indices nil))

    ;; Find all of the `don't care' slots, and keep track of the bits which
    ;; have homes to go to.
    (dotimes (i n)
      (let ((x (aref p i)))
	(cond ((minusp x)
	       (push i indices)
	       (incf nmissing))
	      (t (setf (gethash x seen) t)))))

    ;; Fill in numbers of the input bits which don't have fixed places to go.
    (setf missing (make-array nmissing :element-type 'fixnum))
    (let ((j 0))
      (dotimes (i n)
	(unless (gethash i seen)
	  (setf (aref missing j) i)
	  (incf j)))
      (assert (= j nmissing)))

    ;; Run the search, printing successes as we find them to keep the user
    ;; amused.
    (let ((best nil) (best-length nil))
      (loop
	(cond ((eql attempts 0) (return best))
	      (attempts (shuffle missing) (decf attempts))
	      ((null (next-permutation missing)) (return best)))
	(do ((idx indices (cdr idx))
	     (i 0 (1+ i)))
	    ((endp idx))
	  (setf (aref p (car idx)) (aref missing i)))
	(let* ((benes (benes-search p)) (len (length benes)))
	  (when (or (null best-length)
		    (< len best-length))
	    (setf best-length len
		  best benes)
	    (print-permutation-network benes)
	    (force-output)))))))

;;;--------------------------------------------------------------------------
;;; Examples and useful runes.

#+example
(let* ((ip #(58 50 42 34 26 18 10  2
	     60 52 44 36 28 20 12  4
	     62 54 46 38 30 22 14  6
	     64 56 48 40 32 24 16  8
	     57 49 41 33 25 17  9  1
	     59 51 43 35 27 19 11  3
	     61 53 45 37 29 21 13  5
	     63 55 47 39 31 23 15  7))
       (fixed-ip (map '(vector fixnum)
		      (lambda (x) (- 64 x))
		      (reverse ip)))

       ;; The traditional network.
       (trad-network
	(make-permutation-network
	 64				;  5  4  3  2  1  0
	 '((:exchange-invert 2 5)	; ~2  4  3 ~5  1  0
	   (:exchange-invert 1 4)	; ~2 ~1  3 ~5 ~4  0
	   (:exchange-invert 0 3)	; ~2 ~1 ~0 ~5 ~4 ~3
	   (:exchange-invert 3 4)	; ~2  0  1 ~5 ~4 ~3
	   (:exchange-invert 4 5))))	; ~0  2  1 ~5 ~4 ~3

       ;; The new twizzle-optimized network.
       (new-network
	(make-permutation-network
	 64				;  5  4  3  2  1  0
	 '((:exchange-invert 2 5)	; ~2  4  3 ~5  1  0
	   (:exchange-invert 4 5)	; ~4  2  3 ~5  1  0
	   (:exchange	     1 5)	;  1  2  3 ~5 ~4  0
	   (:exchange	     3 5)	;  3  2  1 ~5 ~4  0
	   (:exchange-invert 0 5)))))	; ~0  2  1 ~5 ~4 ~3

  (fresh-line)

  (let ((benes-network (benes-search fixed-ip)))
    (print-permutation-network benes-network)
    (demonstrate-permutation-network 64 benes-network :reference fixed-ip))
  (terpri)
  (print-permutation-network trad-network)
  (demonstrate-permutation-network 64 trad-network :reference fixed-ip)
  (terpri)
  (print-permutation-network new-network)
  (demonstrate-permutation-network 64 new-network :reference fixed-ip))

#+example
(benes-search-des #( 0  0  0  0
		    57 49 41 33 25 17  9  1
		    58 50 42 34 26 18 10  2
		    59 51 43 35 27 19 11  3
		    60 52 44 36
		     0  0  0  0
		    63 55 47 39 31 23 15  7
		    62 54 46 38 30 22 14  6
		    61 53 45 37 29 21 13  5
				28 20 12  4))

#+example
(let ((pc2 (make-array '(8 6)
		       :element-type 'fixnum
		       :initial-contents '((14 17 11 24  1  5)
					   ( 3 28 15  6 21 10)
					   (23 19 12  4 26  8)
					   (16  7 27 20 13  2)
					   (41 52 31 37 47 55)
					   (30 40 51 45 33 48)
					   (44 49 39 56 34 53)
					   (46 42 50 36 29 32)))))
  (benes-search-des
   (make-array 64
	       :element-type 'fixnum
	       :initial-contents
		 (loop for i in '(2 4 6 8 1 3 5 7)
		       nconc (list 0 0)
		       nconc (loop for j below 6
				   for x = (aref pc2 (1- i) j)
				   collect (if (<= x 32) (+ x 4) (+ x 8)))))
   1000))