From 6afe4c83448d7bb36d767a352f3692f528e6783b Mon Sep 17 00:00:00 2001
From: "Juan M. Bello-Rivas" <jbellorivas@rigetti.com>
Date: Wed, 11 Sep 2019 11:20:41 -0700
Subject: [PATCH 1/6] Accelerate APPLY-QUBIT-PERMUTATION.

---
 dqvm/src/global-addresses.lisp        |   2 +-
 dqvm/src/permutation.lisp             | 296 ++++++++++++++++----------
 dqvm/tests/distributed-qvm-tests.lisp |   2 -
 dqvm/tests/permutation-tests.lisp     |  51 +++++
 4 files changed, 236 insertions(+), 115 deletions(-)

diff --git a/dqvm/src/global-addresses.lisp b/dqvm/src/global-addresses.lisp
index 1c045068..64bfccc5 100644
--- a/dqvm/src/global-addresses.lisp
+++ b/dqvm/src/global-addresses.lisp
@@ -26,7 +26,7 @@
     :reader permutation
     :writer update-permutation
     :initarg :permutation
-    :type (or null permutation)
+    :type (or null permutation permutation-transposition permutation-general)
     :documentation "Last qubit permutation evaluated, stored in a format suitable for use by APPLY-QUBIT-PERMUTATION.")
 
    ;; The following attributes are calculated during instantiation.
diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp
index 3704c70f..127d0b1a 100644
--- a/dqvm/src/permutation.lisp
+++ b/dqvm/src/permutation.lisp
@@ -4,33 +4,63 @@
 
 (in-package #:dqvm2)
 
-;;; A simple implementation of a permutation data structure.
-
-;;; Note that (make-permutation) and NIL both represent the identity.
+;;; Permutation classes for permuting sets of qubits.
+;;;
+;;; The value NIL represents the identity permutation. General permutations
+;;; are embodied by the PERMUTATION-GENERAL class. Permutations involving a
+;;; single transposition swapping 0 with another qubit are represented by the
+;;; PERMUTATION-TRANSPOSITION class.
+;;;
+;;; The generic function APPLY-QUBIT-PERMUTATION does the heavy lifting of
+;;; permuting addresses and the class hierarchy laid out here allows us to
+;;; accomplish significant speed-ups (applying a PERMUTATION-TRANSPOSITION is
+;;; more than three times faster than the equivalent application of a
+;;; PERMUTATION-GENERAL object).
 
 (deftype transposition ()
   '(or null (cons alexandria:non-negative-fixnum
              alexandria:non-negative-fixnum)))
 
 (defclass permutation ()
+  ()
+  (:documentation "Base class for permutations."))
+
+(defclass permutation-general (permutation)
   ((number-of-transpositions
-    :initarg :number-of-transpositions
     :type alexandria:non-negative-integer
+    :initarg :number-of-transpositions
     :documentation "Number of transpositions defining the permutation.")
    (transpositions
-    :initarg :transpositions
     :type list
+    :initarg :transpositions
     :reader permutation-transpositions
     :documentation "Bijective map determined by transpositions, stored as an association list sorted by CAR."))
   (:default-initargs
    :transpositions nil)
-  (:documentation "Permutation acting on sets of qubit indices."))
+  (:documentation "Arbitrary permutation acting on sets of qubit indices."))
+
+(defclass permutation-transposition ()
+  ((tau
+    :type (unsigned-byte 6)            ; Implies a maximum of 2⁶ = 64 qubits.
+    :initarg :tau
+    :initform (error-missing-initform :tau)
+    :documentation "Positive value of τ in π = (0 τ)."))
+  (:documentation "Specialized permutation involving a single transposition of the form π = (0 τ) where τ ≠ 0."))
+
+(defmethod permutation-transpositions ((permutation permutation-transposition))
+  (let ((tau (slot-value permutation 'tau)))
+    (list (cons 0 tau) (cons tau 0))))
 
-(defmethod print-object ((permutation permutation) stream)
+(defmethod print-object ((permutation permutation-general) stream)
   (print-unreadable-object (permutation stream :type t :identity t)
     (let ((transpositions (permutation-transpositions permutation)))
       (format stream "~:[~:A~;~{~A~^ ~}~]" transpositions transpositions))))
 
+(defmethod print-object ((permutation permutation-transposition) stream)
+  (print-unreadable-object (permutation stream :type t :identity t)
+    (let ((tau (slot-value permutation 'tau)))
+      (format stream "(0 . ~D) (~D . 0)" tau tau))))
+
 (defun-inlinable make-permutation (&optional transpositions)
   "Allocate a permutation defined by TRANSPOSITIONS.
 
@@ -43,11 +73,10 @@ DQVM2> (make-permutation '((2 . 1) (1 . 0)))
 #<permutation (0 . 2) (1 . 0) (2 . 1) {10086BB8B3}>
 
 Note that in the example above, the transposition (0 2) was automatically added."
-  (declare (optimize (speed 3) (safety 0))
+  (declare #.qvm::*optimize-dangerously-fast*
            (type list transpositions))
 
-  (let ((permutation (make-instance 'permutation))
-        (transpositions* nil)
+  (let ((transpositions* nil)
         (domain nil)
         (codomain nil))
 
@@ -59,6 +88,8 @@ Note that in the example above, the transposition (0 2) was automatically added.
                  (error "Malformed permutation. A mapping ~D ↦ ~D already existed."
                         (first z) (rest z))))))
 
+      (declare (inline check-transposition))
+
       (loop :for (a . b) :in transpositions :do
         (check-transposition a b)
         (unless (= a b)
@@ -69,93 +100,85 @@ Note that in the example above, the transposition (0 2) was automatically added.
     (loop :for a :of-type alexandria:non-negative-fixnum
             :in (set-difference codomain domain)
           :for b :of-type alexandria:non-negative-fixnum
-            :in (nset-difference domain codomain)
+            :in (set-difference domain codomain)
           :unless (= a b) :do
             (pushnew (cons a b) transpositions* :test #'equal))
 
-    (setf (slot-value permutation 'number-of-transpositions) (length transpositions*)
-          (slot-value permutation 'transpositions) (sort transpositions* #'< :key #'first))
-
-    permutation))
-
-(defun-inlinable inverse-permutation (permutation)
-  "Return the inverse of PERMUTATION."
-  (declare (optimize (speed 3) (safety 0)))
-  (when permutation
-
-    (let ((inverse-permutation (make-instance 'permutation))
-          (transpositions (permutation-transpositions permutation)))
-
-      (setf (slot-value inverse-permutation 'transpositions) (loop :for (a . b) :in transpositions :collect (cons b a))
-            (slot-value inverse-permutation 'number-of-transpositions) (slot-value permutation 'number-of-transpositions))
-
-      inverse-permutation)))
-
-(defun is-identity-permutation-p (permutation)
-  "Return T if PERMUTATION is the identity, NIL otherwise."
-  (if (or (null permutation) (null (permutation-transpositions permutation)))
-      t
-      nil))
-
-(defun-inlinable apply-permutation (permutation item)
-  "Apply PERMUTATION to ITEM.
-
-Examples
---------
-
-DQVM2> (apply-permutation (make-permutation) 42)
-42
-
-DQVM2> (apply-permutation (make-permutation '((2 . 0))) 2)
-0
-
-DQVM2> (apply-permutation (make-permutation '((2 . 1) (1 . 0))) 2)
-1"
-  (declare (optimize (speed 3) (safety 0))
-           (type (or null permutation) permutation)
-           (type alexandria:non-negative-fixnum item))
-  (the alexandria:non-negative-fixnum
-       (if permutation
-           (alexandria:if-let ((transposition (assoc item (permutation-transpositions permutation))))
-             (rest transposition)
-             item)
-           item)))
-
-(defun-inlinable apply-inverse-permutation (permutation item)
-  "Apply PERMUTATION⁻¹ to ITEM."
-  (apply-permutation (inverse-permutation permutation) item))
+    (cond
+      ((and (null domain) (null codomain)) nil)
+      ((and (= 1 (length domain))
+            (zerop (min (the qvm:amplitude-address (first domain))
+                        (the qvm:amplitude-address (first codomain)))))
+       (make-instance 'permutation-transposition
+                      :tau (max (the qvm:amplitude-address (first domain))
+                                (the qvm:amplitude-address (first codomain)))))
+      ((and (= 2 (length domain))
+            (null (set-difference domain codomain))
+            (zerop (the qvm:amplitude-address (apply #'min domain))))
+       (make-instance 'permutation-transposition :tau (apply #'max domain)))
+      (t
+       (make-instance 'permutation-general :number-of-transpositions (length transpositions*)
+                                           :transpositions (sort transpositions* #'< :key #'first))))))
+
+(defgeneric inverse-permutation (permutation)
+  (:documentation "Return the inverse of PERMUTATION.")
+  (declare #.qvm::*optimize-dangerously-fast*))
+
+(defmethod inverse-permutation ((permutation (eql nil)))
+  nil)
+
+(defmethod inverse-permutation ((permutation permutation-transposition))
+  permutation)
+
+(defmethod inverse-permutation ((permutation permutation-general))
+  (make-instance 'permutation-general
+                 :transpositions (loop :for (a . b) :in (permutation-transpositions permutation) :collect (cons b a))
+                 :number-of-transpositions (slot-value permutation 'number-of-transpositions)))
+
+(defgeneric is-identity-permutation-p (permutation)
+  (:documentation "Return T if PERMUTATION is the identity, NIL otherwise."))
+
+(defmethod is-identity-permutation-p ((permutation (eql nil)))
+  t)
+
+(defmethod is-identity-permutation-p ((permutation permutation-transposition))
+  nil) ; By construction PERMUTATION-TRANSPOSITION objects cannot be the identity.
+
+(defmethod is-identity-permutation-p ((permutation permutation-general))
+  (null (permutation-transpositions permutation)))
 
 (defun compose-permutations (&rest permutations)
   "Return a new permutation that is the composition of PERMUTATIONS.
 
 If PERMUTATIONS is the list π₁, π₂, ..., πₛ, then the result is the composition π₁ ∘ π₂ ∘ ... ∘ πₛ. In other words, the composition starts from right to left as in standard mathematical notation."
-  (let (transpositions)
-
-    (let (domain)
-      ;; Aggregate the domain of the composed permutation to get a list of
-      ;; all possible relevant inputs.
-      (loop :for permutation :in permutations :when permutation :do
-        (loop :for transposition :of-type transposition :in (permutation-transpositions permutation) :do
-          (let ((a (first transposition)))
-            (declare (type alexandria:non-negative-fixnum a))
-            (pushnew a domain))))
-
-      ;; Now map each domain element to obtain transpositions.
-      (loop :with codomain := (coerce domain 'vector)
-            :for permutation :in (nreverse permutations) :when permutation :do
-              (loop :for i :from 0 :for b :across codomain :do
-                (setf (aref codomain i)
-                      (apply-permutation permutation (aref codomain i))))
-            :finally
-               (loop :for a :of-type alexandria:non-negative-fixnum :in domain
-                     :for b :of-type alexandria:non-negative-fixnum :across codomain
-                     :unless (= a b) :do
-                       (pushnew (cons a b) transpositions :test #'equal))))
+  (let ((transpositions nil)
+        (domain nil))
+
+    ;; Aggregate the domain of the composed permutation to get a list of
+    ;; all possible relevant inputs.
+    (loop :for permutation :in permutations :when permutation :do
+      (loop :for transposition :of-type transposition :in (permutation-transpositions permutation) :do
+        (let ((a (first transposition)))
+          (declare (type alexandria:non-negative-fixnum a))
+          (pushnew a domain))))
+
+    ;; Now map each domain element to obtain transpositions.
+    (loop :with codomain := (coerce domain 'vector)
+          :for permutation :in (nreverse permutations) :when permutation :do
+            (loop :for i :from 0 :for b :across codomain :do
+              (setf (aref codomain i)
+                    (apply-permutation permutation (aref codomain i))))
+          :finally
+             (loop :for a :of-type alexandria:non-negative-fixnum :in domain
+                   :for b :of-type alexandria:non-negative-fixnum :across codomain
+                   :unless (= a b) :do
+                     (pushnew (cons a b) transpositions :test #'equal)))
 
     (make-permutation transpositions)))
 
-(defun-inlinable apply-qubit-permutation (permutation address)
-  "Apply PERMUTATION to an index ADDRESS within a wavefunction.
+(defgeneric apply-qubit-permutation (permutation address)
+  (:documentation
+   "Apply PERMUTATION to an index ADDRESS within a wavefunction.
 
 Examples
 --------
@@ -165,36 +188,56 @@ DQVM2> (apply-qubit-permutation (make-permutation '((2 . 0))) #b100)
 
 DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :base 2)
 100
-4"
+4")
+  (declare #.qvm::*optimize-dangerously-fast*))
+
+(defmethod apply-qubit-permutation ((permutation (eql nil)) address)
+  address)
+
+(defmethod apply-qubit-permutation ((permutation permutation-transposition) address)
+  (declare #.qvm::*optimize-dangerously-fast*
+           (type (or null permutation) permutation)
+           ;; (type qvm:amplitude-address address)
+           (type (unsigned-byte 64) address) ; Imposed maximum number of qubits.
+           (values qvm:amplitude-address))
+
+  (let ((tau (slot-value permutation 'tau)))
+    (declare (type (unsigned-byte 6) tau))
+
+    (rotatef (ldb (byte 1 0) address) (ldb (byte 1 tau) address))
+    address))
+
+(defmethod apply-qubit-permutation ((permutation permutation) address)
   ;; Alternatively, in-place permutations could be implemented following:
   ;;
   ;; F. Fich, J. Munro, and P. Poblete, “Permuting in Place,” SIAM
   ;; J. Comput., vol. 24, no. 2, pp. 266–278, Apr. 1995.
 
-  (declare (optimize (speed 3) (safety 0))
+  (declare #.qvm::*optimize-dangerously-fast*
            (type (or null permutation) permutation)
-           (type qvm:amplitude-address address))
-
-  (the qvm:amplitude-address
-       (if permutation
-           (let* ((transpositions (slot-value permutation 'transpositions))
-                  (number-of-transpositions (slot-value permutation 'number-of-transpositions))
-                  (bit-vector (make-array number-of-transpositions :element-type 'bit)))
-             ;; (declare (dynamic-extent bit-vector))
-
-             (loop :for index :from 0
-                   :for transposition :in transpositions :do
-                     (setf (bit bit-vector index) (ldb (byte 1 (first transposition))
-                                                       address)))
-
-             (loop :for index :from 0
-                   :for transposition :of-type transposition :in transpositions :do
-                     (setf address (dpb (bit bit-vector index)
-                                        ;; (byte 1 (the (unsigned-byte 6) (rest transposition))) ; Enable this for speed (assumes a maximum of 64 qubits).
-                                        (byte 1 (rest transposition))
-                                        address))
-                   :finally (return address)))
-           address)))
+           ;; (type qvm:amplitude-address address)
+           (type (unsigned-byte 64) address) ; Imposed maximum number of qubits.
+           (values qvm:amplitude-address))
+
+  (let* ((transpositions (slot-value permutation 'transpositions))
+         (number-of-transpositions (slot-value permutation 'number-of-transpositions))
+         (bit-vector (make-array number-of-transpositions :element-type 'bit)))
+    (declare (type (integer 0 128) number-of-transpositions)
+             (dynamic-extent bit-vector))
+
+    (loop :for index :from 0
+          :for transposition :in transpositions :do
+            (setf (bit bit-vector index) (ldb (byte 1 (first transposition))
+                                              address)))
+
+    (loop :for index :from 0
+          :for transposition :of-type transposition :in transpositions :do
+            (setf address (the qvm:amplitude-address
+                               (dpb (bit bit-vector index)
+                                    (byte 1 (the (unsigned-byte 6) (rest transposition))) ; Enable this for speed (assumes a maximum of 64 qubits).
+                                    ;; (byte 1 (rest transposition))
+                                    address)))
+          :finally (return address))))
 
 (defun-inlinable apply-inverse-qubit-permutation (permutation address)
   (apply-qubit-permutation (inverse-permutation permutation) address))
@@ -203,8 +246,8 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
                                               (stream *standard-output*))
   "Print the address permutation induced by PERMUTATION (possibly using up to NUMBER-OF-QUBITS) in STREAM."
   (let* ((n (or number-of-qubits
-                (1+ (loop :for transposition :in (permutation-transpositions permutation)
-                          :maximizing (max transposition)))))
+                (1+ (loop :for (a . b) :in (permutation-transpositions permutation)
+                          :maximizing (max a b)))))
          (max-value (expt 2 n))
          (aux-control-string (format nil "~~~DD |~~~D,'0B>"
                                      (ceiling (log max-value 10)) n))
@@ -213,3 +256,32 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
     (dotimes (i1 max-value (values))
       (let ((i2 (apply-qubit-permutation permutation i1)))
         (format stream control-string i1 i1 i2 i2)))))
+
+(defun-inlinable apply-permutation (permutation item)
+  "Apply PERMUTATION to ITEM.
+
+Examples
+--------
+
+DQVM2> (apply-permutation (make-permutation) 42)
+42
+
+DQVM2> (apply-permutation (make-permutation '((2 . 0))) 2)
+0
+
+DQVM2> (apply-permutation (make-permutation '((2 . 1) (1 . 0))) 2)
+1"
+  (declare #.qvm::*optimize-dangerously-fast*
+           (type (or null permutation) permutation)
+           (type alexandria:non-negative-fixnum item)
+           (values alexandria:non-negative-fixnum))
+
+  (if permutation
+      (alexandria:if-let ((transposition (assoc item (permutation-transpositions permutation))))
+        (rest transposition)
+        item)
+      item))
+
+(defun-inlinable apply-inverse-permutation (permutation item)
+  "Apply PERMUTATION⁻¹ to ITEM."
+  (apply-permutation (inverse-permutation permutation) item))
diff --git a/dqvm/tests/distributed-qvm-tests.lisp b/dqvm/tests/distributed-qvm-tests.lisp
index e56adf2c..17a7c473 100644
--- a/dqvm/tests/distributed-qvm-tests.lisp
+++ b/dqvm/tests/distributed-qvm-tests.lisp
@@ -25,8 +25,6 @@
   "Find amplitude addresses to exchange when applying NEXT-PERMUTATION and the rank where the amplitudes are located.
 
 Returns four sequences: current addresses, new addresses, and the source and target addresses."
-  (check-type next-permutation permutation)
-
   (let ((permutation (permutation addresses))
         (effective-permutation
           (dqvm2::get-effective-permutation addresses next-permutation))
diff --git a/dqvm/tests/permutation-tests.lisp b/dqvm/tests/permutation-tests.lisp
index 50bf6c99..82a1028f 100644
--- a/dqvm/tests/permutation-tests.lisp
+++ b/dqvm/tests/permutation-tests.lisp
@@ -26,6 +26,13 @@
     (is (= (apply-inverse-permutation permutation 1) 1))
     (is (= (apply-inverse-permutation permutation 2) 2)))
 
+  (let ((permutation (make-permutation '((2 . 1)))))
+    (is (eq (type-of permutation) 'dqvm2::permutation-general)))
+
+  (let ((permutation (make-permutation '((2 . 0)))))
+    (is (eq (type-of permutation) 'dqvm2::permutation-transposition))
+    (is (eq permutation (inverse-permutation permutation))))
+
   (let ((permutation (make-permutation '((2 . 1) (1 . 0)))))
     (is (not (is-identity-permutation-p permutation)))
 
@@ -72,3 +79,47 @@
     (is (= (apply-permutation composition 0) 2))
     (is (= (apply-permutation composition 1) 1))
     (is (= (apply-permutation composition 2) 0))))
+
+(deftest benchmark-apply-qubit-permutation ()
+  (labels ((get-elapsed-time-in-seconds (start stop)
+             "Compute elapsed time in seconds between START and STOP."
+             (float (/ (- stop start) internal-time-units-per-second)))
+
+           (time-apply-qubit-permutation (permutation number-of-qubits)
+             "Measure the time taken by calls to APPLY-QUBIT-PERMUTATION on addresses from 0 to 2^NUMBER-OF-QUBITS."
+             (let ((start (get-internal-real-time)))
+               (dotimes (x (expt 2 number-of-qubits))
+                 (let ((y (apply-qubit-permutation permutation x)))
+                   (values x y)))
+               (get-elapsed-time-in-seconds start (get-internal-real-time))))
+
+           ;; (time-map-reordered-amplitudes (permutation number-of-qubits)
+           ;;   "Measure the time taken by calls to MAP-REORDERED-AMPLITUDES on addresses from 0 to 2^NUMBER-OF-QUBITS."
+           ;;   (let ((nat-tuple (apply #'qvm::nat-tuple
+           ;;                           (mapcar (lambda (x)
+           ;;                                     (apply-permutation permutation x))
+           ;;                                   (loop :for i :below number-of-qubits :collect i))))
+           ;;         (start (get-internal-real-time)))
+           ;;     (qvm::map-reordered-amplitudes 0 (lambda (x y) (values x y)) nat-tuple)
+           ;;     (get-elapsed-time-in-seconds start (get-internal-real-time))))
+           )
+
+    (let* ((tau 4)
+           (number-of-qubits 24)
+           (permutation-0 (make-instance 'dqvm2::permutation-transposition :tau tau))
+           (permutation-1 (make-instance 'dqvm2::permutation-general
+                                         :number-of-transpositions 2
+                                         :transpositions (list (cons 0 tau) (cons tau 0)))))
+
+      (loop :for x :below (expt 2 (1+ tau)) :do
+        (is (= (apply-qubit-permutation permutation-0 x)
+               (apply-qubit-permutation permutation-1 x))))
+
+      (is (> (/ (time-apply-qubit-permutation permutation-1 number-of-qubits)
+                (time-apply-qubit-permutation permutation-0 number-of-qubits))
+             3))
+
+      ;; (is (> (/ (time-map-reordered-amplitudes permutation-0 number-of-qubits)
+      ;;           (time-apply-qubit-permutation permutation-0 number-of-qubits))
+      ;;        7))
+      )))

From 00c4b6277193b55edde3e0247ab3a78e4c59fc3a Mon Sep 17 00:00:00 2001
From: "Juan M. Bello-Rivas" <jbellorivas@rigetti.com>
Date: Thu, 12 Sep 2019 06:58:13 -0700
Subject: [PATCH 2/6] Generate code for applying permutations.

---
 dqvm/src/permutation.lisp | 60 +++++++++++++++++++++++++++++++++++----
 1 file changed, 55 insertions(+), 5 deletions(-)

diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp
index 127d0b1a..628754b9 100644
--- a/dqvm/src/permutation.lisp
+++ b/dqvm/src/permutation.lisp
@@ -90,7 +90,7 @@ Note that in the example above, the transposition (0 2) was automatically added.
 
       (declare (inline check-transposition))
 
-      (loop :for (a . b) :in transpositions :do
+      (loop :for (a . b) :of-type alexandria:non-negative-fixnum :in transpositions :do
         (check-transposition a b)
         (unless (= a b)
           (pushnew (cons a b) transpositions*)
@@ -204,10 +204,11 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
   (let ((tau (slot-value permutation 'tau)))
     (declare (type (unsigned-byte 6) tau))
 
-    (rotatef (ldb (byte 1 0) address) (ldb (byte 1 tau) address))
-    address))
+    ;; Swap bits 0 and TAU in ADDRESS.
+    (let ((x (logxor (logand address 1) (logand (ash address (- tau)) 1))))
+      (logxor address (logior x (ash x tau))))))
 
-(defmethod apply-qubit-permutation ((permutation permutation) address)
+(defmethod apply-qubit-permutation ((permutation permutation-general) address)
   ;; Alternatively, in-place permutations could be implemented following:
   ;;
   ;; F. Fich, J. Munro, and P. Poblete, “Permuting in Place,” SIAM
@@ -226,7 +227,7 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
              (dynamic-extent bit-vector))
 
     (loop :for index :from 0
-          :for transposition :in transpositions :do
+          :for transposition :of-type transposition :in transpositions :do
             (setf (bit bit-vector index) (ldb (byte 1 (first transposition))
                                               address)))
 
@@ -239,6 +240,55 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
                                     address)))
           :finally (return address))))
 
+(defgeneric generate-qubit-permutation-code (permutation)
+  (:documentation "Generate lambda function equivalent to APPLY-QUBIT-PERMUTATION suitable to be compiled.")
+  (declare #.qvm::*optimize-dangerously-fast*))
+
+(defmethod generate-qubit-permutation-code ((permutation (eql nil)))
+  (let ((address (gensym "ADDRESS-")))
+    `(lambda (,address)
+       (declare #.qvm::*optimize-dangerously-fast*)
+       ,address)))
+
+(defmethod generate-qubit-permutation-code ((permutation permutation-transposition))
+  (let* ((address (gensym "ADDRESS-"))
+         (tau (slot-value permutation 'tau))
+         (minus-tau (- tau)))
+    `(lambda (,address)
+       (declare #.qvm::*optimize-dangerously-fast*
+                (type (unsigned-byte 64) ,address) ; Imposed maximum number of qubits.
+                (values qvm:amplitude-address))
+
+       ;; Swap bits 0 and TAU in ADDRESS.
+       (let ((x (logxor (logand ,address 1) (logand (ash ,address ,minus-tau) 1))))
+         (logxor ,address (logior x (ash x ,tau)))))))
+
+(defmethod generate-qubit-permutation-code ((permutation permutation-general))
+  (let ((address (gensym "ADDRESS-"))
+        (transpositions (slot-value permutation 'transpositions))
+        (number-of-transpositions (slot-value permutation 'number-of-transpositions)))
+    `(lambda (,address)
+       (declare #.qvm::*optimize-dangerously-fast*
+                (type (or null permutation) permutation)
+                (type (unsigned-byte 64) ,address)
+                (values qvm:amplitude-address))
+
+       (let ((bit-vector (make-array ,number-of-transpositions :element-type 'bit)))
+         (declare (dynamic-extent bit-vector))
+
+         ,@(loop :for index :from 0
+                 :for transposition :of-type transposition :in transpositions
+                 :collect `(setf (bit bit-vector ,index) (ldb (byte 1 ,(first transposition))
+                                                              ,address)))
+
+         ,@(loop :for index :from 0
+                 :for transposition :of-type transposition :in transpositions
+                 :collect `(setf ,address (the qvm:amplitude-address
+                                               (dpb (bit bit-vector ,index)
+                                                    (byte 1 (the (unsigned-byte 6) ,(rest transposition)))
+                                                    ,address))))
+         ,address))))
+
 (defun-inlinable apply-inverse-qubit-permutation (permutation address)
   (apply-qubit-permutation (inverse-permutation permutation) address))
 

From f2619b83dcc01dfb1d6d1920084c29f9f6a24b21 Mon Sep 17 00:00:00 2001
From: "Juan M. Bello-Rivas" <jbellorivas@rigetti.com>
Date: Thu, 12 Sep 2019 09:32:38 -0700
Subject: [PATCH 3/6] Define meaningful types.

---
 dqvm/src/permutation.lisp | 58 +++++++++++++++++++++------------------
 1 file changed, 31 insertions(+), 27 deletions(-)

diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp
index 628754b9..fefa27e9 100644
--- a/dqvm/src/permutation.lisp
+++ b/dqvm/src/permutation.lisp
@@ -17,9 +17,17 @@
 ;;; more than three times faster than the equivalent application of a
 ;;; PERMUTATION-GENERAL object).
 
+(defconstant +qubit-index-length+ (ceiling (log qvm::+max-nat-tuple-cardinality+ 2))
+  "Number of bits required to represent a qubit index.")
+
+(defconstant +max-number-of-transpositions+ (* 2 #.qvm::+max-nat-tuple-cardinality+)
+  "Upper bound on the number of transpositions defining an arbitrary permutation.")
+
+(deftype qubit-index ()
+  '(unsigned-byte #.+qubit-index-length+))
+
 (deftype transposition ()
-  '(or null (cons alexandria:non-negative-fixnum
-             alexandria:non-negative-fixnum)))
+  '(or null (cons qubit-index qubit-index)))
 
 (defclass permutation ()
   ()
@@ -39,9 +47,9 @@
    :transpositions nil)
   (:documentation "Arbitrary permutation acting on sets of qubit indices."))
 
-(defclass permutation-transposition ()
+(defclass permutation-transposition (permutation)
   ((tau
-    :type (unsigned-byte 6)            ; Implies a maximum of 2⁶ = 64 qubits.
+    :type qubit-index
     :initarg :tau
     :initform (error-missing-initform :tau)
     :documentation "Positive value of τ in π = (0 τ)."))
@@ -81,7 +89,7 @@ Note that in the example above, the transposition (0 2) was automatically added.
         (codomain nil))
 
     (flet ((check-transposition (a b)
-             (declare (type alexandria:non-negative-fixnum a b))
+             (declare (type qubit-index a b))
              (let ((x (assoc a transpositions*))
                    (y (rassoc b transpositions*)))
                (alexandria:when-let ((z (or x y)))
@@ -90,17 +98,15 @@ Note that in the example above, the transposition (0 2) was automatically added.
 
       (declare (inline check-transposition))
 
-      (loop :for (a . b) :of-type alexandria:non-negative-fixnum :in transpositions :do
+      (loop :for (a . b) :of-type qubit-index :in transpositions :do
         (check-transposition a b)
         (unless (= a b)
           (pushnew (cons a b) transpositions*)
           (pushnew a domain)
           (pushnew b codomain))))
 
-    (loop :for a :of-type alexandria:non-negative-fixnum
-            :in (set-difference codomain domain)
-          :for b :of-type alexandria:non-negative-fixnum
-            :in (set-difference domain codomain)
+    (loop :for a :of-type qubit-index :in (set-difference codomain domain)
+          :for b :of-type qubit-index :in (set-difference domain codomain)
           :unless (= a b) :do
             (pushnew (cons a b) transpositions* :test #'equal))
 
@@ -159,7 +165,7 @@ If PERMUTATIONS is the list π₁, π₂, ..., πₛ, then the result is the com
     (loop :for permutation :in permutations :when permutation :do
       (loop :for transposition :of-type transposition :in (permutation-transpositions permutation) :do
         (let ((a (first transposition)))
-          (declare (type alexandria:non-negative-fixnum a))
+          (declare (type qubit-index a))
           (pushnew a domain))))
 
     ;; Now map each domain element to obtain transpositions.
@@ -169,8 +175,8 @@ If PERMUTATIONS is the list π₁, π₂, ..., πₛ, then the result is the com
               (setf (aref codomain i)
                     (apply-permutation permutation (aref codomain i))))
           :finally
-             (loop :for a :of-type alexandria:non-negative-fixnum :in domain
-                   :for b :of-type alexandria:non-negative-fixnum :across codomain
+             (loop :for a :of-type qubit-index :in domain
+                   :for b :of-type qubit-index :across codomain
                    :unless (= a b) :do
                      (pushnew (cons a b) transpositions :test #'equal)))
 
@@ -196,13 +202,12 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
 
 (defmethod apply-qubit-permutation ((permutation permutation-transposition) address)
   (declare #.qvm::*optimize-dangerously-fast*
-           (type (or null permutation) permutation)
-           ;; (type qvm:amplitude-address address)
-           (type (unsigned-byte 64) address) ; Imposed maximum number of qubits.
+           (type permutation permutation)
+           (type qvm:amplitude-address address)
            (values qvm:amplitude-address))
 
   (let ((tau (slot-value permutation 'tau)))
-    (declare (type (unsigned-byte 6) tau))
+    (declare (type qubit-index tau))
 
     ;; Swap bits 0 and TAU in ADDRESS.
     (let ((x (logxor (logand address 1) (logand (ash address (- tau)) 1))))
@@ -215,15 +220,14 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
   ;; J. Comput., vol. 24, no. 2, pp. 266–278, Apr. 1995.
 
   (declare #.qvm::*optimize-dangerously-fast*
-           (type (or null permutation) permutation)
-           ;; (type qvm:amplitude-address address)
-           (type (unsigned-byte 64) address) ; Imposed maximum number of qubits.
+           (type permutation permutation)
+           (type qvm:amplitude-address address)
            (values qvm:amplitude-address))
 
   (let* ((transpositions (slot-value permutation 'transpositions))
          (number-of-transpositions (slot-value permutation 'number-of-transpositions))
          (bit-vector (make-array number-of-transpositions :element-type 'bit)))
-    (declare (type (integer 0 128) number-of-transpositions)
+    (declare (type (integer 0 #.+max-number-of-transpositions+) number-of-transpositions)
              (dynamic-extent bit-vector))
 
     (loop :for index :from 0
@@ -235,8 +239,7 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
           :for transposition :of-type transposition :in transpositions :do
             (setf address (the qvm:amplitude-address
                                (dpb (bit bit-vector index)
-                                    (byte 1 (the (unsigned-byte 6) (rest transposition))) ; Enable this for speed (assumes a maximum of 64 qubits).
-                                    ;; (byte 1 (rest transposition))
+                                    (byte 1 (the qubit-index (rest transposition)))
                                     address)))
           :finally (return address))))
 
@@ -256,7 +259,8 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
          (minus-tau (- tau)))
     `(lambda (,address)
        (declare #.qvm::*optimize-dangerously-fast*
-                (type (unsigned-byte 64) ,address) ; Imposed maximum number of qubits.
+                (type permutation permutation)
+                (type qvm:amplitude-address ,address)
                 (values qvm:amplitude-address))
 
        ;; Swap bits 0 and TAU in ADDRESS.
@@ -269,8 +273,8 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
         (number-of-transpositions (slot-value permutation 'number-of-transpositions)))
     `(lambda (,address)
        (declare #.qvm::*optimize-dangerously-fast*
-                (type (or null permutation) permutation)
-                (type (unsigned-byte 64) ,address)
+                (type permutation permutation)
+                (type qvm:amplitude-address ,address)
                 (values qvm:amplitude-address))
 
        (let ((bit-vector (make-array ,number-of-transpositions :element-type 'bit)))
@@ -285,7 +289,7 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
                  :for transposition :of-type transposition :in transpositions
                  :collect `(setf ,address (the qvm:amplitude-address
                                                (dpb (bit bit-vector ,index)
-                                                    (byte 1 (the (unsigned-byte 6) ,(rest transposition)))
+                                                    (byte 1 (the qubit-index ,(rest transposition)))
                                                     ,address))))
          ,address))))
 

From a2b7d2678c0cb9605183612cd915c487c06c23cc Mon Sep 17 00:00:00 2001
From: "Juan M. Bello-Rivas" <jbellorivas@rigetti.com>
Date: Thu, 12 Sep 2019 10:52:40 -0700
Subject: [PATCH 4/6] Remove spurious type declaration.

---
 dqvm/src/permutation.lisp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp
index fefa27e9..cdf7ae7f 100644
--- a/dqvm/src/permutation.lisp
+++ b/dqvm/src/permutation.lisp
@@ -259,7 +259,6 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
          (minus-tau (- tau)))
     `(lambda (,address)
        (declare #.qvm::*optimize-dangerously-fast*
-                (type permutation permutation)
                 (type qvm:amplitude-address ,address)
                 (values qvm:amplitude-address))
 

From 19502eaa7991a49290df430d5763e853d4ade70f Mon Sep 17 00:00:00 2001
From: "Juan M. Bello-Rivas" <jbellorivas@rigetti.com>
Date: Thu, 12 Sep 2019 13:16:10 -0700
Subject: [PATCH 5/6] Address reviewers' remarks.

---
 dqvm/src/package.lisp             |  2 +-
 dqvm/src/permutation.lisp         | 46 +++++++++++++++----------------
 dqvm/tests/permutation-tests.lisp | 14 +++++-----
 3 files changed, 30 insertions(+), 32 deletions(-)

diff --git a/dqvm/src/package.lisp b/dqvm/src/package.lisp
index ff30e3b7..561314ba 100644
--- a/dqvm/src/package.lisp
+++ b/dqvm/src/package.lisp
@@ -38,7 +38,7 @@
            #:global-addresses
            #:global-addresses=
            #:inverse-permutation
-           #:is-identity-permutation-p
+           #:identity-permutation-p
            #:make-addresses
            #:make-addresses-like
            #:make-distributed-qvm
diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp
index cdf7ae7f..7fd1983d 100644
--- a/dqvm/src/permutation.lisp
+++ b/dqvm/src/permutation.lisp
@@ -29,11 +29,7 @@
 (deftype transposition ()
   '(or null (cons qubit-index qubit-index)))
 
-(defclass permutation ()
-  ()
-  (:documentation "Base class for permutations."))
-
-(defclass permutation-general (permutation)
+(defclass permutation-general ()
   ((number-of-transpositions
     :type alexandria:non-negative-integer
     :initarg :number-of-transpositions
@@ -47,7 +43,7 @@
    :transpositions nil)
   (:documentation "Arbitrary permutation acting on sets of qubit indices."))
 
-(defclass permutation-transposition (permutation)
+(defclass permutation-transposition ()
   ((tau
     :type qubit-index
     :initarg :tau
@@ -55,6 +51,9 @@
     :documentation "Positive value of τ in π = (0 τ)."))
   (:documentation "Specialized permutation involving a single transposition of the form π = (0 τ) where τ ≠ 0."))
 
+(deftype permutation ()
+  '(or null permutation-general permutation-transposition))
+
 (defmethod permutation-transpositions ((permutation permutation-transposition))
   (let ((tau (slot-value permutation 'tau)))
     (list (cons 0 tau) (cons tau 0))))
@@ -120,8 +119,8 @@ Note that in the example above, the transposition (0 2) was automatically added.
                                 (the qvm:amplitude-address (first codomain)))))
       ((and (= 2 (length domain))
             (null (set-difference domain codomain))
-            (zerop (the qvm:amplitude-address (apply #'min domain))))
-       (make-instance 'permutation-transposition :tau (apply #'max domain)))
+            (zerop (the qvm:amplitude-address (alexandria:extremum domain #'<))))
+       (make-instance 'permutation-transposition :tau (alexandria:extremum domain #'>)))
       (t
        (make-instance 'permutation-general :number-of-transpositions (length transpositions*)
                                            :transpositions (sort transpositions* #'< :key #'first))))))
@@ -130,7 +129,7 @@ Note that in the example above, the transposition (0 2) was automatically added.
   (:documentation "Return the inverse of PERMUTATION.")
   (declare #.qvm::*optimize-dangerously-fast*))
 
-(defmethod inverse-permutation ((permutation (eql nil)))
+(defmethod inverse-permutation ((permutation null))
   nil)
 
 (defmethod inverse-permutation ((permutation permutation-transposition))
@@ -141,16 +140,16 @@ Note that in the example above, the transposition (0 2) was automatically added.
                  :transpositions (loop :for (a . b) :in (permutation-transpositions permutation) :collect (cons b a))
                  :number-of-transpositions (slot-value permutation 'number-of-transpositions)))
 
-(defgeneric is-identity-permutation-p (permutation)
+(defgeneric identity-permutation-p (permutation)
   (:documentation "Return T if PERMUTATION is the identity, NIL otherwise."))
 
-(defmethod is-identity-permutation-p ((permutation (eql nil)))
+(defmethod identity-permutation-p ((permutation null))
   t)
 
-(defmethod is-identity-permutation-p ((permutation permutation-transposition))
+(defmethod identity-permutation-p ((permutation permutation-transposition))
   nil) ; By construction PERMUTATION-TRANSPOSITION objects cannot be the identity.
 
-(defmethod is-identity-permutation-p ((permutation permutation-general))
+(defmethod identity-permutation-p ((permutation permutation-general))
   (null (permutation-transpositions permutation)))
 
 (defun compose-permutations (&rest permutations)
@@ -197,12 +196,12 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
 4")
   (declare #.qvm::*optimize-dangerously-fast*))
 
-(defmethod apply-qubit-permutation ((permutation (eql nil)) address)
+(defmethod apply-qubit-permutation ((permutation null) address)
   address)
 
 (defmethod apply-qubit-permutation ((permutation permutation-transposition) address)
   (declare #.qvm::*optimize-dangerously-fast*
-           (type permutation permutation)
+           (type permutation-transposition)
            (type qvm:amplitude-address address)
            (values qvm:amplitude-address))
 
@@ -220,7 +219,7 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
   ;; J. Comput., vol. 24, no. 2, pp. 266–278, Apr. 1995.
 
   (declare #.qvm::*optimize-dangerously-fast*
-           (type permutation permutation)
+           (type permutation-general permutation)
            (type qvm:amplitude-address address)
            (values qvm:amplitude-address))
 
@@ -232,13 +231,13 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
 
     (loop :for index :from 0
           :for transposition :of-type transposition :in transpositions :do
-            (setf (bit bit-vector index) (ldb (byte 1 (first transposition))
-                                              address)))
+            (setf (sbit bit-vector index) (ldb (byte 1 (first transposition))
+                                               address)))
 
     (loop :for index :from 0
           :for transposition :of-type transposition :in transpositions :do
             (setf address (the qvm:amplitude-address
-                               (dpb (bit bit-vector index)
+                               (dpb (sbit bit-vector index)
                                     (byte 1 (the qubit-index (rest transposition)))
                                     address)))
           :finally (return address))))
@@ -247,7 +246,7 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
   (:documentation "Generate lambda function equivalent to APPLY-QUBIT-PERMUTATION suitable to be compiled.")
   (declare #.qvm::*optimize-dangerously-fast*))
 
-(defmethod generate-qubit-permutation-code ((permutation (eql nil)))
+(defmethod generate-qubit-permutation-code ((permutation null))
   (let ((address (gensym "ADDRESS-")))
     `(lambda (,address)
        (declare #.qvm::*optimize-dangerously-fast*)
@@ -272,7 +271,6 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
         (number-of-transpositions (slot-value permutation 'number-of-transpositions)))
     `(lambda (,address)
        (declare #.qvm::*optimize-dangerously-fast*
-                (type permutation permutation)
                 (type qvm:amplitude-address ,address)
                 (values qvm:amplitude-address))
 
@@ -281,13 +279,13 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
 
          ,@(loop :for index :from 0
                  :for transposition :of-type transposition :in transpositions
-                 :collect `(setf (bit bit-vector ,index) (ldb (byte 1 ,(first transposition))
-                                                              ,address)))
+                 :collect `(setf (sbit bit-vector ,index) (ldb (byte 1 ,(first transposition))
+                                                               ,address)))
 
          ,@(loop :for index :from 0
                  :for transposition :of-type transposition :in transpositions
                  :collect `(setf ,address (the qvm:amplitude-address
-                                               (dpb (bit bit-vector ,index)
+                                               (dpb (sbit bit-vector ,index)
                                                     (byte 1 (the qubit-index ,(rest transposition)))
                                                     ,address))))
          ,address))))
diff --git a/dqvm/tests/permutation-tests.lisp b/dqvm/tests/permutation-tests.lisp
index 82a1028f..55d0c042 100644
--- a/dqvm/tests/permutation-tests.lisp
+++ b/dqvm/tests/permutation-tests.lisp
@@ -9,14 +9,14 @@
   (signals error (make-permutation '((1 . 2) (3 . 2))))
 
   (let ((permutation '()))
-    (is (is-identity-permutation-p permutation))
+    (is (identity-permutation-p permutation))
 
     (is (= (apply-permutation permutation 0) 0))
     (is (= (apply-permutation permutation 1) 1))
     (is (= (apply-permutation permutation 2) 2)))
 
   (let ((permutation (make-permutation)))
-    (is (is-identity-permutation-p permutation))
+    (is (identity-permutation-p permutation))
 
     (is (= (apply-permutation permutation 0) 0))
     (is (= (apply-permutation permutation 1) 1))
@@ -34,7 +34,7 @@
     (is (eq permutation (inverse-permutation permutation))))
 
   (let ((permutation (make-permutation '((2 . 1) (1 . 0)))))
-    (is (not (is-identity-permutation-p permutation)))
+    (is (not (identity-permutation-p permutation)))
 
     (is (= (apply-permutation permutation 2) 1))
     (is (= (apply-permutation permutation 1) 0))
@@ -65,14 +65,14 @@
   (let* ((permutation1 (make-permutation '((2 . 1))))
          (permutation2 (make-permutation '((0 . 1) (1 . 2))))
          (composition (compose-permutations permutation1 permutation2)))
-    (is (not (is-identity-permutation-p permutation1)))
-    (is (not (is-identity-permutation-p permutation2)))
+    (is (not (identity-permutation-p permutation1)))
+    (is (not (identity-permutation-p permutation2)))
 
-    (is (is-identity-permutation-p
+    (is (identity-permutation-p
          (compose-permutations permutation1
                                (inverse-permutation permutation1))))
 
-    (is (is-identity-permutation-p
+    (is (identity-permutation-p
          (compose-permutations permutation2
                                (inverse-permutation permutation2))))
 

From b63eed331285446ea0f54453c8c62c8946e55c61 Mon Sep 17 00:00:00 2001
From: "Juan M. Bello-Rivas" <jbellorivas@rigetti.com>
Date: Sun, 15 Sep 2019 21:56:02 -0700
Subject: [PATCH 6/6] Compile qubit permutations using look-up tables.

---
 dqvm/src/global-addresses.lisp    |  2 +-
 dqvm/src/permutation.lisp         | 87 +++++++++++++++++++++++--------
 dqvm/tests/permutation-tests.lisp | 23 +++++++-
 3 files changed, 88 insertions(+), 24 deletions(-)

diff --git a/dqvm/src/global-addresses.lisp b/dqvm/src/global-addresses.lisp
index 64bfccc5..1c045068 100644
--- a/dqvm/src/global-addresses.lisp
+++ b/dqvm/src/global-addresses.lisp
@@ -26,7 +26,7 @@
     :reader permutation
     :writer update-permutation
     :initarg :permutation
-    :type (or null permutation permutation-transposition permutation-general)
+    :type (or null permutation)
     :documentation "Last qubit permutation evaluated, stored in a format suitable for use by APPLY-QUBIT-PERMUTATION.")
 
    ;; The following attributes are calculated during instantiation.
diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp
index 7fd1983d..c71023bf 100644
--- a/dqvm/src/permutation.lisp
+++ b/dqvm/src/permutation.lisp
@@ -29,7 +29,17 @@
 (deftype transposition ()
   '(or null (cons qubit-index qubit-index)))
 
-(defclass permutation-general ()
+(defclass permutation ()
+  ((size
+    :type (integer 0 +qubit-index-length+)
+    :reader permutation-size
+    :initarg :size
+    :documentation "Maximum number of bits on which the permutation acts."))
+  (:default-initargs
+   :size 0)
+  (:documentation "Permutation of qubits."))
+
+(defclass permutation-general (permutation)
   ((number-of-transpositions
     :type alexandria:non-negative-integer
     :initarg :number-of-transpositions
@@ -40,10 +50,11 @@
     :reader permutation-transpositions
     :documentation "Bijective map determined by transpositions, stored as an association list sorted by CAR."))
   (:default-initargs
+   :number-of-transpositions 0
    :transpositions nil)
   (:documentation "Arbitrary permutation acting on sets of qubit indices."))
 
-(defclass permutation-transposition ()
+(defclass permutation-transposition (permutation)
   ((tau
     :type qubit-index
     :initarg :tau
@@ -51,9 +62,6 @@
     :documentation "Positive value of τ in π = (0 τ)."))
   (:documentation "Specialized permutation involving a single transposition of the form π = (0 τ) where τ ≠ 0."))
 
-(deftype permutation ()
-  '(or null permutation-general permutation-transposition))
-
 (defmethod permutation-transpositions ((permutation permutation-transposition))
   (let ((tau (slot-value permutation 'tau)))
     (list (cons 0 tau) (cons tau 0))))
@@ -85,7 +93,10 @@ Note that in the example above, the transposition (0 2) was automatically added.
 
   (let ((transpositions* nil)
         (domain nil)
-        (codomain nil))
+        (codomain nil)
+        (max-index 0))
+
+    (declare (type qubit-index max-index))
 
     (flet ((check-transposition (a b)
              (declare (type qubit-index a b))
@@ -97,12 +108,14 @@ Note that in the example above, the transposition (0 2) was automatically added.
 
       (declare (inline check-transposition))
 
-      (loop :for (a . b) :of-type qubit-index :in transpositions :do
-        (check-transposition a b)
-        (unless (= a b)
-          (pushnew (cons a b) transpositions*)
-          (pushnew a domain)
-          (pushnew b codomain))))
+      (loop :for (a . b) :of-type qubit-index :in transpositions
+            :maximize (max a b) :into max :do
+              (check-transposition a b)
+              (unless (= a b)
+                (pushnew (cons a b) transpositions*)
+                (pushnew a domain)
+                (pushnew b codomain))
+            :finally (setf max-index max)))
 
     (loop :for a :of-type qubit-index :in (set-difference codomain domain)
           :for b :of-type qubit-index :in (set-difference domain codomain)
@@ -114,16 +127,15 @@ Note that in the example above, the transposition (0 2) was automatically added.
       ((and (= 1 (length domain))
             (zerop (min (the qvm:amplitude-address (first domain))
                         (the qvm:amplitude-address (first codomain)))))
-       (make-instance 'permutation-transposition
-                      :tau (max (the qvm:amplitude-address (first domain))
-                                (the qvm:amplitude-address (first codomain)))))
+       (make-instance 'permutation-transposition :tau max-index :size (1+ max-index)))
       ((and (= 2 (length domain))
             (null (set-difference domain codomain))
             (zerop (the qvm:amplitude-address (alexandria:extremum domain #'<))))
-       (make-instance 'permutation-transposition :tau (alexandria:extremum domain #'>)))
+       (make-instance 'permutation-transposition :tau max-index :size (1+ max-index)))
       (t
        (make-instance 'permutation-general :number-of-transpositions (length transpositions*)
-                                           :transpositions (sort transpositions* #'< :key #'first))))))
+                                           :transpositions (sort transpositions* #'< :key #'first)
+                                           :size (1+ max-index))))))
 
 (defgeneric inverse-permutation (permutation)
   (:documentation "Return the inverse of PERMUTATION.")
@@ -249,7 +261,9 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
 (defmethod generate-qubit-permutation-code ((permutation null))
   (let ((address (gensym "ADDRESS-")))
     `(lambda (,address)
-       (declare #.qvm::*optimize-dangerously-fast*)
+       (declare #.qvm::*optimize-dangerously-fast*
+                (type qvm:amplitude-address ,address)
+                (values qvm:amplitude-address))
        ,address)))
 
 (defmethod generate-qubit-permutation-code ((permutation permutation-transposition))
@@ -290,15 +304,46 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas
                                                     ,address))))
          ,address))))
 
+(defmethod generate-qubit-permutation-code-with-look-up-table ((permutation null))
+  (generate-qubit-permutation-code permutation))
+
+(defmethod generate-qubit-permutation-code-with-look-up-table ((permutation permutation))
+  (let* ((num-bits (slot-value permutation 'size))
+         (num-entries (expt 2 num-bits))
+         (table (make-array num-entries :element-type `(unsigned-byte ,+qubit-index-length+)
+                                        :initial-contents (loop :for i :below num-entries :collect (apply-qubit-permutation permutation i)))))
+    (let ((address (gensym "ADDRESS-")))
+      `(lambda (,address)
+         (declare #.qvm::*optimize-dangerously-fast*
+                  (type qvm:amplitude-address ,address)
+                  (values qvm:amplitude-address))
+         (dpb (aref ,table (ldb (byte ,num-bits 0) ,address))
+              (byte ,num-bits 0)
+              ,address)))))
+
+(defun compile-qubit-permutation (permutation)
+  "Compile PERMUTATION and return a compiled function equivalent to (LAMBDA (ADDRESS) (APPLY-QUBIT-PERMUTATION PERMUTATION ADDRESS))."
+  (declare #.qvm::*optimize-dangerously-fast*
+           (type (or null permutation) permutation))
+  ;; Try the fastest method first. Namely, compilation with a look-up
+  ;; table. If that fails, fall back to single transposition or loop
+  ;; unrolling.
+  (let* ((size (if permutation
+                   (slot-value permutation 'size)
+                   0))
+         (function (if (<= size +qubit-index-length+)
+                       (generate-qubit-permutation-code-with-look-up-table permutation)
+                       (generate-qubit-permutation-code permutation))))
+    (declare (type qubit-index size))
+    (qvm::compile-lambda function)))
+
 (defun-inlinable apply-inverse-qubit-permutation (permutation address)
   (apply-qubit-permutation (inverse-permutation permutation) address))
 
 (defun print-qubit-permutation (permutation &optional number-of-qubits
                                               (stream *standard-output*))
   "Print the address permutation induced by PERMUTATION (possibly using up to NUMBER-OF-QUBITS) in STREAM."
-  (let* ((n (or number-of-qubits
-                (1+ (loop :for (a . b) :in (permutation-transpositions permutation)
-                          :maximizing (max a b)))))
+  (let* ((n (or number-of-qubits (permutation-size permutation)))
          (max-value (expt 2 n))
          (aux-control-string (format nil "~~~DD |~~~D,'0B>"
                                      (ceiling (log max-value 10)) n))
diff --git a/dqvm/tests/permutation-tests.lisp b/dqvm/tests/permutation-tests.lisp
index 55d0c042..51a9404c 100644
--- a/dqvm/tests/permutation-tests.lisp
+++ b/dqvm/tests/permutation-tests.lisp
@@ -86,13 +86,22 @@
              (float (/ (- stop start) internal-time-units-per-second)))
 
            (time-apply-qubit-permutation (permutation number-of-qubits)
-             "Measure the time taken by calls to APPLY-QUBIT-PERMUTATION on addresses from 0 to 2^NUMBER-OF-QUBITS."
+             "Measure the time taken by calls to APPLY-QUBIT-PERMUTATION on addresses from 0 to 2^NUMBER-OF-QUBITS - 1."
              (let ((start (get-internal-real-time)))
                (dotimes (x (expt 2 number-of-qubits))
                  (let ((y (apply-qubit-permutation permutation x)))
                    (values x y)))
                (get-elapsed-time-in-seconds start (get-internal-real-time))))
 
+           (time-compiled-qubit-permutation (permutation number-of-qubits)
+             "Measure the time taken by calls to COMPILE-QUBIT-PERMUTATION and its result on addresses from 0 to 2^NUMBER-OF-QUBITS - 1."
+             (let ((start (get-internal-real-time))
+                   (permute (dqvm2::compile-qubit-permutation permutation)))
+               (dotimes (x (expt 2 number-of-qubits))
+                 (let ((y (funcall permute x)))
+                   (values x y)))
+               (get-elapsed-time-in-seconds start (get-internal-real-time))))
+
            ;; (time-map-reordered-amplitudes (permutation number-of-qubits)
            ;;   "Measure the time taken by calls to MAP-REORDERED-AMPLITUDES on addresses from 0 to 2^NUMBER-OF-QUBITS."
            ;;   (let ((nat-tuple (apply #'qvm::nat-tuple
@@ -109,7 +118,8 @@
            (permutation-0 (make-instance 'dqvm2::permutation-transposition :tau tau))
            (permutation-1 (make-instance 'dqvm2::permutation-general
                                          :number-of-transpositions 2
-                                         :transpositions (list (cons 0 tau) (cons tau 0)))))
+                                         :transpositions (list (cons 0 tau) (cons tau 0))))
+           (permutation-2 (make-permutation '((6 . 0) (3 . 1) (4 . 2)))))
 
       (loop :for x :below (expt 2 (1+ tau)) :do
         (is (= (apply-qubit-permutation permutation-0 x)
@@ -119,6 +129,15 @@
                 (time-apply-qubit-permutation permutation-0 number-of-qubits))
              3))
 
+
+      (is (> (/ (time-apply-qubit-permutation permutation-0 number-of-qubits)
+                (time-compiled-qubit-permutation permutation-0 number-of-qubits))
+             3.5))
+
+      (is (> (/ (time-apply-qubit-permutation permutation-2 number-of-qubits)
+                (time-compiled-qubit-permutation permutation-2 number-of-qubits))
+             4))
+
       ;; (is (> (/ (time-map-reordered-amplitudes permutation-0 number-of-qubits)
       ;;           (time-apply-qubit-permutation permutation-0 number-of-qubits))
       ;;        7))