QuCoSi/Gate

Go to the documentation of this file.

// QuCoSi - Quantum Computer Simulation
// Copyright © 2009 Frank S. Thomas <frank@blue-dwarf.de>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.

#ifndef QUCOSI_GATE_H
#define QUCOSI_GATE_H

#include <bitset>
#include <cassert>
#include <cmath>
#include <limits>
#include <vector>

#include "Aux"
#include "Vector"

namespace QuCoSi {

/** \class Gate
  *
  * \brief Single- and multi-qubit gate of varying complexity
  *
  * The Gate class is the base for all quantum circuits in QuCoSi. It
  * provides common quantum gates for one to three qubits, more complex gates
  * for an arbitrary number of qubits, and methods to combine or extend them.
  * Since gates are represented by unitary matrices, this class is derived
  * of Eigen's dynamic size matrix class that uses complex numbers. Therefore
  * composition of multiple gates is easily accomplished by simple matrix
  * multiplication of single gates.
  */
class Gate : public MatrixXc
{
  public:
    /** \brief Constructs the 2 × 2 zero matrix
      */
    inline Gate() : MatrixXc(2,2) {}

    /** \brief Constructs the \p rows × \p cols zero matrix
      *
      * \param rows the number of rows of this gate
      * \param cols the number of columns of this gate
      */
    inline Gate(const int rows, const int cols) : MatrixXc(rows,cols) {}

    inline Gate& operator=(const MatrixXc& m)
    {
      MatrixXc::operator=(m);
      return *this;
    }

    /** \brief Computes the tensor product of this gate with \p m
      *
      * The tensor product \f$\mathbf{A} \otimes \mathbf{B}\f$ of two gates
      * \f$\mathbf{A} \in K^{m \times n}\f$ and \f$\mathbf{B} \in K^{p \times
      * q}\f$ is defined as:
      * \f[
      *   \mathbf{A}\, \otimes\, \mathbf{B} =
      *     \left(\begin{array}{ccc}
      *       a_{11} \mathbf{B} & \cdots & a_{1n} \mathbf{B} \\
      *       \vdots & \ddots & \vdots \\
      *       a_{m1} \mathbf{B} & \cdots & a_{mn} \mathbf{B}
      *     \end{array}\right)
      *     =
      *     \left(\begin{array}{cccccccccc}
      *       a_{11} b_{11} & \cdots & a_{11} b_{1q} & \cdots & \cdots &
      *         a_{1n} b_{11} & \cdots & a_{1n} b_{1q} \\
      *       \vdots & \ddots & \vdots & & & \vdots & \ddots & \vdots \\
      *       a_{11} b_{p1} & \cdots & a_{11} b_{pq} & \cdots & \cdots &
      *         a_{1n} b_{p1} & \cdots & a_{1n} b_{pq} \\
      *       \vdots & & \vdots & \ddots & & \vdots & & \vdots \\
      *       \vdots & & \vdots & & \ddots & \vdots & & \vdots \\
      *       a_{m1} b_{11} & \cdots & a_{m1} b_{1q} & \cdots & \cdots &
      *         a_{mn} b_{11} & \cdots & a_{mn} b_{1q} \\
      *       \vdots & \ddots & \vdots & & & \vdots & \ddots & \vdots \\
      *       a_{m1} b_{p1} & \cdots & a_{m1} b_{pq} & \cdots & \cdots &
      *         a_{mn} b_{p1} & \cdots & a_{mn} b_{pq}
      *     \end{array}\right) \in K^{mp \times nq} \ .
      * \f]
      *
      * \param m the right hand side operand of the tensor product
      * \return the tensor product of this gate with Gate \p m
      * \sa http://en.wikipedia.org/wiki/Kronecker_product
      */
    inline Gate tensorDot(const Gate& m) const
    {
      const int r1 = rows();
      const int c1 = cols();
      const int r2 = m.rows();
      const int c2 = m.cols();
      Gate x(r1*r2, c1*c2);

      for (int c = 0; c < c1; ++c) {
        for (int r = 0; r < r1; ++r) {
          x.block(r*r2, c*c2, r2, c2) = (*this)(r,c)*m;
        }
      }
      return x;
    }

    /** \brief Sets the tensor product of this gate and \p m as this gate
      *
      * \param m the right hand side operand of the tensor product
      * \return a reference to \c *this
      * \sa tensorDot()
      */
    inline Gate& tensorDotSet(const Gate& m)
    {
      *this = tensorDot(m);
      return *this;
    }

