cl-quantum/state.lisp

;;;; Functions for creating and applying quantum gates
(in-package :cl-quantum/state)

;; Fix SBCL's strict handling of `defconstant'
(defmacro define-constant (name value &optional doc)
  "Define NAME to be a constant with value VALUE. If NAME is already defined, do
nothing."
  `(defconstant ,name (if (boundp ',name)
                          (symbol-value ',name)
                          ,value)
     ,@(when doc (list doc))))

;;; Gates and Operators:
(define-constant +unset-projector+
  #2A((1 0)
      (0 0))
  "|0><0| projector matrix.")

(define-constant +set-projector+
  #2A((0 0)
      (0 1))
  "|1><1| projector matrix.")

(define-constant +identity-2x2+
    (make-identity-matrix 2)
  "A 2 by 2 identity matrix.")

(define-constant +pauli-x-gate+
  #2A((0 1)
      (1 0))
  "A 2 by 2 Pauli X gate that can be applied to a single qbit.")

(define-constant +pauli-y-gate+
  #2A((0 #C(0 -1))
      (#C(0 -1) 0))
  "A 2 by 2 Pauli Y gate that can be applied to a single qbit.")

(define-constant +pauli-z-gate+
  #2A((1 0)
      (0 -1))
  "A 2 by 2 Pauli Z gate that can be applied to a single qbit.")

(define-constant +hadamard-gate+
    (let ((oort (/ (sqrt 2))))
      (make-array '(2 2) :initial-contents
                  `((,oort ,oort)
                    (,oort ,(- oort)))))
  "A 2 by 2 Hadamard gate that can be applied to a single qbit.")

