Accelerate APPLY-QUBIT-PERMUTATION.

jmbr · jmbr · commit 6afe4c83448d · 2019-09-11T14:13:27.000-07:00
diff --git 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
@@ -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
@@ -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