    /** \brief Computes the <tt>n</tt>th tensor power of this gate
      *
      * The <tt>n</tt>th tensor power of a gate \f$\mathbf{A}\f$ is the
      * <tt>n</tt>-fold tensor product of \f$\mathbf{A}\f$ with itself:
      * \f[
      *   \mathbf{A}^{\otimes n} = \underbrace{\mathbf{A} \otimes \ldots
      *     \otimes \mathbf{A}}_n \ .
      * \f]
      *
      * \param n the exponent of the tensor power
      * \return this gate raised to the <tt>n</tt>th power
      * \sa tensorDot()
      */
    inline Gate tensorPow(const int n) const
    {
      Gate x = *this;
      for (int i = 1; i < n; ++i) {
        x.tensorDotSet(*this);
      }
      return x;
    }

    /** \brief Sets the <tt>n</tt>th tensor power of this gate as this gate
      *
      * \param n the exponent of the tensor power
      * \return a reference to \c *this
      * \sa tensorPow()
      */
    inline Gate& tensorPowSet(const int n)
    {
      *this = tensorPow(n);
      return *this;
    }

    /** \brief Extends this gate to a <tt>n</tt>-qubits gate
      *
      * This method constructs a gate that acts on \p n qubits so that the
      * original gate acts on the qubit(s) at position \p j and all other
      * qubits are left unchanged. This is accomplished by tensor
      * multiplication of an appropriate number of identity gates from the
      * left and the right to the original gate.
      *
      * \param j the position of the qubit(s) the original gate acts on
      * \param n the number of qubits the returned gate acts on
      * \return the for \p n qubits extended gate
      * \sa tensorDot()
      */
    inline Gate applyTo(const int j, const int n) const
    {
      const int k = n-j-log2(rows());
      Gate id, x = *this;

      if (j > 0) {
        id.resize(std::pow(2,j), std::pow(2,j));
        id.setIdentity();
        x = id.tensorDot(x);
      }
      if (j >= 0 && k > 0) {
        id.resize(std::pow(2,k), std::pow(2,k));
        id.setIdentity();
        x.tensorDotSet(id);
      }
      return x;
    }

    /** \brief Sets the to <tt>n</tt>-qubits extended gate as this gate
      *
      * \param j the position of the qubit(s) the original gate acts on
      * \param n the number of qubits the extended gate acts on
      * \return a reference to \c *this
      * \sa applyTo()
      */
    inline Gate& applyToSet(const int j, const int n)
    {
      *this = applyTo(j,n);
      return *this;
    }

    /** \brief \b X gate (\f$\sigma_1\f$ Pauli matrix, NOT gate)
      *
      * \f[\mathbf{X} =
      *   \left(\begin{array}{cc}
      *     0 & 1\\
      *     1 & 0
      *   \end{array}\right)
      * \f]
      *
      * \return a reference to \c *this
      */
    inline Gate& X()
    {
      resize(2,2);
      *this << 0, 1,
               1, 0;
      return *this;
    }

    /** \brief \b Y gate (\f$\sigma_2\f$ Pauli matrix)
      *
      * \f[\mathbf{Y} =
      *   \left(\begin{array}{cc}
      *     0 & -i\\
      *     i & 0
      *   \end{array}\right)
      * \f]
      *
      * \return a reference to \c *this
      */
    inline Gate& Y()
    {
      resize(2,2);
      *this << 0, field(0,-1),
               field(0,1), 0;
      return *this;
    }

    /** \brief \b Z gate (\f$\sigma_3\f$ Pauli matrix)
      *
      * \f[\mathbf{Z} =
      *   \left(\begin{array}{cc}
      *     1 & 0\\
      *     0 & -1
      *   \end{array}\right)
      *   = \mathbf{R}(2)
      * \f]
      *
      * \return a reference to \c *this
      * \sa R()
      */
    inline Gate& Z()
    {
      resize(2,2);
      *this << 1, 0,
               0, -1;
      return *this;
    }

