symm/des-base.h: Improve the IP permutation network.

[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 &optional reference (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."

  (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)))

	(dotimes (i n) (push (car (aref v i)) (aref v i)))

	(dotimes (i n)
	  (when (logbitp i mask)
	    (rotatef (car (aref v i))
		     (car (aref v (+ i shift))))))))

    ;; Print the result.
    (let ((ok (not (null reference))))
      (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)))

;;;--------------------------------------------------------------------------
;;; 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)))
       (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
       (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)

  (print-permutation-network trad-network)
  (demonstrate-permutation-network 64 trad-network fixed-ip)
  (terpri)
  (print-permutation-network new-network)
  (demonstrate-permutation-network 64 new-network fixed-ip))