Quantum++  v1.1
A modern C++11 quantum computing library
gates.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2019 Vlad Gheorghiu (vgheorgh@gmail.com)
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef CLASSES_GATES_H_
33 #define CLASSES_GATES_H_
34 
35 namespace qpp {
40 class Gates final : public internal::Singleton<const Gates> // const Singleton
41 {
42  friend class internal::Singleton<const Gates>;
43 
44  public:
45  // One qubit gates
46  cmat Id2{cmat::Identity(2, 2)};
47  cmat H{cmat::Zero(2, 2)};
48  cmat X{cmat::Zero(2, 2)};
49  cmat Y{cmat::Zero(2, 2)};
50  cmat Z{cmat::Zero(2, 2)};
51  cmat S{cmat::Zero(2, 2)};
52  cmat T{cmat::Zero(2, 2)};
53 
54  // two qubit gates
55  cmat CNOT{cmat::Identity(4, 4)};
56  cmat CZ{cmat::Identity(4, 4)};
57  cmat CNOTba{cmat::Zero(4, 4)};
58  cmat SWAP{cmat::Identity(4, 4)};
59 
60  // three qubit gates
61  cmat TOF{cmat::Identity(8, 8)};
62  cmat FRED{cmat::Identity(8, 8)};
63  private:
67  Gates() {
68  H << 1 / std::sqrt(2.), 1 / std::sqrt(2.), 1 / std::sqrt(2.),
69  -1 / std::sqrt(2.);
70  X << 0, 1, 1, 0;
71  Z << 1, 0, 0, -1;
72  Y << 0, -1_i, 1_i, 0;
73  S << 1, 0, 0, 1_i;
74  T << 1, 0, 0, std::exp(1_i * pi / 4.0);
75  CNOT.block(2, 2, 2, 2) = X;
76  CNOTba(0, 0) = 1;
77  CNOTba(1, 3) = 1;
78  CNOTba(2, 2) = 1;
79  CNOTba(3, 1) = 1;
80  CZ(3, 3) = -1;
81 
82  SWAP.block(1, 1, 2, 2) = X;
83  TOF.block(6, 6, 2, 2) = X;
84  FRED.block(4, 4, 4, 4) = SWAP;
85  }
86 
90  ~Gates() = default;
91 
92  public:
93  // variable gates
94 
95  // one qubit gates
96 
105  cmat Rn(double theta, const std::vector<double>& n) const {
106  // EXCEPTION CHECKS
107 
108  // check 3-dimensional vector
109  if (n.size() != 3)
111  "qpp::Gates::Rn()", "n is not a 3-dimensional vector!");
112  // END EXCEPTION CHECKS
113 
114  cmat result(2, 2);
115  result = std::cos(theta / 2) * Id2 -
116  1_i * std::sin(theta / 2) * (n[0] * X + n[1] * Y + n[2] * Z);
117 
118  return result;
119  }
120 
127  cmat RX(double theta) const {
128  // EXCEPTION CHECKS
129 
130  // END EXCEPTION CHECKS
131 
132  return Rn(theta, {1, 0, 0});
133  }
134 
141  cmat RY(double theta) const {
142  // EXCEPTION CHECKS
143 
144  // END EXCEPTION CHECKS
145 
146  return Rn(theta, {0, 1, 0});
147  }
148 
155  cmat RZ(double theta) const {
156  // EXCEPTION CHECKS
157 
158  // END EXCEPTION CHECKS
159 
160  return Rn(theta, {0, 0, 1});
161  }
162 
163  // one quDit gates
164 
174  cmat Zd(idx D = 2) const {
175  // EXCEPTION CHECKS
176 
177  if (D == 0)
178  throw exception::DimsInvalid("qpp::Gates::Zd()");
179  // END EXCEPTION CHECKS
180 
181  cmat result = cmat::Zero(D, D);
182  for (idx i = 0; i < D; ++i)
183  result(i, i) = std::pow(omega(D), static_cast<double>(i));
184 
185  return result;
186  }
187 
194  cmat SWAPd(idx D = 2) const {
195  // EXCEPTION CHECKS
196 
197  if (D == 0)
198  throw exception::DimsInvalid("qpp::Gates::SWAPd()");
199  // END EXCEPTION CHECKS
200 
201  cmat result = cmat::Zero(D * D, D * D);
202 
203 #ifdef WITH_OPENMP_
204 #pragma omp parallel for collapse(2)
205 #endif // WITH_OPENMP_
206  // column major order for speed
207  for (idx j = 0; j < D; ++j)
208  for (idx i = 0; i < D; ++i)
209  result(D * i + j, i + D * j) = 1;
210 
211  return result;
212  }
213 
224  cmat Fd(idx D = 2) const {
225  // EXCEPTION CHECKS
226 
227  if (D == 0)
228  throw exception::DimsInvalid("qpp::Gates::Fd()");
229  // END EXCEPTION CHECKS
230 
231  cmat result(D, D);
232 
233 #ifdef WITH_OPENMP_
234 #pragma omp parallel for collapse(2)
235 #endif // WITH_OPENMP_
236  // column major order for speed
237  for (idx j = 0; j < D; ++j)
238  for (idx i = 0; i < D; ++i)
239  result(i, j) = 1 / std::sqrt(D) *
240  std::pow(omega(D), static_cast<double>(i * j));
241 
242  return result;
243  }
244 
254  cmat Xd(idx D = 2) const {
255  // EXCEPTION CHECKS
256 
257  if (D == 0)
258  throw exception::DimsInvalid("qpp::Gates::Xd()");
259  // END EXCEPTION CHECKS
260 
261  return Fd(D).inverse() * Zd(D) * Fd(D);
262  }
263 
273  template <typename Derived = Eigen::MatrixXcd>
274  Derived Id(idx D = 2) const {
275  // EXCEPTION CHECKS
276 
277  if (D == 0)
278  throw exception::DimsInvalid("qpp::Gates::Id()");
279  // END EXCEPTION CHECKS
280 
281  return Derived::Identity(D, D);
282  }
283 
299  template <typename Derived>
301  CTRL(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& ctrl,
302  const std::vector<idx>& subsys, idx n, idx d = 2) const {
303  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
304 
305  // EXCEPTION CHECKS
306 
307  // check matrix zero size
309  throw exception::ZeroSize("qpp::Gates::CTRL()");
310 
311  // check square matrix
313  throw exception::MatrixNotSquare("qpp::Gates::CTRL()");
314 
315  // check lists zero size
316  if (ctrl.size() == 0)
317  throw exception::ZeroSize("qpp::Gates::CTRL()");
318  if (subsys.size() == 0)
319  throw exception::ZeroSize("qpp::Gates::CTRL()");
320 
321  // check out of range
322  if (n == 0)
323  throw exception::OutOfRange("qpp::Gates::CTRL()");
324 
325  // check valid local dimension
326  if (d == 0)
327  throw exception::DimsInvalid("qpp::Gates::CTRL()");
328 
329  // ctrl + gate subsystem vector
330  std::vector<idx> ctrlgate = ctrl;
331  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys),
332  std::end(subsys));
333  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
334 
335  std::vector<idx> dims(n, d); // local dimensions vector
336 
337  // check that ctrl + gate subsystem is valid
338  // with respect to local dimensions
339  if (!internal::check_subsys_match_dims(ctrlgate, dims))
340  throw exception::SubsysMismatchDims("qpp::Gates::CTRL()");
341 
342  // check that subsys list match the dimension of the matrix
343  using Index = typename dyn_mat<typename Derived::Scalar>::Index;
344  if (rA.rows() !=
345  static_cast<Index>(std::llround(std::pow(d, subsys.size()))))
346  throw exception::DimsMismatchMatrix("qpp::Gates::CTRL()");
347  // END EXCEPTION CHECKS
348 
349  // Use static allocation for speed!
350  idx Cdims[maxn];
351  idx midx_row[maxn];
352  idx midx_col[maxn];
353 
354  idx CdimsA[maxn];
355  idx midxA_row[maxn];
356  idx midxA_col[maxn];
357 
358  idx Cdims_bar[maxn];
359  idx Csubsys_bar[maxn];
360  idx midx_bar[maxn];
361 
362  idx n_gate = subsys.size();
363  idx n_ctrl = ctrl.size();
364  idx n_subsys_bar = n - ctrlgate.size();
365  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
366  idx DA = static_cast<idx>(rA.rows());
367  idx Dsubsys_bar =
368  static_cast<idx>(std::llround(std::pow(d, n_subsys_bar)));
369 
370  // compute the complementary subsystem of ctrlgate w.r.t. dims
371  std::vector<idx> subsys_bar = complement(ctrlgate, n);
372  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
373  std::begin(Csubsys_bar));
374 
375  for (idx k = 0; k < n; ++k) {
376  midx_row[k] = midx_col[k] = 0;
377  Cdims[k] = d;
378  }
379 
380  for (idx k = 0; k < n_subsys_bar; ++k) {
381  Cdims_bar[k] = d;
382  midx_bar[k] = 0;
383  }
384 
385  for (idx k = 0; k < n_gate; ++k) {
386  midxA_row[k] = midxA_col[k] = 0;
387  CdimsA[k] = d;
388  }
389 
393 
394  // run over the complement indexes
395  for (idx i = 0; i < Dsubsys_bar; ++i) {
396  // get the complement row multi-index
397  internal::n2multiidx(i, n_subsys_bar, Cdims_bar, midx_bar);
398  for (idx k = 0; k < d; ++k) {
399  Ak = powm(rA, k); // compute rA^k
400  // run over the subsys row multi-index
401  for (idx a = 0; a < DA; ++a) {
402  // get the subsys row multi-index
403  internal::n2multiidx(a, n_gate, CdimsA, midxA_row);
404 
405  // construct the result row multi-index
406 
407  // first the ctrl part (equal for both row and column)
408  for (idx c = 0; c < n_ctrl; ++c)
409  midx_row[ctrl[c]] = midx_col[ctrl[c]] = k;
410 
411  // then the complement part (equal for column)
412  for (idx c = 0; c < n_subsys_bar; ++c)
413  midx_row[Csubsys_bar[c]] = midx_col[Csubsys_bar[c]] =
414  midx_bar[c];
415 
416  // then the subsys part
417  for (idx c = 0; c < n_gate; ++c)
418  midx_row[subsys[c]] = midxA_row[c];
419 
420  // run over the subsys column multi-index
421  for (idx b = 0; b < DA; ++b) {
422  // get the subsys column multi-index
423  internal::n2multiidx(b, n_gate, CdimsA, midxA_col);
424 
425  // construct the result column multi-index
426  for (idx c = 0; c < n_gate; ++c)
427  midx_col[subsys[c]] = midxA_col[c];
428 
429  // finally write the values
430  result(internal::multiidx2n(midx_row, n, Cdims),
431  internal::multiidx2n(midx_col, n, Cdims)) =
432  Ak(a, b);
433  }
434  }
435  }
436  }
437 
438  return result;
439  }
440 
456  template <typename Derived>
458  expandout(const Eigen::MatrixBase<Derived>& A, idx pos,
459  const std::vector<idx>& dims) const {
460  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
461 
462  // EXCEPTION CHECKS
463 
464  // check zero-size
466  throw exception::ZeroSize("qpp::Gates::expandout()");
467 
468  // check that dims is a valid dimension vector
469  if (!internal::check_dims(dims))
470  throw exception::DimsInvalid("qpp::Gates::expandout()");
471 
472  // check square matrix
474  throw exception::MatrixNotSquare("qpp::Gates::expandout()");
475 
476  // check that position is valid
477  if (pos > dims.size() - 1)
478  throw exception::OutOfRange("qpp::Gates::expandout()");
479 
480  // check that dims[pos] match the dimension of A
481  if (static_cast<idx>(rA.rows()) != dims[pos])
482  throw exception::DimsMismatchMatrix("qpp::Gates::expandout()");
483  // END EXCEPTION CHECKS
484 
485  idx D = std::accumulate(std::begin(dims), std::end(dims),
486  static_cast<idx>(1), std::multiplies<idx>());
489 
490  idx Cdims[maxn];
491  idx midx_row[maxn];
492  idx midx_col[maxn];
493 
494  for (idx k = 0; k < dims.size(); ++k) {
495  midx_row[k] = midx_col[k] = 0;
496  Cdims[k] = dims[k];
497  }
498 
499  // run over the main diagonal multi-indexes
500  for (idx i = 0; i < D; ++i) {
501  // get row multi_index
502  internal::n2multiidx(i, dims.size(), Cdims, midx_row);
503  // get column multi_index (same as row)
504  internal::n2multiidx(i, dims.size(), Cdims, midx_col);
505  // run over the gate row multi-index
506  for (idx a = 0; a < static_cast<idx>(rA.rows()); ++a) {
507  // construct the total row multi-index
508  midx_row[pos] = a;
509 
510  // run over the gate column multi-index
511  for (idx b = 0; b < static_cast<idx>(rA.cols()); ++b) {
512  // construct the total column multi-index
513  midx_col[pos] = b;
514 
515  // finally write the values
516  result(internal::multiidx2n(midx_row, dims.size(), Cdims),
517  internal::multiidx2n(midx_col, dims.size(), Cdims)) =
518  rA(a, b);
519  }
520  }
521  }
522 
523  return result;
524  }
525 
547  template <typename Derived>
549  expandout(const Eigen::MatrixBase<Derived>& A, idx pos,
550  const std::initializer_list<idx>& dims) const {
551  return this->expandout(A, pos, std::vector<idx>(dims));
552  }
553 
570  template <typename Derived>
572  expandout(const Eigen::MatrixBase<Derived>& A, idx pos, idx n,
573  idx d = 2) const {
574  // EXCEPTION CHECKS
575 
576  // check zero size
578  throw exception::ZeroSize("qpp::Gates::expandout()");
579 
580  // check valid dims
581  if (d == 0)
582  throw exception::DimsInvalid("qpp::Gates::expandout()");
583  // END EXCEPTION CHECKS
584 
585  std::vector<idx> dims(n, d); // local dimensions vector
586 
587  return this->expandout(A, pos, dims);
588  }
589 }; /* class Gates */
590 
591 } /* namespace qpp */
592 
593 #endif /* CLASSES_GATES_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
Dimension(s) mismatch matrix size exception.
Definition: exception.h:300
cmat SWAP
SWAP gate.
Definition: gates.h:58
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:75
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:210
Derived Id(idx D=2) const
Identity gate.
Definition: gates.h:274
~Gates()=default
Default destructor.
cplx omega(idx D)
D-th root of unity.
Definition: constants.h:97
Custom exception.
Definition: exception.h:571
cmat RZ(double theta) const
Qubit rotation of theta about the Z axis.
Definition: gates.h:155
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
dyn_mat< typename Derived::Scalar > expandout(const Eigen::MatrixBase< Derived > &A, idx pos, const std::initializer_list< idx > &dims) const
Expands out.
Definition: gates.h:549
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:43
Singleton policy class, used internally to implement the singleton pattern via CRTP (Curiously recurr...
Definition: singleton.h:80
Subsystems mismatch dimensions exception.
Definition: exception.h:365
const Singleton class that implements most commonly used gates
Definition: gates.h:40
Quantum++ main namespace.
Definition: codes.h:35
cmat Xd(idx D=2) const
Generalized X gate for qudits.
Definition: gates.h:254
Invalid dimension(s) exception.
Definition: exception.h:269
cmat S
S gate.
Definition: gates.h:51
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1733
cmat T
T gate.
Definition: gates.h:52
cmat FRED
Fredkin gate.
Definition: gates.h:62
cmat CNOTba
Controlled-NOT target control gate.
Definition: gates.h:57
cmat Zd(idx D=2) const
Generalized Z gate for qudits.
Definition: gates.h:174
cmat H
Hadamard gate.
Definition: gates.h:47
cmat Rn(double theta, const std::vector< double > &n) const
Qubit rotation of theta about the 3-dimensional real (unit) vector n.
Definition: gates.h:105
cmat Id2
Identity gate.
Definition: gates.h:46
cmat X
Pauli Sigma-X gate.
Definition: gates.h:48
dyn_mat< typename Derived::Scalar > powm(const Eigen::MatrixBase< Derived > &A, idx n)
Fast matrix power based on the SQUARE-AND-MULTIPLY algorithm.
Definition: functions.h:751
cmat Fd(idx D=2) const
Quantum Fourier transform gate for qudits.
Definition: gates.h:224
dyn_mat< typename Derived::Scalar > expandout(const Eigen::MatrixBase< Derived > &A, idx pos, const std::vector< idx > &dims) const
Expands out.
Definition: gates.h:458
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
cmat RX(double theta) const
Qubit rotation of theta about the X axis.
Definition: gates.h:127
constexpr double pi
Definition: constants.h:80
cmat SWAPd(idx D=2) const
SWAP gate for qudits.
Definition: gates.h:194
Parameter out of range exception.
Definition: exception.h:515
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
cmat CNOT
Controlled-NOT control target gate.
Definition: gates.h:55
std::size_t idx
Non-negative integer index.
Definition: types.h:39
cmat Z
Pauli Sigma-Z gate.
Definition: gates.h:50
cmat CZ
Controlled-Phase gate.
Definition: gates.h:56
Matrix is not square exception.
Definition: exception.h:149
cmat RY(double theta) const
Qubit rotation of theta about the Y axis.
Definition: gates.h:141
cmat Y
Pauli Sigma-Y gate.
Definition: gates.h:49
Gates()
Initializes the gates.
Definition: gates.h:67
dyn_mat< typename Derived::Scalar > expandout(const Eigen::MatrixBase< Derived > &A, idx pos, idx n, idx d=2) const
Expands out.
Definition: gates.h:572
cmat TOF
Toffoli gate.
Definition: gates.h:61
Object has zero size exception.
Definition: exception.h:134
dyn_mat< typename Derived::Scalar > CTRL(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &ctrl, const std::vector< idx > &subsys, idx n, idx d=2) const
Generates the multi-partite multiple-controlled-A gate in matrix form.
Definition: gates.h:301