    /** \brief \b H gate (Hadamard gate)
      *
      * \f[\mathbf{H} = \frac{1}{\sqrt{2}}
      *   \left(\begin{array}{cc}
      *     1 & 1\\
      *     1 & -1
      *   \end{array}\right)
      *   = \frac{1}{\sqrt{2}} (\mathbf{X} + \mathbf{Z})
      * \f]
      *
      * \return a reference to \c *this
      */
    inline Gate& H()
    {
      resize(2,2);
      *this << c_sqrt1_2, c_sqrt1_2,
               c_sqrt1_2, -c_sqrt1_2;
      return *this;
    }

    /** \brief \b P gate (phase shift gate)
      *
      * \f[\mathbf{P} =
      *   \left(\begin{array}{cc}
      *     1 & 0\\
      *     0 & i
      *   \end{array}\right)
      *   = \mathbf{R}(4)
      * \f]
      *
      * \return a reference to \c *this
      * \sa R()
      */
    inline Gate& P()
    {
      resize(2,2);
      *this << 1, 0,
               0, field(0,1);
      return *this;
    }

    /** \brief \b T gate (\f$\pi/4\f$ phase shift gate)
      *
      * \f[\mathbf{T} =
      *   \left(\begin{array}{cc}
      *     1 & 0\\
      *     0 & e^{\pi i/4}
      *   \end{array}\right)
      *   = \mathbf{R}(8)
      * \f]
      *
      * \return a reference to \c *this
      * \sa R()
      */
    inline Gate& T()
    {
      resize(2,2);
      *this << 1, 0,
               0, field(c_sqrt1_2,c_sqrt1_2);
      return *this;
    }

    /** \brief <b>R</b>(\p k) gate (general phase shift gate)
      *
      * \f[\mathbf{R}(k) =
      *   \left(\begin{array}{cc}
      *     1 & 0\\
      *     0 & e^{2 \pi i/k}
      *   \end{array}\right)
      * \f]
      *
      * \param k the phase shift of this gate
      * \return a reference to \c *this
      */
    inline Gate& R(const fptype k)
    {
      resize(2,2);
      *this << 1, 0,
               0, std::exp(2*c_pi/k*field(0,1));
      return *this;
    }

    /** \brief \b I gate (identity gate)
      *
      * \f[\mathbf{I} =
      *   \left(\begin{array}{cc}
      *     1 & 0\\
      *     0 & 1
      *   \end{array}\right)
      *   = \mathbf{R}(1)
      * \f]
      *
      * \return a reference to \c *this
      */
    inline Gate& I()
    {
      resize(2,2);
      setIdentity();
      return *this;
    }

    /** \brief <b>R</b><sub>x</sub>(\p theta) gate
      *
      * Rotates a qubit about the x-axis on the Bloch sphere by the angle
      * \p theta.
      *
      * \f[\mathbf{R_x}(\theta) =
      *   \exp \left( -i \frac{\theta}{2} \mathbf{X} \right) =
      *   \cos \frac{\theta}{2} \mathbf{I} - i \sin \frac{\theta}{2}
      *     \mathbf{X} =
      *   \left(\begin{array}{cc}
      *     \cos \frac{\theta}{2} & -i \sin \frac{\theta}{2} \\
      *     -i \sin \frac{\theta}{2} & \cos \frac{\theta}{2}
      *   \end{array}\right)
      * \f]
      *
      * \param theta the angle of rotation about the x-axis
      * \return a reference to \c *this
      */
    inline Gate& Rx(const fptype theta)
    {
      resize(2,2);
      *this << std::cos(theta/2.), std::sin(theta/2.)*field(0,-1),
               std::sin(theta/2.)*field(0,-1), std::cos(theta/2.);
      return *this;
    }

    /** \brief <b>R</b><sub>y</sub>(\p theta) gate
      *
      * Rotates a qubit about the y-axis on the Bloch sphere by the angle
      * \p theta.
      *
      * \f[\mathbf{R_y}(\theta) =
      *   \exp \left( -i \frac{\theta}{2} \mathbf{Y} \right) =
      *   \cos \frac{\theta}{2} \mathbf{I} - i \sin \frac{\theta}{2}
      *     \mathbf{Y} =
      *   \left(\begin{array}{cc}
      *     \cos \frac{\theta}{2} & -\sin \frac{\theta}{2} \\
      *     \sin \frac{\theta}{2} & \cos \frac{\theta}{2}
      *   \end{array}\right)
      * \f]
      *
      * \param theta the angle of rotation about the y-axis
      * \return a reference to \c *this
      */
    inline Gate& Ry(const fptype theta)
    {
      resize(2,2);
      *this << std::cos(theta/2.), -std::sin(theta/2.),
               std::sin(theta/2.), std::cos(theta/2.);
      return *this;
    }

