Quantum++  v1.2
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  // check valid dimension
178  if (D == 0)
179  throw exception::DimsInvalid("qpp::Gates::Zd()");
180  // END EXCEPTION CHECKS
181 
182  cmat result = cmat::Zero(D, D);
183  for (idx i = 0; i < D; ++i)
184  result(i, i) = std::pow(omega(D), static_cast<double>(i));
185 
186  return result;
187  }
188 
195  cmat SWAPd(idx D = 2) const {
196  // EXCEPTION CHECKS
197 
198  // check valid dimension
199  if (D == 0)
200  throw exception::DimsInvalid("qpp::Gates::SWAPd()");
201  // END EXCEPTION CHECKS
202 
203  cmat result = cmat::Zero(D * D, D * D);
204 
205 #ifdef WITH_OPENMP_
206 #pragma omp parallel for collapse(2)
207 #endif // WITH_OPENMP_
208  // column major order for speed
209  for (idx j = 0; j < D; ++j)
210  for (idx i = 0; i < D; ++i)
211  result(D * i + j, i + D * j) = 1;
212 
213  return result;
214  }
215 
226  cmat Fd(idx D = 2) const {
227  // EXCEPTION CHECKS
228 
229  // check valid dimension
230  if (D == 0)
231  throw exception::DimsInvalid("qpp::Gates::Fd()");
232  // END EXCEPTION CHECKS
233 
234  cmat result(D, D);
235 
236 #ifdef WITH_OPENMP_
237 #pragma omp parallel for collapse(2)
238 #endif // WITH_OPENMP_
239  // column major order for speed
240  for (idx j = 0; j < D; ++j)
241  for (idx i = 0; i < D; ++i)
242  result(i, j) = 1 / std::sqrt(D) *
243  std::pow(omega(D), static_cast<double>(i * j));
244 
245  return result;
246  }
247 
263  cmat MODMUL(idx a, idx N, idx n) const {
264 
265 // check co-primality (unitarity) only in DEBUG version
266 #ifndef NDEBUG
267  assert(gcd(a, N) == 1);
268 #endif
269  // EXCEPTION CHECKS
270 
271  // check valid arguments
272  if (N < 3 || a >= N) {
273  throw exception::OutOfRange("qpp::Gates::MODMUL()");
274  }
275 
276  // check enough qubits
277  if (n < static_cast<idx>(std::ceil(std::log2(N)))) {
278  throw exception::OutOfRange("qpp::Gates::MODMUL()");
279  }
280  // END EXCEPTION CHECKS
281 
282  // minimum number of qubits required to implement the gate
283  idx D = static_cast<idx>(std::llround(std::pow(2, n)));
284 
285  cmat result = cmat::Zero(D, D);
286 
287 #ifdef WITH_OPENMP_
288 #pragma omp parallel for collapse(2)
289 #endif // WITH_OPENMP_
290  // column major order for speed
291  for (idx j = 0; j < N; ++j)
292  for (idx i = 0; i < N; ++i)
293  if (static_cast<idx>(modmul(j, a, N)) == i)
294  result(i, j) = 1;
295 
296 #ifdef WITH_OPENMP_
297 #pragma omp parallel for
298 #endif // WITH_OPENMP_
299  // complete the matrix
300  for (idx i = N; i < D; ++i)
301  result(i, i) = 1;
302 
303  return result;
304  }
305 
315  cmat Xd(idx D = 2) const {
316  // EXCEPTION CHECKS
317 
318  // check valid dimension
319  if (D == 0)
320  throw exception::DimsInvalid("qpp::Gates::Xd()");
321  // END EXCEPTION CHECKS
322 
323  return Fd(D).inverse() * Zd(D) * Fd(D);
324  }
325 
335  template <typename Derived = Eigen::MatrixXcd>
336  Derived Id(idx D = 2) const {
337  // EXCEPTION CHECKS
338 
339  // check valid dimension
340  if (D == 0)
341  throw exception::DimsInvalid("qpp::Gates::Id()");
342  // END EXCEPTION CHECKS
343 
344  return Derived::Identity(D, D);
345  }
346 
362  template <typename Derived>
364  CTRL(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& ctrl,
365  const std::vector<idx>& target, idx n, idx d = 2) const {
366  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
367 
368  // EXCEPTION CHECKS
369 
370  // check matrix zero-size
372  throw exception::ZeroSize("qpp::Gates::CTRL()");
373 
374  // check square matrix
376  throw exception::MatrixNotSquare("qpp::Gates::CTRL()");
377 
378  // check lists zero-size
379  if (ctrl.size() == 0)
380  throw exception::ZeroSize("qpp::Gates::CTRL()");
381  if (target.size() == 0)
382  throw exception::ZeroSize("qpp::Gates::CTRL()");
383 
384  // check out of range
385  if (n == 0)
386  throw exception::OutOfRange("qpp::Gates::CTRL()");
387 
388  // check valid local dimension
389  if (d == 0)
390  throw exception::DimsInvalid("qpp::Gates::CTRL()");
391 
392  // ctrl + gate subsystem vector
393  std::vector<idx> ctrlgate = ctrl;
394  ctrlgate.insert(std::end(ctrlgate), std::begin(target),
395  std::end(target));
396  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
397 
398  std::vector<idx> dims(n, d); // local dimensions vector
399 
400  // check that ctrl + gate subsystem is valid
401  // with respect to local dimensions
402  if (!internal::check_subsys_match_dims(ctrlgate, dims))
403  throw exception::SubsysMismatchDims("qpp::Gates::CTRL()");
404 
405  // check that target list match the dimension of the matrix
406  using Index = typename dyn_mat<typename Derived::Scalar>::Index;
407  if (rA.rows() !=
408  static_cast<Index>(std::llround(std::pow(d, target.size()))))
409  throw exception::DimsMismatchMatrix("qpp::Gates::CTRL()");
410  // END EXCEPTION CHECKS
411 
412  // Use static allocation for speed!
413  idx Cdims[maxn];
414  idx midx_row[maxn];
415  idx midx_col[maxn];
416 
417  idx CdimsA[maxn];
418  idx midxA_row[maxn];
419  idx midxA_col[maxn];
420 
421  idx Cdims_bar[maxn];
422  idx Csubsys_bar[maxn];
423  idx midx_bar[maxn];
424 
425  idx n_gate = target.size();
426  idx n_ctrl = ctrl.size();
427  idx n_subsys_bar = n - ctrlgate.size();
428  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
429  idx DA = static_cast<idx>(rA.rows());
430  idx Dsubsys_bar =
431  static_cast<idx>(std::llround(std::pow(d, n_subsys_bar)));
432 
433  // compute the complementary subsystem of ctrlgate w.r.t. dims
434  std::vector<idx> subsys_bar = complement(ctrlgate, n);
435  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
436  std::begin(Csubsys_bar));
437 
438  for (idx k = 0; k < n; ++k) {
439  midx_row[k] = midx_col[k] = 0;
440  Cdims[k] = d;
441  }
442 
443  for (idx k = 0; k < n_subsys_bar; ++k) {
444  Cdims_bar[k] = d;
445  midx_bar[k] = 0;
446  }
447 
448  for (idx k = 0; k < n_gate; ++k) {
449  midxA_row[k] = midxA_col[k] = 0;
450  CdimsA[k] = d;
451  }
452 
456 
457  // run over the complement indexes
458  for (idx i = 0; i < Dsubsys_bar; ++i) {
459  // get the complement row multi-index
460  internal::n2multiidx(i, n_subsys_bar, Cdims_bar, midx_bar);
461  for (idx k = 0; k < d; ++k) {
462  Ak = powm(rA, k); // compute rA^k
463  // run over the target row multi-index
464  for (idx a = 0; a < DA; ++a) {
465  // get the target row multi-index
466  internal::n2multiidx(a, n_gate, CdimsA, midxA_row);
467 
468  // construct the result row multi-index
469 
470  // first the ctrl part (equal for both row and column)
471  for (idx c = 0; c < n_ctrl; ++c)
472  midx_row[ctrl[c]] = midx_col[ctrl[c]] = k;
473 
474  // then the complement part (equal for column)
475  for (idx c = 0; c < n_subsys_bar; ++c)
476  midx_row[Csubsys_bar[c]] = midx_col[Csubsys_bar[c]] =
477  midx_bar[c];
478 
479  // then the target part
480  for (idx c = 0; c < n_gate; ++c)
481  midx_row[target[c]] = midxA_row[c];
482 
483  // run over the target column multi-index
484  for (idx b = 0; b < DA; ++b) {
485  // get the target column multi-index
486  internal::n2multiidx(b, n_gate, CdimsA, midxA_col);
487 
488  // construct the result column multi-index
489  for (idx c = 0; c < n_gate; ++c)
490  midx_col[target[c]] = midxA_col[c];
491 
492  // finally write the values
493  result(internal::multiidx2n(midx_row, n, Cdims),
494  internal::multiidx2n(midx_col, n, Cdims)) =
495  Ak(a, b);
496  }
497  }
498  }
499  }
500 
501  return result;
502  }
503 
519  template <typename Derived>
521  expandout(const Eigen::MatrixBase<Derived>& A, idx pos,
522  const std::vector<idx>& dims) const {
523  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
524 
525  // EXCEPTION CHECKS
526 
527  // check zero-size
529  throw exception::ZeroSize("qpp::Gates::expandout()");
530 
531  // check that dims is a valid dimension vector
532  if (!internal::check_dims(dims))
533  throw exception::DimsInvalid("qpp::Gates::expandout()");
534 
535  // check square matrix
537  throw exception::MatrixNotSquare("qpp::Gates::expandout()");
538 
539  // check that position is valid
540  if (pos + 1 > dims.size())
541  throw exception::OutOfRange("qpp::Gates::expandout()");
542 
543  // check that dims[pos] match the dimension of A
544  if (static_cast<idx>(rA.rows()) != dims[pos])
545  throw exception::DimsMismatchMatrix("qpp::Gates::expandout()");
546  // END EXCEPTION CHECKS
547 
548  idx D = std::accumulate(std::begin(dims), std::end(dims),
549  static_cast<idx>(1), std::multiplies<idx>());
552 
553  idx Cdims[maxn];
554  idx midx_row[maxn];
555  idx midx_col[maxn];
556 
557  for (idx k = 0; k < dims.size(); ++k) {
558  midx_row[k] = midx_col[k] = 0;
559  Cdims[k] = dims[k];
560  }
561 
562  // run over the main diagonal multi-indexes
563  for (idx i = 0; i < D; ++i) {
564  // get row multi_index
565  internal::n2multiidx(i, dims.size(), Cdims, midx_row);
566  // get column multi_index (same as row)
567  internal::n2multiidx(i, dims.size(), Cdims, midx_col);
568  // run over the gate row multi-index
569  for (idx a = 0; a < static_cast<idx>(rA.rows()); ++a) {
570  // construct the total row multi-index
571  midx_row[pos] = a;
572 
573  // run over the gate column multi-index
574  for (idx b = 0; b < static_cast<idx>(rA.cols()); ++b) {
575  // construct the total column multi-index
576  midx_col[pos] = b;
577 
578  // finally write the values
579  result(internal::multiidx2n(midx_row, dims.size(), Cdims),
580  internal::multiidx2n(midx_col, dims.size(), Cdims)) =
581  rA(a, b);
582  }
583  }
584  }
585 
586  return result;
587  }
588 
610  template <typename Derived>
612  expandout(const Eigen::MatrixBase<Derived>& A, idx pos,
613  const std::initializer_list<idx>& dims) const {
614  return this->expandout(A, pos, std::vector<idx>(dims));
615  }
616 
633  template <typename Derived>
635  expandout(const Eigen::MatrixBase<Derived>& A, idx pos, idx n,
636  idx d = 2) const {
637  // EXCEPTION CHECKS
638 
639  // check zero size
641  throw exception::ZeroSize("qpp::Gates::expandout()");
642 
643  // check valid dims
644  if (d == 0)
645  throw exception::DimsInvalid("qpp::Gates::expandout()");
646  // END EXCEPTION CHECKS
647 
648  std::vector<idx> dims(n, d); // local dimensions vector
649 
650  return this->expandout(A, pos, dims);
651  }
652 
653  // getters
654 
663  std::string get_name(const cmat& U) const {
664  // EXCEPTION CHECKS
665 
666  // check zero size
668  throw exception::ZeroSize("qpp::Gates::get_name()");
669 
670  // check square matrix
672  return "";
673 
674  // END EXCEPTION CHECKS
675 
676  const idx D = static_cast<idx>(U.rows());
677 
678  switch (D) {
679  // 1 qubit gates
680  case 2:
681  if (U == this->Id2)
682  return "Id2";
683  else if (U == this->H)
684  return "H";
685  else if (U == this->X)
686  return "X";
687  else if (U == this->Y)
688  return "Y";
689  else if (U == this->Z)
690  return "Z";
691  else if (U == this->S)
692  return "S";
693  else if (U == this->T)
694  return "T";
695  else
696  return "";
697  break;
698  // 2 qubit gates
699  case 4:
700  if (U == this->CNOT)
701  return "CNOT";
702  else if (U == this->CZ)
703  return "CZ";
704  else if (U == this->CNOTba)
705  return "CNOTba";
706  else if (U == this->SWAP)
707  return "SWAP";
708  else
709  return "";
710  break;
711  // 3 qubit gates
712  case 8:
713  if (U == this->TOF)
714  return "TOF";
715  else if (U == this->FRED)
716  return "FRED";
717  else
718  return "";
719  break;
720 
721  default:
722  return "";
723  }
724  }
725  // end getters
726 }; /* class Gates */
727 
728 } /* namespace qpp */
729 
730 #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:67
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:219
Derived Id(idx D=2) const
Identity gate.
Definition: gates.h:336
~Gates()=default
Default destructor.
cplx omega(idx D)
D-th root of unity.
Definition: constants.h:89
Custom exception.
Definition: exception.h:600
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:612
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: circuits.h:35
cmat Xd(idx D=2) const
Generalized X gate for qudits.
Definition: gates.h:315
Invalid dimension(s) exception.
Definition: exception.h:269
cmat S
S gate.
Definition: gates.h:51
cmat MODMUL(idx a, idx N, idx n) const
Modular multiplication gate for qubits Implements .
Definition: gates.h:263
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:800
cmat Fd(idx D=2) const
Quantum Fourier transform gate for qudits.
Definition: gates.h:226
dyn_mat< typename Derived::Scalar > expandout(const Eigen::MatrixBase< Derived > &A, idx pos, const std::vector< idx > &dims) const
Expands out.
Definition: gates.h:521
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
std::string get_name(const cmat &U) const
Get the name of the most common qubit gates.
Definition: gates.h:663
constexpr double pi
Definition: constants.h:72
cmat SWAPd(idx D=2) const
SWAP gate for qudits.
Definition: gates.h:195
Argument 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
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:304
cmat CNOT
Controlled-NOT control target gate.
Definition: gates.h:55
std::size_t idx
Non-negative integer index, make sure you use an unsigned type.
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
std::vector< idx > complement(std::vector< idx > subsys, idx n)
Constructs the complement of a subsystem vector.
Definition: functions.h:1780
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:635
cmat TOF
Toffoli gate.
Definition: gates.h:61
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:115
dyn_mat< typename Derived::Scalar > CTRL(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &ctrl, const std::vector< idx > &target, idx n, idx d=2) const
Generates the multi-partite multiple-controlled-A gate in matrix form.
Definition: gates.h:364
Object has zero size exception.
Definition: exception.h:134