XACC
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 - 2020 Vlad Gheorghiu.
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)
110  throw exception::CustomException(
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 
365  template <typename Derived>
366  dyn_mat<typename Derived::Scalar>
367  CTRL(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& ctrl,
368  const std::vector<idx>& target, idx n, idx d = 2,
369  std::vector<idx> shift = {}) const {
370  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
371 
372  // EXCEPTION CHECKS
373 
374  // check matrix zero-size
375  if (!internal::check_nonzero_size(rA))
376  throw exception::ZeroSize("qpp::Gates::CTRL()");
377 
378  // check square matrix
379  if (!internal::check_square_mat(rA))
380  throw exception::MatrixNotSquare("qpp::Gates::CTRL()");
381 
382  // check lists zero-size
383  if (ctrl.empty())
384  throw exception::ZeroSize("qpp::Gates::CTRL()");
385  if (target.empty())
386  throw exception::ZeroSize("qpp::Gates::CTRL()");
387 
388  // check out of range
389  if (n == 0)
390  throw exception::OutOfRange("qpp::Gates::CTRL()");
391 
392  // check valid local dimension
393  if (d == 0)
394  throw exception::DimsInvalid("qpp::Gates::CTRL()");
395 
396  // ctrl + gate subsystem vector
397  std::vector<idx> ctrlgate = ctrl;
398  ctrlgate.insert(std::end(ctrlgate), std::begin(target),
399  std::end(target));
400  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
401 
402  std::vector<idx> dims(n, d); // local dimensions vector
403 
404  // check that ctrl + gate subsystem is valid
405  // with respect to local dimensions
406  if (!internal::check_subsys_match_dims(ctrlgate, dims))
407  throw exception::SubsysMismatchDims("qpp::Gates::CTRL()");
408 
409  // check that target list match the dimension of the matrix
410  using Index = typename dyn_mat<typename Derived::Scalar>::Index;
411  if (rA.rows() !=
412  static_cast<Index>(std::llround(std::pow(d, target.size()))))
413  throw exception::DimsMismatchMatrix("qpp::Gates::CTRL()");
414 
415  // check shift
416  if (!shift.empty() && (shift.size() != ctrl.size()))
417  throw exception::SizeMismatch("qpp::Gates::CTRL()");
418  if (!shift.empty())
419  for (auto&& elem : shift)
420  if (elem >= d)
421  throw exception::OutOfRange("qpp::Gates::CTRL()");
422  // END EXCEPTION CHECKS
423 
424  if (shift.empty())
425  shift = std::vector<idx>(ctrl.size(), 0);
426 
427  // Use static allocation for speed!
428  idx Cdims[maxn];
429  idx midx_row[maxn];
430  idx midx_col[maxn];
431 
432  idx CdimsA[maxn];
433  idx midxA_row[maxn];
434  idx midxA_col[maxn];
435 
436  idx Cdims_bar[maxn];
437  idx Csubsys_bar[maxn];
438  idx midx_bar[maxn];
439 
440  idx n_gate = target.size();
441  idx n_ctrl = ctrl.size();
442  idx n_subsys_bar = n - ctrlgate.size();
443  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
444  idx DA = static_cast<idx>(rA.rows());
445  idx Dsubsys_bar =
446  static_cast<idx>(std::llround(std::pow(d, n_subsys_bar)));
447 
448  // compute the complementary subsystem of ctrlgate w.r.t. dims
449  std::vector<idx> subsys_bar = complement(ctrlgate, n);
450  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
451  std::begin(Csubsys_bar));
452 
453  for (idx k = 0; k < n; ++k) {
454  midx_row[k] = midx_col[k] = 0;
455  Cdims[k] = d;
456  }
457 
458  for (idx k = 0; k < n_subsys_bar; ++k) {
459  Cdims_bar[k] = d;
460  midx_bar[k] = 0;
461  }
462 
463  for (idx k = 0; k < n_gate; ++k) {
464  midxA_row[k] = midxA_col[k] = 0;
465  CdimsA[k] = d;
466  }
467 
468  dyn_mat<typename Derived::Scalar> result =
469  dyn_mat<typename Derived::Scalar>::Identity(D, D);
470  dyn_mat<typename Derived::Scalar> Ak;
471 
472  // run over the complement indexes
473  for (idx i = 0; i < Dsubsys_bar; ++i) {
474  // get the complement row multi-index
475  internal::n2multiidx(i, n_subsys_bar, Cdims_bar, midx_bar);
476  for (idx k = 0; k < d; ++k) {
477  Ak = powm(rA, k); // compute rA^k
478  // run over the target row multi-index
479  for (idx a = 0; a < DA; ++a) {
480  // get the target row multi-index
481  internal::n2multiidx(a, n_gate, CdimsA, midxA_row);
482 
483  // construct the result row multi-index
484 
485  // first the ctrl part (equal for both row and column)
486  for (idx c = 0; c < n_ctrl; ++c)
487  midx_row[ctrl[c]] = midx_col[ctrl[c]] =
488  (k + d - shift[c]) % d;
489 
490  // then the complement part (equal for column)
491  for (idx c = 0; c < n_subsys_bar; ++c)
492  midx_row[Csubsys_bar[c]] = midx_col[Csubsys_bar[c]] =
493  midx_bar[c];
494 
495  // then the target part
496  for (idx c = 0; c < n_gate; ++c)
497  midx_row[target[c]] = midxA_row[c];
498 
499  // run over the target column multi-index
500  for (idx b = 0; b < DA; ++b) {
501  // get the target column multi-index
502  internal::n2multiidx(b, n_gate, CdimsA, midxA_col);
503 
504  // construct the result column multi-index
505  for (idx c = 0; c < n_gate; ++c)
506  midx_col[target[c]] = midxA_col[c];
507 
508  // finally write the values
509  result(internal::multiidx2n(midx_row, n, Cdims),
510  internal::multiidx2n(midx_col, n, Cdims)) =
511  Ak(a, b);
512  }
513  }
514  }
515  }
516 
517  return result;
518  }
519 
535  template <typename Derived>
536  dyn_mat<typename Derived::Scalar>
537  expandout(const Eigen::MatrixBase<Derived>& A, idx pos,
538  const std::vector<idx>& dims) const {
539  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
540 
541  // EXCEPTION CHECKS
542 
543  // check zero-size
544  if (!internal::check_nonzero_size(rA))
545  throw exception::ZeroSize("qpp::Gates::expandout()");
546 
547  // check that dims is a valid dimension vector
548  if (!internal::check_dims(dims))
549  throw exception::DimsInvalid("qpp::Gates::expandout()");
550 
551  // check square matrix
552  if (!internal::check_square_mat(rA))
553  throw exception::MatrixNotSquare("qpp::Gates::expandout()");
554 
555  // check that position is valid
556  if (pos + 1 > dims.size())
557  throw exception::OutOfRange("qpp::Gates::expandout()");
558 
559  // check that dims[pos] match the dimension of A
560  if (static_cast<idx>(rA.rows()) != dims[pos])
561  throw exception::DimsMismatchMatrix("qpp::Gates::expandout()");
562  // END EXCEPTION CHECKS
563 
564  idx D = std::accumulate(std::begin(dims), std::end(dims),
565  static_cast<idx>(1), std::multiplies<idx>());
566  dyn_mat<typename Derived::Scalar> result =
567  dyn_mat<typename Derived::Scalar>::Identity(D, D);
568 
569  idx Cdims[maxn];
570  idx midx_row[maxn];
571  idx midx_col[maxn];
572 
573  for (idx k = 0; k < dims.size(); ++k) {
574  midx_row[k] = midx_col[k] = 0;
575  Cdims[k] = dims[k];
576  }
577 
578  // run over the main diagonal multi-indexes
579  for (idx i = 0; i < D; ++i) {
580  // get row multi_index
581  internal::n2multiidx(i, dims.size(), Cdims, midx_row);
582  // get column multi_index (same as row)
583  internal::n2multiidx(i, dims.size(), Cdims, midx_col);
584  // run over the gate row multi-index
585  for (idx a = 0; a < static_cast<idx>(rA.rows()); ++a) {
586  // construct the total row multi-index
587  midx_row[pos] = a;
588 
589  // run over the gate column multi-index
590  for (idx b = 0; b < static_cast<idx>(rA.cols()); ++b) {
591  // construct the total column multi-index
592  midx_col[pos] = b;
593 
594  // finally write the values
595  result(internal::multiidx2n(midx_row, dims.size(), Cdims),
596  internal::multiidx2n(midx_col, dims.size(), Cdims)) =
597  rA(a, b);
598  }
599  }
600  }
601 
602  return result;
603  }
604 
626  template <typename Derived>
627  dyn_mat<typename Derived::Scalar>
628  expandout(const Eigen::MatrixBase<Derived>& A, idx pos,
629  const std::initializer_list<idx>& dims) const {
630  return expandout(A, pos, std::vector<idx>(dims));
631  }
632 
649  template <typename Derived>
650  dyn_mat<typename Derived::Scalar>
651  expandout(const Eigen::MatrixBase<Derived>& A, idx pos, idx n,
652  idx d = 2) const {
653  // EXCEPTION CHECKS
654 
655  // check zero size
656  if (!internal::check_nonzero_size(A))
657  throw exception::ZeroSize("qpp::Gates::expandout()");
658 
659  // check valid dims
660  if (d == 0)
661  throw exception::DimsInvalid("qpp::Gates::expandout()");
662  // END EXCEPTION CHECKS
663 
664  std::vector<idx> dims(n, d); // local dimensions vector
665 
666  return expandout(A, pos, dims);
667  }
668 
669  // getters
670 
680  std::string get_name(const cmat& U) const {
681  // EXCEPTION CHECKS
682 
683  // check zero size
684  if (!internal::check_nonzero_size(U))
685  throw exception::ZeroSize("qpp::Gates::get_name()");
686 
687  // check square matrix
688  if (!internal::check_square_mat(U))
689  return "";
690 
691  // END EXCEPTION CHECKS
692 
693  const idx D = static_cast<idx>(U.rows());
694 
695  switch (D) {
696  // 1 qubit gates
697  case 2:
698  if (U == Id2)
699  return "Id2";
700  else if (U == H)
701  return "H";
702  else if (U == X)
703  return "X";
704  else if (U == Y)
705  return "Y";
706  else if (U == Z)
707  return "Z";
708  else if (U == S)
709  return "S";
710  else if (U == T)
711  return "T";
712  else
713  return "";
714  break;
715  // 2 qubit gates
716  case 4:
717  if (U == CNOT)
718  return "CNOT";
719  else if (U == CZ)
720  return "CZ";
721  else if (U == CNOTba)
722  return "CNOTba";
723  else if (U == SWAP)
724  return "SWAP";
725  else
726  return "";
727  break;
728  // 3 qubit gates
729  case 8:
730  if (U == TOF)
731  return "TOF";
732  else if (U == FRED)
733  return "FRED";
734  else
735  return "";
736  break;
737 
738  default:
739  return "";
740  }
741  }
742  // end getters
743 }; /* class Gates */
744 
745 } /* namespace qpp */
746 
747 #endif /* CLASSES_GATES_H_ */
qpp
Quantum++ main namespace.
Definition: circuits.h:35
Generated by   doxygen 1.8.11