    /** \brief <b>R</b><sub>z</sub>(\p theta) gate
      *
      * Rotates a qubit about the z-axis on the Bloch sphere by the angle
      * \p theta.
      *
      * \f[\mathbf{R_z}(\theta) =
      *   \exp \left( -i \frac{\theta}{2} \mathbf{Z} \right) =
      *   \cos \frac{\theta}{2} \mathbf{I} - i \sin \frac{\theta}{2}
      *     \mathbf{Z} =
      *   \left(\begin{array}{cc}
      *     \exp \left(-i \frac{\theta}{2} \right) & 0 \\
      *     0 & \exp \left(i \frac{\theta}{2} \right)
      *   \end{array}\right)
      * \f]
      *
      * \param theta the angle of rotation about the z-axis
      * \return a reference to \c *this
      */
    inline Gate& Rz(const fptype theta)
    {
      resize(2,2);
      *this << std::exp(theta/2.*field(0,-1)), 0,
               0, std::exp(theta/2.*field(0,1));
      return *this;
    }

    /** \brief \b CNOT gate (controlled NOT gate)
      *
      * \f[\mathbf{CNOT} =
      *   \left(\begin{array}{cccc}
      *     1 & 0 & 0 & 0\\
      *     0 & 1 & 0 & 0\\
      *     0 & 0 & 0 & 1\\
      *     0 & 0 & 1 & 0
      *   \end{array}\right)
      *   = \mathbf{C}_{102}(\mathbf{X})
      * \f]
      *
      * \return a reference to \c *this
      * \sa C(), X()
      */
    inline Gate& CNOT()
    {
      resize(4,4);
      *this << 1, 0, 0, 0,
               0, 1, 0, 0,
               0, 0, 0, 1,
               0, 0, 1, 0;
      return *this;
    }

    /** \brief \b CCNOT gate (Toffoli gate, controlled \b CNOT gate)
      *
      * \f[\mathbf{CCNOT} =
      *   \left(\begin{array}{cccccccc}
      *     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
      *   \end{array}\right)
      *   = \mathbf{C}_{103}(\mathbf{CNOT})
      * \f]
      *
      * \return a reference to \c *this
      * \sa C(), CNOT()
      */
    inline Gate& CCNOT()
    {
      resize(8,8);
      setIdentity();
      (*this)(6,6) = 0;
      (*this)(7,7) = 0;
      (*this)(7,6) = 1;
      (*this)(6,7) = 1;
      return *this;
    }

    /** \brief \b CSWAP gate (Fredkin gate, controlled \b SWAP gate)
      *
      * \f[\mathbf{CSWAP} =
      *   \left(\begin{array}{cccccccc}
      *     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 & 0\\
      *     0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
      *   \end{array}\right)
      *   = \mathbf{C}_{103}(\mathbf{SWAP})
      * \f]
      *
      * \return a reference to \c *this
      * \sa C(), SWAP()
      */
    inline Gate& CSWAP()
    {
      resize(8,8);
      setIdentity();
      (*this)(5,5) = 0;
      (*this)(6,6) = 0;
      (*this)(6,5) = 1;
      (*this)(5,6) = 1;
      return *this;
    }

    /** \brief <b>C</b><sub>\p tcn</sub>(\p U) gate (controlled \p U gate)
      *
      * This method constructs a gate that acts on \p n qubits so that the
      * gate \p U is applied to the qubit(s) at position \p t (the target
      * qubit(s)) if the qubit at position \p c (the control qubit) is 1.
      *
      * \param t the position of the target qubit(s)
      * \param c the position of the control qubit
      * \param n the number of qubits this gate acts on
      * \param U the gate that acts on the target qubit and that is
      *          controlled by the control qubit
      * \return a reference to \c *this
      */
    inline Gate& C(const int t, const int c, const int n, const Gate& U)
    {
      assert(t < n || c < n || t == c);

      // Construct the controlled U gate with the first qubit as control and
      // the second qubit as target.
      const int d = U.rows();
      Gate cu(2*d,2*d);
      cu.setIdentity();
      cu.block(d,d,d,d) = U;

      // Extend the controlled U gate for n qubits if it does not already
      // operate on n qubits.
      if (std::pow(2,n) > 2*d) {
        cu = cu.applyTo(0,n);
      }

      // Construct a cyclic permutation that displaces t to position 1.
      std::vector<int> sigma(n);
      for (int i = 0; i < n; ++i) {
        sigma[(n-t+i+1)%n] = i;
      }
      // Permute c to position 0.
      const int newc = (n-t+c+1)%n, tmp = sigma[0];
      sigma[0] = sigma[newc];
      sigma[newc] = tmp;

      // Construct an S gate from the permutation sigma.
      Gate s;
      s.S(sigma);

      // Now combine the S and the controlled U gate.
      *this = s.transpose()*cu*s;
      return *this;
    }

