X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/d3f33b9a37f45bc5809af314b1b436b712d59a80..e7ee4000d06cadd8355404c8ddfb3d16265d24ca:/utils/permute.lisp

diff --git a/utils/permute.lisp b/utils/permute.lisp
index 372f5cfc..884068ac 100644
--- a/utils/permute.lisp
+++ b/utils/permute.lisp
@@ -205,50 +205,484 @@
     (terpri stream)))
 
 (defun demonstrate-permutation-network
-    (n steps &optional reference (stream *standard-output*))
+    (n steps
+     &key reference
+	  offsets
+	  (stream *standard-output*))
   "Print, on STREAM, a demonstration of the permutation STEPS.
 
-   Begin, on the left, with the integers from 0 up to N - 1.  For each
-   (SHIFT . MASK) element in STEPS, print an additional column showing the
-   effect of that step on the vector.  If REFERENCE is not nil, then it
-   should be a vector of length at least N: on the right, print the REFERENCE
-   vector, showing where the result of the permutation STEPS differs from the
-   REFERENCE.  Return non-nil if the output matches the reference; return nil
-   if the output doesn't match, or no reference was supplied."
+   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)
 
-  (let ((v (make-array n)))
-
-    ;; Initialize a vector of lists which will record, for each step in the
-    ;; permutation network, which value is in that position.  The lists are
-    ;; reversed, so the `current' value is at the front.
-    (dotimes (i n) (setf (aref v i) (cons i nil)))
-
-    ;; Work through the permutation steps updating the vector.
-    (dolist (step steps)
-      (let ((shift (car step)) (mask (cdr step)))
+	(dolist (step steps)
+	  (let ((shift (car step)) (mask (cdr step)))
+	    (apply-step shift mask v)
+	    (format stream ";;~%")
+	    (show-stage shift mask v)))
 
-	(dotimes (i n) (push (car (aref v i)) (aref v i)))
+	(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)))))
 
-	(dotimes (i n)
-	  (when (logbitp i mask)
-	    (rotatef (car (aref v i))
-		     (car (aref v (+ i shift))))))))
+;;;--------------------------------------------------------------------------
+;;; 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)))))))))
 
-    ;; Print the result.
-    (let ((ok (not (null reference))))
+;;;--------------------------------------------------------------------------
+;;; 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)
-	(let* ((entry (aref v i))
-	       (final (car entry)))
-	  (format stream "~{ ~7D~}" (reverse entry))
-	  (when reference
-	    (let* ((want (aref reference i))
-		   (match (eql final want)))
-	      (format stream " ~:[/=~;==~] ~7D" match want)
-	      (unless match (setf ok nil))))
-	  (terpri stream)))
-      (when reference
-	(format stream "~:[FAIL~;pass~]~%" ok))
-      ok)))
+	(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.
@@ -265,6 +699,8 @@
        (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
@@ -272,9 +708,60 @@
 	   (: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
+	   (: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 fixed-ip))
+  (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))