(define-constant +phase-gate+
  #2A((1 0)
      (0 #C(0 1)))
  "A 2 by 2 phase gate that can be applied to a single qbit.")

(define-constant +pi/8-gate+
    (make-array '(2 2) :initial-contents
                `((1 0)
                  (0 ,(exp (complex 0 (/ pi 4))))))
  "A 2 by 2 phase gate that can be applied to a single qbit.")

(define-constant +cnot-gate+
  #2A((1 0 0 0)
      (0 1 0 0)
      (0 0 0 1)
      (0 0 1 0))
  "A 4 by 4 controlled NOT (Pauli X) gate that applies to a 2 qbit state. The 0
indexed qbit is the target and the 1 indexed qbit is the control.")

(define-constant +cz-gate+
  #2A((1 0 0 0)
      (0 1 0 0)
      (0 0 1 0)
      (0 0 0 -1))
  "A 4 by 4 controlled Pauli Z gate that applies to a 2 qbit state.")

(define-constant +swap-gate+
  #2A((1 0 0 0)
      (0 0 1 0)
      (0 1 0 0)
      (0 0 0 1))
  "A 4 by 4 swap gate that applies to a 2 qbit state.")

(define-constant +ccnot-gate+
  #2A((1 0 0 0 0 0 0 0)
      (0 1 0 0 0 0 0 0)
      (0 0 1 0 0 0 0 0)
      (0 0 0 1 0 0 0 0)
      (0 0 0 0 1 0 0 0)
      (0 0 0 0 0 1 0 0)
      (0 0 0 0 0 0 0 1)
      (0 0 0 0 0 0 1 0))
  "An 8 by 8 2-controlled NOT (Pauli X) gate that applies to a 3 qbit sate. The
0 indexed qbit is the target and the 1 and 2 indexed qbits are the controls.")

;;; State Functions:
(defun normal-state-p (state &key (places 5))
  "Return non-nil if state is normalized. PLACES is the number of places to
round the norm of STATE before checking."
  (= (round-to-place (dot state state) places) 1))

(defun normalize-state (state)
  "Return a copy of STATE that is normalized."
  (/vs state (norm state)))

(defun nnormalize-state (state)
  "Normalize STATE by mutating it."
  (let ((norm (norm state)))
    (dotimes (i (length state) state)
      (setf (aref state i) (/ (aref state i) norm)))))

(defun state-bits (state)
  "Return the number of bits in STATE."
  (values (ceiling (log (length state) 2))))

(defun make-zero-state (bits)
  "Make a normalized state with a value of zero (probability of 1 in the zero
place)."
  (let ((arr (make-array (ash 1 bits) :initial-element 0)))
    (setf (aref arr 0) 1)
    arr))

(defun make-normal-state (probabilities)
  "Create a new normalized state with the probability of each state
corresponding to an element in PROBABILITIES."
  (map 'vector 'sqrt probabilities))

(defun make-uniform-normal-state (bits)
  "Make a uniform normalized quantum state of BITS qbits."
  (let ((size (ash 1 bits)))
    (make-array (ash 1 bits) :initial-element (/ (sqrt size)))))

(defun bit-unset-index (bit n &key (period (ash 1 bit)))
  "Return the Nth index in a state in which BIT is 0."
  (multiple-value-bind (quo rem)
      (floor n period)
    (+ (* 2 period quo) rem)))

(defun bit-set-index (bit n &key (period (ash 1 bit)))
  "Return the Nth index in a state in which BIT is 1."
  (+ period (bit-unset-index bit n :period period)))

(defun bit-set-p (bit n)
  "Return non-nil if BIT would be set at index N."
  (let ((period (ash 1 bit)))
    (>= (rem n (* 2 period)) period)))

(defun bit-probability (state bit)
  "Return the probability that BIT is set in STATE."
  (setq state (normalize-state state))
  (loop with period = (ash 1 bit)
        for i below (/ (length state) 2)
        for index = (bit-set-index bit i :period period)
        for coef = (aref state index)
        summing (* coef coef)))

(defun nmeasure (state bit &key (places 5))
  "Collapse BIT in STATE by measuring it. This will return t or nil depending on
the state the bit collapsed to. Note that this will also modify STATE. The
probability will be rounded to PLACES before calculations are carried out."
  (loop with prob = (round-to-place (bit-probability state bit) places)
        ;; Theocratically there are multiple numbers that could be equivalent
        ;; to 100% here, but the chance of those failing is so low anyway that
        ;; it doesn't matter.
        with limit = (* most-positive-fixnum prob)
        with rnum = (random most-positive-fixnum)
        with result = (< rnum limit)
        with period = (ash 1 bit)
        with 1-prob = (- 1 prob)
        for i below (/ (length state) 2)
        for unset-index = (bit-unset-index bit i :period period)
        for set-index = (+ period unset-index)
        for unset-coef = (aref state unset-index)
        for set-coef = (aref state set-index)
        if result
          do (setf (aref state unset-index) 0
                   (aref state set-index)
                   (sqrt (/ (* set-coef set-coef) prob)))
        else
          do (setf (aref state unset-index)
                   (sqrt (/ (* unset-coef unset-coef) 1-prob))
                   (aref state set-index) 0)
        finally (return (values result state))))

(defun measure (state bit &key (places 5))
  "Like `nmeasure', but don't modify STATE. Note that the new state is returned
as the second value, with the first value being the result of the
measurement. You will probably want to save the second value as the first value
too as the second value alone means pretty much nothing."
  (nmeasure (copy-seq state) bit :places places))

(defun ncollapse (state &key (places 5))
  "Collapse STATE into a single quantum state. After this call, STATE will be a
state with a zero coefficient in all but one place. This returns two values, the
first is the index into which STATE collapsed and the second is STATE."
  (nnormalize-state state)
  (loop with rval = (random most-positive-fixnum)
        with did-find = nil
        for i below (length state)
        for coef = (aref state i)
        for plimit = 0 then limit
        for prob = (round-to-place (* coef coef) places)
        for limit = (+ plimit (* prob (1- most-positive-fixnum)))
        when (and (not did-find) (< rval limit))
          do (setf (aref state i) 1
                   did-find i)
        else
          do (setf (aref state i) 0)
        finally (return (values did-find state))))

(defun collapse (state &key (places 5))
  "Collapse STATE into a single quantum state. This is like `ncollapse', except
that it does not modify STATE. Thus, you will probably want to keep both
returned values."
  (ncollapse (copy-seq state) :places places))

(defun make-operator (bits operator target)
  "Create an operator matrix that can act on a state with BITS bits and will
apply OPERATOR to TARGET."
  (loop with out = (if (= (1- bits) target)
                       operator
                       +identity-2x2+)
        for i from (- bits 2) downto 0
        do (setq out (tensor-mm out (if (= i target)
                                        operator
                                        +identity-2x2+)))
        finally (return out)))

(defun make-controlled-operator (bits operator target control)
  "Create an operator matrix that can act on a state with BITS bits and will
apply OPERATOR to TARGET if CONTROL is set."
  (labels ((matrix-for (bit target-operator control-operator)
             (cond
               ((= bit target) target-operator)
               ((= bit control) control-operator)
               (t +identity-2x2+)))
           (tensor-chain (target-operator control-operator)
             (loop with out = (matrix-for (1- bits) target-operator
                                          control-operator)
                   for i from (- bits 2) downto 0
                   do (setq out (tensor-mm out (matrix-for i target-operator
                                                           control-operator)))
                   finally (return out))))
    (+mm (tensor-chain +identity-2x2+ +unset-projector+)
         (tensor-chain operator +set-projector+))))

(defun make-n-toffoli-operator (bits target controls)
  "Create an n-Toffoli gate for a state of BITS bits that targets TARGET and is
controlled by CONTROLS."
  (assert (not (member target controls :test '=))
          (target controls)
          "Target ~s cannot be in control set ~s" target controls)
  (let* ((state-size (ash 1 bits))
         (mat (make-identity-matrix state-size))
         (offset (ash 1 target)))
    (loop for i below (/ state-size 2)
          for index = (bit-unset-index target i)
          when (every (lambda (elt)
                        (bit-set-p elt index))
                      controls)
            do (nswap-rows mat index (+ index offset)))
    mat))

(defun make-n-controlled-z-operator (bits target controls)
  "Create an n-controlled z gate for a state of BITS bits that targets TARGET
and is controlled by CONTROLS."
  (assert (not (member target controls :test '=))
          (target controls)
          "Target ~s cannot be in control set ~s" target controls)
  (let* ((state-size (ash 1 bits))
         (mat (make-identity-matrix state-size)))
    (loop for i below (/ state-size 2)
          for index = (bit-set-index target i)
          when (every (lambda (elt)
                        (bit-set-p elt index))
                      controls)
            do (setf (aref mat index index) -1))
    mat))

(defun replace-state (target template)
  "Replace each element of TARGET with the corresponding element in TEMPLATE."
  (assert (= (length target) (length template))
          (target template)
          "Vectors must be of the same length: ~s and ~s" target template)
  (dotimes (i (length target) target)
    (setf (aref target i) (aref template i))))

(defun apply-gate (state operator target)
  "Apply OPERATOR to the bit numbered TARGET in STATE."
  (*mv (make-operator (state-bits state) operator target) state))

(defun napply-gate (state operator target)
  "Apply OPERATOR to the bit numbered TARGET in STATE. This modifies state."
  (replace-state state (apply-gate state operator target)))

(defun apply-controlled-gate (state operator target control)
  "Apply OPERATOR to the bit numbered TARGET in STATE if CONTROL is set."
  (*mv (make-controlled-operator (state-bits state) operator target control)
       state))

(defun napply-controlled-gate (state operator target control)
  "Apply OPERATOR to the bit numbered TARGET in STATE if CONTROL is set. This
modified STATE."
  (replace-state state (apply-controlled-gate state operator target control)))