chiark - git - mdw - catacomb/blob - 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           (stream *standard-output*))
 211   "Print, on STREAM, a demonstration of the permutation STEPS.
 212
 213    Begin, on the left, with the integers from 0 up to N - 1.  For each
 214    (SHIFT . MASK) element in STEPS, print an additional column showing the
 215    effect of that step on the vector.  If REFERENCE is not nil, then it
 216    should be a vector of length at least N: on the right, print the REFERENCE
 217    vector, showing where the result of the permutation STEPS differs from the
 218    REFERENCE.  Return non-nil if the output matches the reference; return nil
 219    if the output doesn't match, or no reference was supplied."
 220
 221   (let ((v (make-array n)))
 222
 223     ;; Initialize a vector of lists which will record, for each step in the
 224     ;; permutation network, which value is in that position.  The lists are
 225     ;; reversed, so the `current' value is at the front.
 226     (dotimes (i n) (setf (aref v i) (cons i nil)))
 227
 228     ;; Work through the permutation steps updating the vector.
 229     (dolist (step steps)
 230       (let ((shift (car step)) (mask (cdr step)))
 231
 232         (dotimes (i n) (push (car (aref v i)) (aref v i)))
 233
 234         (dotimes (i n)
 235           (when (logbitp i mask)
 236             (rotatef (car (aref v i))
 237                      (car (aref v (+ i shift))))))))
 238
 239     ;; Print the result.
 240     (let ((ok (not (null reference))))
 241       (dotimes (i n)
 242         (let* ((entry (aref v i))
 243                (final (car entry)))
 244           (format stream "~{ ~7D~}" (reverse entry))
 245           (when reference
 246             (let* ((want (aref reference i))
 247                    (match (eql final want)))
 248               (format stream " ~:[/=~;==~] ~7D" match want)
 249               (unless match (setf ok nil))))
 250           (terpri stream)))
 251       (when reference
 252         (format stream "~:[FAIL~;pass~]~%" ok))
 253       ok)))
 254
 255 ;;;--------------------------------------------------------------------------
 256 ;;; Beneš networks.
 257
 258 (defun compute-benes-step (n p p-inv bit clear-input)
 259   "Compute a single layer of a Beneš network.
 260
 261    N is a fixnum.  P is a vector of fixnums defining a permutation: for each
 262    output bit position i (numbering the least significant bit 0), element i
 263    gives the number of the input which should end up in that position; P-INV
 264    gives the inverse permutation in the same form.  BIT is a power of 2 which
 265    gives the distance between bits we should consider.  CLEAR-INPUT is
 266    a (generalized) boolean: if true, we attempt to do no work on the input
 267    side; if false, we try to do no work on the output side.  The length of P
 268    must be at least (logior N BIT).
 269
 270    The output consists of a pair of masks M0 and M1, to be used on the input
 271    and output sides respectively.  The permutation stage, applied to an input
 272    X, is as follows:
 273
 274         (let ((tmp (logand (logxor x (ash x (- bit))) mask)))
 275           (logxor x tmp (ash tmp bit)))
 276
 277    The critical property of the masks is that it's possible to compute an
 278    inner permutation, mapping the output of the M0 stage to the input of the
 279    M1 stage, such that (a) the overall composition of the three permutations
 280    is the given permutation P, and (b) the distances that the bits need to
 281    be moved by the inner permutation all have BIT clear.
 282
 283    The resulting permutation will attempt to avoid touching elements at
 284    indices greater than N.  This attempt will succeed if all such elements
 285    contain values equal to their indices.
 286
 287    By appropriately composing layers computed by this function, then, it's
 288    possible to perform an arbitrary permutation of 2^n bits in 2 n - 1 simple
 289    steps like the one above."
 290
 291   ;; Consider the following problem.  You're given two equally-sized
 292   ;; collections of things, called `left' and `right'.  Each left thing is
 293   ;; attached to exactly one right thing with a string, and /vice versa/.
 294   ;; Furthermore, the left things, and the right things, are each linked
 295   ;; together in pairs, so each pair has two strings coming out of it.  Our
 296   ;; job is to paint the strings so that each linked pair of things has one
 297   ;; red string and one blue string.
 298   ;;
 299   ;; This is quite straightforward.  Pick a pair whose strings aren't yet
 300   ;; coloured, and colour one of its strings, chosen arbitrarily, red.  Find
 301   ;; the pair at the other end of this red string.  If the two other things
 302   ;; in these two pairs are connected, then paint that string blue and move
 303   ;; on.  Otherwise, both things have an uncoloured string, so paint both of
 304   ;; them blue and trace along these now blue strings to find two more thing
 305   ;; pairs.  Again, the other thing in each pair has an uncoloured string.
 306   ;; If these things share the /same/ string, paint it red and move on.
 307   ;; Otherwise, paint both strings red, trace, and repeat.  Eventually, we'll
 308   ;; find two things joined by the same string, each paired with another
 309   ;; thing whose strings we've just painted the same colour.  Once we're
 310   ;; done, we'll have painted a cycle's worth of strings, and each pair of
 311   ;; things will have either both of its strings painted different colours,
 312   ;; or neither of them.
 313   ;;
 314   ;; The right things are the integers 0, 1, ..., n - 1, where n is some
 315   ;; power of 2, paired according to whether they differ only in BIT.  The
 316   ;; left things are the same integers, connected to the right things
 317   ;; according to the permutation P: the right thing labelled i is connected
 318   ;; to the left thing labelled P(i).  Similarly, two left things are paired
 319   ;; if their labels P(i) and P(j) differ only in BIT.  We're going to paint
 320   ;; a string red if we're going to arrange to clear BIT in the labels at
 321   ;; both ends, possibly by swapping the two labels, and paint it red if
 322   ;; we're going to arrange to set BIT.  Once we've done this, later stages
 323   ;; of the filter will permute the red- and blue-painted things
 324   ;; independently.
 325
 326   (let ((m0 0) (m1 0) (done 0))
 327
 328     ;; Now work through the permutation cycles.
 329     (do ((i (1- n) (1- i)))
 330         ((minusp i))
 331
 332       ;; Skip if we've done this one already.
 333       (unless (or (plusp (logand i bit))
 334                   (logbitp i done))
 335
 336         ;; Find the other associated values.
 337         (let* ((i0 i) (i1 (aref p-inv i))
 338                (sense (cond ((>= (logior i0 bit) n) 0)
 339                             (clear-input (logand i0 bit))
 340                             (t (logand i1 bit)))))
 341
 342           #+noise
 343           (format t ";; new cycle: i0 = ~D, j0 = ~D; i1 = ~D, j1 = ~D~%"
 344                   i0 (logxor i0 bit)
 345                   i1 (logxor i1 bit))
 346
 347           ;; Mark this index as done.
 348           (setf (ldb (byte 1 i0) done) 1)
 349           #+noise (format t ";;   done = #x~2,'0X~%" done)
 350
 351           ;; Now trace round the cycle.
 352           (loop
 353
 354             ;; Mark this index as done.
 355             (setf (ldb (byte 1 (logandc2 i0 bit)) done) 1)
 356             #+noise (format t ";;   done = #x~2,'0X~%" done)
 357
 358             ;; Swap the input and output pairs if necessary.
 359             (unless (= (logand i0 bit) sense)
 360               #+noise
 361               (format t ";;   swap input:  ~D <-> ~D~%"
 362                       (logandc2 i0 bit) (logior i0 bit))
 363               (setf (ldb (byte 1 (logandc2 i0 bit)) m0) 1))
 364             (unless (= (logand i1 bit) sense)
 365               #+noise
 366               (format t ";;   swap output: ~D <-> ~D~%"
 367                       (logandc2 i1 bit) (logior i1 bit))
 368               (setf (ldb (byte 1 (logandc2 i1 bit)) m1) 1))
 369
 370             ;; Advance around the cycle.
 371             (let* ((j0 (logxor i0 bit))
 372                    (j1 (logxor i1 bit))
 373                    (next-i1 (aref p-inv j0))
 374                    (next-i0 (aref p j1)))
 375               (when (= next-i0 j0) (return))
 376               (setf i0 next-i0
 377                     i1 next-i1
 378                     sense (logxor sense bit)))
 379
 380             #+noise
 381             (format t ";; advance: i0 = ~D, j0 = ~D; i1 = ~D, j1 = ~D~%"
 382                     i0 (logxor i0 bit)
 383                     i1 (logxor i1 bit))))))
 384
 385     (values m0 m1)))
 386
 387 (defun compute-final-benes-step (n p p-inv bit)
 388   "Determine the innermost stage of a Beneš network.
 389
 390    N is a fixnum.  P is a vector of fixnums defining a permutation: for each
 391    output bit position i (numbering the least significant bit 0), element i
 392    gives the number of the input which should end up in that position; P-INV
 393    gives the inverse permutation in the same form.  BIT is a power of 2 which
 394    gives the distance between bits we should consider.  The length of P must
 395    be at least (logior N BIT).
 396
 397    Furthermore, the ith element of P must be equal either to i or to
 398    (logxor i BIT); and therefore P-INV must be equal to P.
 399
 400    Return the mask such that
 401
 402         (let ((tmp (logand (logxor x (ash x (- bit))) mask)))
 403           (logxor x tmp (ash tmp bit)))
 404
 405    applies the permutation P to the bits of x."
 406
 407   (declare (ignorable p-inv))
 408
 409   (let ((m 0))
 410     (dotimes (i n)
 411       (unless (plusp (logand i bit))
 412         (let ((x (aref p i)))
 413           #+paranoid
 414           (assert (= (logandc2 x bit) i))
 415           #+paranoid
 416           (assert (= x (aref p-inv i)))
 417
 418           (unless (= x i)
 419             (setf (ldb (byte 1 i) m) 1)))))
 420     m))
 421
 422 (defun apply-benes-step (p p-inv bit m0 m1)
 423   "Apply input and output steps for a Beneš network to a permutation.
 424
 425    Given the permutation P and its inverse, and the distance BIT, as passed
 426    to `compute-benes-step', and the masks M0 and M1 returned, determine and
 427    return the necessary `inner' permutation to be applied between these
 428    steps, and its inverse.
 429
 430    A permutation-network step, and, in particular, a Beneš step, is an
 431    involution, so the change to the vectors P and P-INV can be undone by
 432    calling the function again with the same arguments."
 433
 434   (flet ((swaps (p p-inv mask)
 435            (dotimes (i0 (length p))
 436              (when (logbitp i0 mask)
 437                (let* ((j0 (logior i0 bit))
 438                       (i1 (aref p-inv i0))
 439                       (j1 (aref p-inv j0)))
 440                  (setf (aref p i1) j0
 441                        (aref p j1) i0)
 442                  (rotatef (aref p-inv i0) (aref p-inv j0)))))))
 443     (swaps p p-inv m0)
 444     (swaps p-inv p m1)
 445
 446     #+paranoid
 447     (let* ((n (length p)))
 448       (dotimes (i n)
 449         (assert (= (aref p (aref p-inv i)) i))
 450         (assert (= (aref p-inv (aref p i)) i))))))
 451
 452 (defun benes-search (p)
 453   "Given a bit permutation P, describe a Beneš network implementing P.
 454
 455    P is a sequence of fixnums defining a permutation: for each output bit
 456    position i (numbering the least significant bit 0), element i gives the
 457    number of the input which should end up in that position.
 458
 459    The return value is a list of steps of the form
 460
 461         (BIT MASK (X . Y) (X' . Y') ...)
 462
 463    To implement this permutation step:
 464
 465      * given an input X, compute
 466
 467         (let ((tmp (logand (logxor x (ash x (- bit))) mask)))
 468           (logxor x tmp (ash tmp bit)))
 469
 470        or, equivalently,
 471
 472      * exchange the bits in the positions given in each of the pairs X, Y,
 473        ..., where each Y = X + BIT."
 474
 475   (let* ((n (length p))
 476          (w (ash 1 (integer-length (1- n))))
 477          (p (let ((new (make-array w :element-type 'fixnum)))
 478               (replace new p)
 479               (do ((i n (1+ i)))
 480                   ((>= i w))
 481                 (setf (aref new i) i))
 482               new))
 483          (p-inv (invert-permutation p)))
 484
 485     (labels ((recurse (todo)
 486                ;; Main recursive search.  DONE is a mask of the bits which
 487                ;; have been searched.  Return the number of skipped stages
 488                ;; and a list of steps (BIT M0 M1), indicating that (BIT M0)
 489                ;; should be performed before the following stages, and
 490                ;; (BIT M1) should be performed afterwards.
 491                ;;
 492                ;; The permutation `p' and its inverse `p-inv' will be
 493                ;; modified and restored.
 494
 495                (cond ((zerop (logand todo (1- todo)))
 496                       ;; Only one more bit left.  Use the more efficient
 497                       ;; final-step computation.
 498
 499                       (let ((m (compute-final-benes-step n p p-inv todo)))
 500                         (values (if m 0 1) (list (list todo m 0)))))
 501
 502                      (t
 503                       ;; More searching to go.  We'll keep the result which
 504                       ;; maximizes the number of skipped stages.
 505                       (let ((best-list nil)
 506                             (best-skips -1))
 507
 508                         (flet ((try (bit clear-input)
 509                                  ;; Try a permutation with the given BIT and
 510                                  ;; CLEAR-INPUT arguments to
 511                                  ;; `compute-benes-step'.
 512
 513                                  ;; Compute the next step.
 514                                  (multiple-value-bind (m0 m1)
 515                                      (compute-benes-step n p p-inv
 516                                                          bit clear-input)
 517
 518                                    ;; Apply the step and recursively
 519                                    ;; determine the inner permutation.
 520                                    (apply-benes-step p p-inv bit m0 m1)
 521                                    (multiple-value-bind (nskip tail)
 522                                        (recurse (logandc2 todo bit))
 523                                      (apply-benes-step p p-inv bit m0 m1)
 524
 525                                      ;; Work out how good this network is.
 526                                      ;; Keep it if it improves over the
 527                                      ;; previous attempt.
 528                                      (when (zerop m0) (incf nskip))
 529                                      (when (zerop m1) (incf nskip))
 530                                      (when (> nskip best-skips)
 531                                        (setf best-list
 532                                                (cons (list bit m0 m1)
 533                                                      tail)
 534                                              best-skips
 535                                                nskip))))))
 536
 537                           ;; Work through each bit that we haven't done
 538                           ;; already, and try skipping both the start and end
 539                           ;; steps.
 540                           (do ((bit 1 (ash bit 1)))
 541                               ((>= bit w))
 542                             (when (plusp (logand bit todo))
 543                               (try bit nil)
 544                               (try bit t))))
 545                         (values best-skips best-list))))))
 546
 547       ;; Find the best permutation network.
 548       (multiple-value-bind (nskip list) (recurse (1- w))
 549         (declare (ignore nskip))
 550
 551         ;; Turn the list returned by `recurse' into a list of (SHIFT MASK)
 552         ;; entries as expected by everything else.
 553         (let ((head nil) (tail nil))
 554           (dolist (step list (nconc (nreverse head) tail))
 555             (destructuring-bind (bit m0 m1) step
 556               (when (plusp m0) (push (cons bit m0) head))
 557               (when (plusp m1) (push (cons bit m1) tail)))))))))
 558
 559 ;;;--------------------------------------------------------------------------
 560 ;;; Special functions for DES permutations.
 561
 562 (defun benes-search-des (p &optional attempts)
 563   "Search for a Beneš network for a DES 64-bit permutation.
 564
 565    P must be a sequence of 64 fixnums, each of which is between 0 and 64
 566    inclusive.  In the DES convention, bits are numbered with the most-
 567    significant bit being bit 1, and increasing towards the least-significant
 568    bit, which is bit 64.  Each nonzero number must appear at most once, and
 569    specifies which input bit must appear in that output position.  There may
 570    also be any number of zero entries, which mean `don't care'.
 571
 572    This function searches for and returns a Beneš network which implements a
 573    satisfactory permutation.  If ATTEMPTS is nil or omitted, then search
 574    exhaustively, returning the shortest network.  Otherwise, return the
 575    shortest network found after considering ATTEMPTS randomly chosen
 576    matching permutations."
 577
 578   (let* ((n (length p))
 579          (p (map '(vector fixnum)
 580                  (lambda (x)
 581                    (if (zerop x) -1
 582                        (- 64 x)))
 583                  (reverse p)))
 584          (seen (make-hash-table))
 585          (nmissing 0) (missing nil) (indices nil))
 586
 587     ;; Find all of the `don't care' slots, and keep track of the bits which
 588     ;; have homes to go to.
 589     (dotimes (i n)
 590       (let ((x (aref p i)))
 591         (cond ((minusp x)
 592                (push i indices)
 593                (incf nmissing))
 594               (t (setf (gethash x seen) t)))))
 595
 596     ;; Fill in numbers of the input bits which don't have fixed places to go.
 597     (setf missing (make-array nmissing :element-type 'fixnum))
 598     (let ((j 0))
 599       (dotimes (i n)
 600         (unless (gethash i seen)
 601           (setf (aref missing j) i)
 602           (incf j)))
 603       (assert (= j nmissing)))
 604
 605     ;; Run the search, printing successes as we find them to keep the user
 606     ;; amused.
 607     (let ((best nil) (best-length nil))
 608       (loop
 609         (cond ((eql attempts 0) (return best))
 610               (attempts (shuffle missing) (decf attempts))
 611               ((null (next-permutation missing)) (return best)))
 612         (do ((idx indices (cdr idx))
 613              (i 0 (1+ i)))
 614             ((endp idx))
 615           (setf (aref p (car idx)) (aref missing i)))
 616         (let* ((benes (benes-search p)) (len (length benes)))
 617           (when (or (null best-length)
 618                     (< len best-length))
 619             (setf best-length len
 620                   best benes)
 621             (print-permutation-network benes)
 622             (force-output)))))))
 623
 624 ;;;--------------------------------------------------------------------------
 625 ;;; Examples and useful runes.
 626
 627 #+example
 628 (let* ((ip #(58 50 42 34 26 18 10  2
 629              60 52 44 36 28 20 12  4
 630              62 54 46 38 30 22 14  6
 631              64 56 48 40 32 24 16  8
 632              57 49 41 33 25 17  9  1
 633              59 51 43 35 27 19 11  3
 634              61 53 45 37 29 21 13  5
 635              63 55 47 39 31 23 15  7))
 636        (fixed-ip (map '(vector fixnum)
 637                       (lambda (x) (- 64 x))
 638                       (reverse ip)))
 639
 640        ;; The traditional network.  (Exchange each `*' with the earliest
 641        ;; available `#'.)
 642        ;;
 643        ;;       - - - - - - - -          0  1  2  3  4  5  6  7
 644        ;;       - - - - - - - -          8  9 10 11 12 13 14 15
 645        ;;       - - - - - - - -         16 17 18 19 20 21 22 23
 646        ;;       - - - - - - - -         24 25 26 27 28 29 30 31
 647        ;;       - - - - - - - -         32 33 34 35 36 37 38 39
 648        ;;       - - - - - - - -         40 41 42 43 44 45 46 47
 649        ;;       - - - - - - - -         48 49 50 51 52 53 54 55
 650        ;;       - - - - - - - -         56 57 58 59 60 61 62 63
 651        ;;
 652        ;;       * * * * - - - -         36 37 38 39  4  5  6  7
 653        ;;       * * * * - - - -         44 45 46 47 12 13 14 15
 654        ;;       * * * * - - - -         52 53 54 55 20 21 22 23
 655        ;;       * * * * - - - -         60 61 62 63 28 29 30 31
 656        ;;       - - - - # # # #         32 33 34 35  0  1  2  3
 657        ;;       - - - - # # # #         40 41 42 43  8  9 10 11
 658        ;;       - - - - # # # #         48 49 50 51 16 17 18 19
 659        ;;       - - - - # # # #         56 57 58 59 24 25 26 27
 660        ;;
 661        ;;       * * - - * * - -         54 55 38 39 22  3 26  7
 662        ;;       * * - - * * - -         62 63 46 47 30 11 34 15
 663        ;;       - - # # - - # #         52 53 36 37 20  1 24  5
 664        ;;       - - # # - - # #         60 61 44 45 28 19 22 13
 665        ;;       * * - - * * - -         50 51 34 35 18  9 12  3
 666        ;;       * * - - * * - -         58 59 42 43 26 17 20 11
 667        ;;       - - # # - - # #         48 49 32 33 16  7 10  1
 668        ;;       - - # # - - # #         56 57 40 41 24  5 28  9
 669        ;;
 670        ;;       * - * - * - * -         63 55 47 39 21 13 35  7
 671        ;;       - # - # - # - #         62 54 46 38 20 12 34  6
 672        ;;       * - * - * - * -         61 53 45 37 29 11 23  5
 673        ;;       - # - # - # - #         60 52 44 36 28 10 22  4
 674        ;;       * - * - * - * -         59 51 43 35 17 19 21  3
 675        ;;       - # - # - # - #         58 50 42 34 16 18 20  2
 676        ;;       * - * - * - * -         57 49 41 33 15  7 29  1
 677        ;;       - # - # - # - #         56 48 40 32 14  6 28  0
 678        ;;
 679        ;;       * * * * * * * *         60 52 44 36 28 20 12  4
 680        ;;       - - - - - - - -         62 54 46 38 30 22 14  6
 681        ;;       - - - - - - - -         61 53 45 37 29 21 13  5
 682        ;;       # # # # # # # #         63 55 47 39 31 23 15  7
 683        ;;       * * * * * * * *         56 48 40 32 24 16  8  0
 684        ;;       - - - - - - - -         58 50 42 34 26 18 10  2
 685        ;;       - - - - - - - -         57 49 41 33 25 17  9  1
 686        ;;       # # # # # # # #         59 51 43 35 27 19 11  3
 687        ;;
 688        ;;       * * * * * * * *         57 49 41 33 25 17  9  1
 689        ;;       * * * * * * * *         59 51 43 35 27 19 11  3
 690        ;;       - - - - - - - -         61 53 45 37 29 21 13  5
 691        ;;       - - - - - - - -         63 55 47 39 31 23 15  7
 692        ;;       - - - - - - - -         56 48 40 32 24 16  8  0
 693        ;;       - - - - - - - -         58 50 42 34 26 18 10  2
 694        ;;       # # # # # # # #         60 52 44 36 28 20 12  4
 695        ;;       # # # # # # # #         62 54 46 38 30 22 14  6
 696        (trad-network
 697         (make-permutation-network
 698          64                             ;  5  4  3  2  1  0
 699          '((:exchange-invert 2 5)       ; ~2  4  3 ~5  1  0
 700            (:exchange-invert 1 4)       ; ~2 ~1  3 ~5 ~4  0
 701            (:exchange-invert 0 3)       ; ~2 ~1 ~0 ~5 ~4 ~3
 702            (:exchange-invert 3 4)       ; ~2  0  1 ~5 ~4 ~3
 703            (:exchange-invert 4 5))))    ; ~0  2  1 ~5 ~4 ~3
 704
 705        ;; The new twizzle-optimized network.  (Exchange each `*' with the
 706        ;; earliest available `#'.)
 707        ;;
 708        ;;       - - - - - - - -          0  1  2  3  4  5  6  7
 709        ;;       - - - - - - - -          8  9 10 11 12 13 14 15
 710        ;;       - - - - - - - -         16 17 18 19 20 21 22 23
 711        ;;       - - - - - - - -         24 25 26 27 28 29 30 31
 712        ;;       - - - - - - - -         32 33 34 35 36 37 38 39
 713        ;;       - - - - - - - -         40 41 42 43 44 45 46 47
 714        ;;       - - - - - - - -         48 49 50 51 52 53 54 55
 715        ;;       - - - - - - - -         56 57 58 59 60 61 62 63
 716        ;;
 717        ;;       * * * * - - - -         36 37 38 39  4  5  6  7
 718        ;;       * * * * - - - -         44 45 46 47 12 13 14 15
 719        ;;       * * * * - - - -         52 53 54 55 20 21 22 23
 720        ;;       * * * * - - - -         60 61 62 63 28 29 30 31
 721        ;;       - - - - # # # #         32 33 34 35  0  1  2  3
 722        ;;       - - - - # # # #         40 41 42 43  8  9 10 11
 723        ;;       - - - - # # # #         48 49 50 51 16 17 18 19
 724        ;;       - - - - # # # #         56 57 58 59 24 25 26 27
 725        ;;
 726        ;;       * * * * * * * *         48 49 50 51 16 17 18 19
 727        ;;       * * * * * * * *         56 57 58 59 24 25 26 27
 728        ;;       - - - - - - - -         52 53 54 55 20 21 22 23
 729        ;;       - - - - - - - -         60 61 62 63 28 29 30 31
 730        ;;       - - - - - - - -         32 33 34 35  0  1  2  3
 731        ;;       - - - - - - - -         40 41 42 43  8  9 10 11
 732        ;;       # # # # # # # #         36 37 38 39  4  5  6  7
 733        ;;       # # # # # # # #         44 45 46 47 12 13 14 15
 734        ;;
 735        ;;       - - * * - - * *         48 49 32 33 16 17  0  1
 736        ;;       - - * * - - * *         56 57 40 41 24 25  8  9
 737        ;;       - - * * - - * *         52 53 36 37 20 21  4  5
 738        ;;       - - * * - - * *         60 61 44 45 28 29 12 13
 739        ;;       # # - - # # - -         50 51 34 35 18 19  2  3
 740        ;;       # # - - # # - -         58 59 42 43 26 27 10 11
 741        ;;       # # - - # # - -         54 55 38 39 22 23  6  7
 742        ;;       # # - - # # - -         62 63 46 47 30 31 14 15
 743        ;;
 744        ;;       - - - - - - - -         48 49 32 33 16 17  0  1
 745        ;;       * * * * * * * *         50 51 34 35 18 19  2  3
 746        ;;       - - - - - - - -         52 53 36 37 20 21  4  5
 747        ;;       * * * * * * * *         54 55 38 39 22 23  6  7
 748        ;;       # # # # # # # #         56 57 40 41 24 25  8  9
 749        ;;       - - - - - - - -         58 59 42 43 26 27 10 11
 750        ;;       # # # # # # # #         60 61 44 45 28 29 12 13
 751        ;;       - - - - - - - -         62 63 46 47 30 31 14 15
 752        ;;
 753        ;;       * - * - * - * -         57 49 41 33 25 17  9  1
 754        ;;       * - * - * - * -         59 51 43 35 27 19 11  3
 755        ;;       * - * - * - * -         61 53 45 37 29 21 13  5
 756        ;;       * - * - * - * -         63 55 47 39 31 23 15  7
 757        ;;       - # - # - # - #         56 48 40 32 24 16  8  0
 758        ;;       - # - # - # - #         58 50 42 34 26 18 10  2
 759        ;;       - # - # - # - #         60 52 44 36 28 20 12  4
 760        ;;       - # - # - # - #         62 54 46 38 30 22 14  6
 761        (new-network
 762         (make-permutation-network
 763          64                             ;  5  4  3  2  1  0
 764          '((:exchange-invert 2 5)       ; ~2  4  3 ~5  1  0
 765            (:exchange-invert 4 5)       ; ~4  2  3 ~5  1  0
 766            (:exchange        1 5)       ;  1  2  3 ~5 ~4  0
 767            (:exchange        3 5)       ;  3  2  1 ~5 ~4  0
 768            (:exchange-invert 0 5)))))   ; ~0  2  1 ~5 ~4 ~3
 769
 770   (fresh-line)
 771
 772   (let ((benes-network (benes-search fixed-ip)))
 773     (print-permutation-network benes-network)
 774     (demonstrate-permutation-network 64 benes-network :reference fixed-ip))
 775   (terpri)
 776   (print-permutation-network trad-network)
 777   (demonstrate-permutation-network 64 trad-network :reference fixed-ip)
 778   (terpri)
 779   (print-permutation-network new-network)
 780   (demonstrate-permutation-network 64 new-network :reference fixed-ip))
 781
 782 #+example
 783 (benes-search-des #( 0  0  0  0
 784                     57 49 41 33 25 17  9  1
 785                     58 50 42 34 26 18 10  2
 786                     59 51 43 35 27 19 11  3
 787                     60 52 44 36
 788                      0  0  0  0
 789                     63 55 47 39 31 23 15  7
 790                     62 54 46 38 30 22 14  6
 791                     61 53 45 37 29 21 13  5
 792                                 28 20 12  4))
 793
 794 #+example
 795 (let ((pc2 (make-array '(8 6)
 796                        :element-type 'fixnum
 797                        :initial-contents '((14 17 11 24  1  5)
 798                                            ( 3 28 15  6 21 10)
 799                                            (23 19 12  4 26  8)
 800                                            (16  7 27 20 13  2)
 801                                            (41 52 31 37 47 55)
 802                                            (30 40 51 45 33 48)
 803                                            (44 49 39 56 34 53)
 804                                            (46 42 50 36 29 32)))))
 805   (benes-search-des
 806    (make-array 64
 807                :element-type 'fixnum
 808                :initial-contents
 809                  (loop for i in '(2 4 6 8 1 3 5 7)
 810                        nconc (list 0 0)
 811                        nconc (loop for j below 6
 812                                    for x = (aref pc2 (1- i) j)
 813                                    collect (if (<= x 32) (+ x 4) (+ x 8)))))
 814    1000))