    /** \brief \b SWAP gate
      *
      * \f[\mathbf{SWAP} =
      *   \left(\begin{array}{cccc}
      *     1 & 0 & 0 & 0\\
      *     0 & 0 & 1 & 0\\
      *     0 & 1 & 0 & 0\\
      *     0 & 0 & 0 & 1
      *   \end{array}\right)
      *   = \mathbf{S}_{102} = \mathbf{S}_{012}
      * \f]
      *
      * \return a reference to \c *this
      * \sa S()
      */
    inline Gate& SWAP()
    {
      resize(4,4);
      *this << 1, 0, 0, 0,
               0, 0, 1, 0,
               0, 1, 0, 0,
               0, 0, 0, 1;
      return *this;
    }

    /** \brief <b>S</b><sub>\p pqn</sub> gate
      *
      * This method constructs a \f$2^n \times 2^n\f$ tensor permutation
      * matrix that permutes the <tt>p</tt>th and <tt>q</tt>th qubits in a
      * tensor product of \p n qubits.
      *
      * \param p the new position of the <tt>q</tt>th qubit
      * \param q the new position of the <tt>p</tt>th qubit
      * \param n the number of qubits this gate acts on
      * \return a reference to \c *this
      */
    inline Gate& S(const int p, const int q, const int n)
    {
      std::vector<int> sigma(n);
      for (int i = 0; i < n; ++i) {
        sigma[i] = i;
      }
      sigma[p] = q;
      sigma[q] = p;
      S(sigma);
      return *this;
    }

    /** \brief <b>S</b>(\p sigma) gate
      *
      * This method constructs a tensor permutation matrix that permutes
      * qubits according to the permutation \p sigma. For example, this
      * matrix acts on the tensor product of the qubits \f$ q_0,\ q_1,\
      * \ldots,\ q_k\f$ as follows:
      * \f[
      *   \mathbf{S}(\sigma) \left(q_0 \otimes q_1 \otimes \ldots \otimes
      *     q_k\right) = q_{\sigma(0)} \otimes q_{\sigma(1)} \otimes \ldots
      *     \otimes q_{\sigma(k)} \ .
      * \f]
      *
      * \note The implementation of this method is based on proposition 6.2
      * in the paper arXiv:math/0508053v2 by Rakotonirina Christian. It takes
      * advantage of the fact that the dimension of single qubits is 2 so
      * that the multiple row and column indices \f$i_1 \ldots i_k\f$ and
      * \f$j_1 \ldots j_k\f$ can be obtained from the row and column indices
      * of the permutation matrix with <tt>std::bitset</tt>s.
      *
      * \param sigma the permutation that will be applied to qubits
      * \return a reference to \c *this
      * \sa http://arxiv.org/abs/math/0508053
      */
    inline Gate& S(const std::vector<int>& sigma)
    {
      const int n = sigma.size();
      const int dim = std::pow(2,n);
      resize(dim,dim);
      setZero();

      // The first and last component of any tensor stays unchanged for all
      // tensor permutations, therefore we ignore the outermost rows and
      // columns.
      (*this)(0,0) = 1;
      (*this)(dim-1,dim-1) = 1;

      for (int c = 1; c < dim-1; ++c) {
        for (int r = 1; r < dim-1; ++r) {
          // Create bitsets of the row and column index. Since the dimension
          // of single qubits is 2, these bitsets can be interpreted as the
          // multiple row and colum indices.
          std::bitset< std::numeric_limits<unsigned long>::digits >
            br(r), bc(c), bcp(0);

          // Permute the bitset bc using the permutation sigma.
          for (int i = 0; i < n; ++i) {
            bcp[n-1-i] = bc[n-1-sigma[i]];
          }

          // If the bitset of the row index is equal to the permuted bitset of
          // the column index, set a 1 at this position. This equals the
          // multiple Kronecker deltas.
          if (br == bcp) {
            (*this)(r,c) = 1;
            break;
          }
        }
      }
      return *this;
    }

