symm/des.c: Replace PC1 and PC2 permutation tables with Beneš networks.
authorMark Wooding <mdw@distorted.org.uk>
Thu, 1 Feb 2024 14:37:27 +0000 (14:37 +0000)
committerMark Wooding <mdw@distorted.org.uk>
Thu, 1 Feb 2024 14:37:27 +0000 (14:37 +0000)
base/permute.h
symm/des.c
utils/permute.lisp

index ef8eb16..21726a4 100644 (file)
  * bits within two different words.  This can be seen as a multiprecision
  * variant of the swizzle.
  *
+ * Finally, we consider general permutations.  These can be implemented using
+ * Beneš networks.  Pick some index bit number %$i$%.  By applying a swizzle
+ * with a shift by %$2^i$% to the inputs, and another to the outputs, we can
+ * reduce the problem to finding two independent permutations, one affecting
+ * bits whose index has bit %$i$% clear, and the other affecting bits whose
+ * index has bit %$i$% set.  This doesn't sound so helpful, except that (a)
+ * the smaller permutations can each be implemented in the same way, and (b)
+ * they can be performed in parallel.  Small Beneš networks can be
+ * constructed by hand, but computer assistance is useful for larger ones;
+ * there are some utilities in `utils/benes.lisp'.
+ *
  * The machinery here expects some parameters to have been defined:
  *
  *   * @regty@ should be an unsigned integer type, and
index e751cf0..19ff667 100644 (file)
@@ -46,61 +46,6 @@ const octet des_keysz[] = { KSZ_SET, 7, 8, 0 };
 
 /*----- Main code ---------------------------------------------------------*/
 
-/* --- @permute@ --- *
- *
- * Arguments:  @const char *p@ = pointer to permutation table
- *             @uint32 a, b@ = source value to permute
- *             @uint32 *d@ = destination for value
- *
- * Returns:    ---
- *
- * Use:                Performs a 64-bit permutation.  The table is given in the
- *             normal (but bizarre) DES bit numbering system.  That's not to
- *             say that the tables in this source file are like the normal
- *             DES tables, because they're not.
- */
-
-static void permute(const char *p, uint32 a, uint32 b, uint32 *d)
-{
-  uint32 x = 0, y = 0;
-  int i;
-
-  for (i = 0; i < 32; i++) {
-    int q = p[i];
-    uint32 t;
-    if (!q)
-      continue;
-    else if (q <= 32)
-      t = a;
-    else {
-      t = b;
-      q -= 32;
-    }
-    if (t & (1 << (32 - q)))
-      x |= (1 << (31 - i));
-  }
-
-  p += 32;
-
-  for (i = 0; i < 32; i++) {
-    int q = p[i];
-    uint32 t;
-    if (!q)
-      continue;
-    else if (q <= 32)
-      t = a;
-    else {
-      t = b;
-      q -= 32;
-    }
-    if (t & (1 << (32 - q)))
-      y |= (1 << (31 - i));
-  }
-
-  d[0] = x;
-  d[1] = y;
-}
-
 /* --- @des_expand@ --- *
  *
  * Arguments:  @const octet *k@ = pointer to key material
@@ -155,64 +100,19 @@ void des_expand(const octet *k, size_t n, uint32 *xx, uint32 *yy)
 
 void des_init(des_ctx *k, const void *buf, size_t sz)
 {
-  uint32 x, y;
+#define REGWD 32
+  typedef uint32 regty;
+
+  uint32 x, y, u, v;
   uint32 *kp = k->k;
-  uint32 ka[2];
   int i;
 
-  /* --- @pc1@ --- *
-   *
-   * This cryptographically useless permutation is used to mangle the key
-   * before it's subjected to the key schedule proper.  I've not actually
-   * messed it about much except for inserting padding at the beginning of
-   * the two halves of the key.
-   */
-
-  static const char pc1[] = {
-     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
-  };
-
-  /* --- @pc2@ --- *
-   *
-   * This irritating but necessary permutation mangles the key between the
-   * simple rotation-based schedule and the actual XOR with which it modifies
-   * the behaviour of the cipher.
-   *
-   * This version of the table doesn't look much like the original.  This is
-   * because some parts of the world have been permuted in order to make
-   * things simpler for the round function.  In particular, everything is
-   * rotated left one place to avoid problems with the wraparound of the
-   * expansion permutation, and the key is split between odd and even S-boxes
-   * rather than high and low ones.  That's without the complication of the
-   * padding bits in the representation of the 56-bit proto-key.
-   */
-
-  static const char pc2[] = {
-     0,         0,  3 + 4, 28 + 4, 15 + 4,  6 + 4, 21 + 4, 10 + 4, /* S-box 2 */
-     0,         0, 16 + 4,  7 + 4, 27 + 4, 20 + 4, 13 + 4,  2 + 4, /* S-box 4 */
-     0,         0, 30 + 8, 40 + 8, 51 + 8, 45 + 8, 33 + 8, 48 + 8, /* S-box 6 */
-     0,         0, 46 + 8, 42 + 8, 50 + 8, 36 + 8, 29 + 8, 32 + 8, /* S-box 8 */
-     0,         0, 14 + 4, 17 + 4, 11 + 4, 24 + 4,  1 + 4,  5 + 4, /* S-box 1 */
-     0,         0, 23 + 4, 19 + 4, 12 + 4,  4 + 4, 26 + 4,  8 + 4, /* S-box 3 */
-     0,         0, 41 + 8, 52 + 8, 31 + 8, 37 + 8, 47 + 8, 55 + 8, /* S-box 5 */
-     0,         0, 44 + 8, 49 + 8, 39 + 8, 56 + 8, 34 + 8, 53 + 8  /* S-box 7 */
-  };
-
-  /* --- @v@ --- *
+  /* --- @r@ --- *
    *
    * Contains the rotation amounts for the key halves.
    */
 
-  static const char v[] = {
+  static const char r[] = {
     1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1
   };
 
@@ -226,24 +126,106 @@ void des_init(des_ctx *k, const void *buf, size_t sz)
   KSZ_ASSERT(des, sz);
   des_expand(buf, sz, &x, &y);
 
-  /* --- Permute using the pointless PC1 --- */
+  /* --- Permute using the pointless PC1 --- *
+   *
+   * For reference, the original PC1 permutation is
+   *
+   *   Left half       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
+   *
+   *   Right half      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
+   *
+   * The network below implements this pretty directly; the two 28-bit halves
+   * end up in the least significant bits of the two output words; the parity
+   * bits, which are formally discarded, end up in the top 4 bits of each
+   * half in some random order, and are finally masked off so that they don't
+   * interfere with the rotation below.  (I did an exhaustive search, and
+   * there are no better Beneš networks.)
+   */
 
-  permute(pc1, x, y, ka);
-  x = ka[0]; y = ka[1];
+  SWIZZLE_2(x, y,  1, 0x55005401, 0x55005500);
+  SWIZZLE_2(x, y,  2, 0x32320101, 0x33330000);
+  TWIZZLE_0(x, y,     0xf0e1f0f1);
+  SWIZZLE_2(x, y,  4, 0x0f0e0f0f, 0x08050201);
+  SWIZZLE_2(x, y,  8, 0x005a003a, 0x005a005a);
+  SWIZZLE_2(x, y, 16, 0x00007c6c, 0x000023cc);
+  TWIZZLE_0(x, y,     0x20f1e0f0);
+  SWIZZLE_2(x, y,  2, 0x10000333, 0x33201013);
+  SWIZZLE_2(x, y,  1, 0x10055005, 0x10455005);
+  x &= 0x0fffffff; y &= 0x0fffffff;
 
   /* --- Now for the key schedule proper --- */
 
   for (i = 0; i < 16; i++) {
-    if (v[i] == 1) {
+    if (r[i] == 1) {
       x = ((x << 1) | (x >> 27)) & 0x0fffffff;
       y = ((y << 1) | (y >> 27)) & 0x0fffffff;
     } else {
       x = ((x << 2) | (x >> 26)) & 0x0fffffff;
       y = ((y << 2) | (y >> 26)) & 0x0fffffff;
     }
-    permute(pc2, x, y, kp);
-    kp += 2;
+
+    /* --- Apply PC2, which is another Beneš network --- *
+     *
+     * The original permutation is described as follows.
+     *
+     *         S-box 1: 14 17 11 24  1  5
+     *         S-box 2:  3 28 15  6 21 10
+     *         S-box 3: 23 19 12  4 26  8
+     *         S-box 4: 16  7 27 20 13  2
+     *         S-box 5: 41 52 31 37 47 55
+     *         S-box 6: 30 40 51 45 33 48
+     *         S-box 7: 44 49 39 56 34 53
+     *         S-box 8: 46 42 50 36 29 32
+     *
+     * Firstly, note the way that the key bits are arranged in the words @x@
+     * and @y@: viewed from the way DES numbers bits from the most-
+     * significant end down, there are four padding bits in positions 1--4,
+     * and another four in positions 33--36.  Because the bits in the left-
+     * hand half of the key all feed into the first four S-boxes, we must
+     * adjust the bit positions by 4; and we must adjust the positions of the
+     * bits in the right-hand half by 8.
+     *
+     * Secondly, this isn't how we want to apply the key.  The formal
+     * description of DES includes an `expansion' %$E$%: essentially, we take
+     * each chunk of four bits in the 32-bit half block, and glue on the
+     * nearest bits from the preceding and following chunk to make a six-bit
+     * chunk, which we then XOR with six bits of key and feed into the S-box
+     * to collapse back down to four bits.  We avoid having to do this in
+     * practice by doing the S-boxes in two steps: first, the even-numbered
+     * ones and then the odd-numbered ones.  Because these two collections of
+     * S-boxes don't involve overlapping input bits, we can just XOR in the
+     * correct key bits and apply the substitution.  There's one more little
+     * problem, which is that the input to the final S-box needs the topmost
+     * bit of the input half-block, which we handle by having previously
+     * rotated the message block left by one position.  And that's the
+     * permutation that we implement here.
+     *
+     * There are too many blank spaces to search exhaustively for an optimal
+     * network.  Based on my experience with PC1, I don't expect the optimal
+     * network to be significantly better than this one.
+     */
+
+    u = x; v = y;
+    SWIZZLE_2(u, v,  1, 0x10551050, 0x05500504);
+    SWIZZLE_2(u, v,  2, 0x12131230, 0x33102201);
+    SWIZZLE_2(u, v,  8, 0x00a200ec, 0x009100ba);
+    SWIZZLE_2(u, v, 16, 0x000012ab, 0x000028e0);
+    SWIZZLE_2(u, v,  4, 0x0a090805, 0x0b040002);
+    TWIZZLE_0(u, v,    0x33856c2a);
+    SWIZZLE_2(u, v, 16, 0x00003385, 0x00004c6a);
+    SWIZZLE_2(u, v,  8, 0x001500c8, 0x004700e8);
+    SWIZZLE_2(u, v,  2, 0x20130212, 0x00310022);
+    SWIZZLE_2(u, v,  1, 0x05404145, 0x54510510);
+    kp[0] = u; kp[1] = v; kp += 2;
   }
+
+#undef REGWD
 }
 
 /* --- @des_eblk@, @des_dblk@ --- *
index a8600fb..7c8bcba 100644 (file)
       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
 
   (fresh-line)
 
+  (let ((benes-network (benes-search fixed-ip)))
+    (print-permutation-network benes-network)
+    (demonstrate-permutation-network 64 benes-network fixed-ip))
+  (terpri)
   (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))
+
+#+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))