    /** \brief <b>U</b><sub>f</sub> gate for one output qubit
      *
      * This method constructs a gate that can be associated with the function
      * \f$f:\,\{0,1,2,\ldots,2^n-1\} \rightarrow \{0,1\}\f$. It is defined so
      * that it acts on a tensor product of n + 1 qubits as
      * \f[
      *   \mathbf{U}_f |x\rangle_n |y\rangle_1 =
      *     |x\rangle_n |y \oplus f(x)\rangle_1 \ ,
      * \f]
      * where \f$\oplus\f$ denotes binary addition.
      *
      * \param f the function associated with this gate
      * \return a reference to \c *this
      */
    inline Gate& U(const std::vector<int>& f)
    {
      const int s = f.size();
      resize(2*s,2*s);
      setIdentity();

      for (int i = 0; i < s; ++i) {
        if (f.at(i) == 1) {
          (*this)(2*i,2*i) = 0;
          (*this)(2*i+1,2*i+1) = 0;
          (*this)(2*i+1,2*i) = 1;
          (*this)(2*i,2*i+1) = 1;
        }
      }
      return *this;
    }

    /** \brief <b>U</b><sub>f</sub> gate for multiple output qubits
      *
      * This method constructs a gate that can be associated with the function
      * \f$f:\,\{0,1,2,\ldots,2^n-1\} \rightarrow \{0,1,2,\ldots,2^m-1\}\f$.
      * It is defined so that it acts on a tensor product of n + \p m qubits
      * as
      * \f[
      *   \mathbf{U}_f |x\rangle_n |y\rangle_m =
      *     |x\rangle_n |y \oplus f(x)\rangle_m \ ,
      * \f]
      * where \f$\oplus\f$ denotes bitwise binary addition.
      *
      * \param f the function associated with this gate
      * \param m the number of output qubits
      * \return a reference to \c *this
      */
    inline Gate& U(const std::vector<int>& f, const int m)
    {
      const int sx = f.size();
      const int sy = std::pow(2,m);
      resize(sx*sy,sx*sy);
      setZero();

      for (int i = 0, j = 0; i < sx; ++i) {
        for (int k = 0; k < sy; ++j, ++k) {
          (*this)(j-k+(k^f.at(i)),j) = 1;
        }
      }
      return *this;
    }

    /** \brief <b>F</b><sub>\p n</sub> gate
      *        (quantum Fourier transform (QFT) gate)
      *
      * This method constructs the quantum Fourier transform gate that acts on
      * \p n qubits. It is defined as:
      * \f[\mathbf{F}_n = 2^{-n/2}
      *   \left(\begin{array}{cccccc}
      *     1 & 1 & 1 & 1 & \cdots & 1 \\
      *     1 & \omega^1 & \omega^2 & \omega^3 & \cdots & \omega^{2^n-1}\\
      *     1 & \omega^2 & \omega^4 & \omega^6 & \cdots & \omega^{2(2^n-1)}\\
      *     1 & \omega^3 & \omega^6 & \omega^9 & \cdots & \omega^{3(2^n-1)}\\
      *     \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
      *     1 & \omega^{2^n-1} & \omega^{2(2^n-1)} & \omega^{3(2^n-1)} &
      *       \cdots & \omega^{(2^n-1)^2}
      *   \end{array}\right),
      *   \quad \omega = \exp\left(2 \pi i/2^n\right) \ .
      * \f]
      *
      * \param n the number of qubits this gate acts on
      * \return a reference to \c *this
      */
    inline Gate& F(const int n)
    {
      const int s = std::pow(2,n);
      resize(s,s);
      for (int x = 0; x < s; ++x) {
        for (int y = x; y < s; ++y) {
          (*this)(x,y) = std::exp(2*c_pi*x*y/s*field(0,1));
        }
      }
      *this *= std::sqrt(1./s);

      MatrixXc t = transpose();
      t.diagonal().setZero();
      *this += t;

      return *this;
    }
};

} // namespace QuCoSi

#endif // QUCOSI_GATE_H

// vim: filetype=cpp shiftwidth=2 textwidth=78