Quantum++  v1.2
A modern C++11 quantum computing library
functions.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 FUNCTIONS_H_
33 #define FUNCTIONS_H_
34 
35 namespace qpp {
36 // Eigen function wrappers
44 template <typename Derived>
45 dyn_mat<typename Derived::Scalar>
46 transpose(const Eigen::MatrixBase<Derived>& A) {
47  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
48 
49  // EXCEPTION CHECKS
50 
51  // check zero-size
53  throw exception::ZeroSize("qpp::transpose()");
54  // END EXCEPTION CHECKS
55 
56  return rA.transpose();
57 }
58 
66 template <typename Derived>
67 dyn_mat<typename Derived::Scalar>
68 conjugate(const Eigen::MatrixBase<Derived>& A) {
69  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
70 
71  // EXCEPTION CHECKS
72 
73  // check zero-size
75  throw exception::ZeroSize("qpp::conjugate()");
76  // END EXCEPTION CHECKS
77 
78  return rA.conjugate();
79 }
80 
88 template <typename Derived>
89 dyn_mat<typename Derived::Scalar> adjoint(const Eigen::MatrixBase<Derived>& A) {
90  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
91 
92  // EXCEPTION CHECKS
93 
94  // check zero-size
96  throw exception::ZeroSize("qpp::adjoint()");
97  // END EXCEPTION CHECKS
98 
99  return rA.adjoint();
100 }
101 
109 template <typename Derived>
110 dyn_mat<typename Derived::Scalar> inverse(const Eigen::MatrixBase<Derived>& A) {
111  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
112 
113  // EXCEPTION CHECKS
114 
115  // check zero-size
117  throw exception::ZeroSize("qpp::inverse()");
118  // END EXCEPTION CHECKS
119 
120  return rA.inverse();
121 }
122 
129 template <typename Derived>
130 typename Derived::Scalar trace(const Eigen::MatrixBase<Derived>& A) {
131  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
132 
133  // EXCEPTION CHECKS
134 
135  // check zero-size
137  throw exception::ZeroSize("qpp::trace()");
138  // END EXCEPTION CHECKS
139 
140  return rA.trace();
141 }
142 
150 template <typename Derived>
151 typename Derived::Scalar det(const Eigen::MatrixBase<Derived>& A) {
152  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
153 
154  // EXCEPTION CHECKS
155 
156  // check zero-size
158  throw exception::ZeroSize("qpp::det()");
159  // END EXCEPTION CHECKS
160 
161  return rA.determinant();
162 }
163 
173 template <typename Derived>
174 typename Derived::Scalar logdet(const Eigen::MatrixBase<Derived>& A) {
175  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
176 
177  // EXCEPTION CHECKS
178 
179  // check zero-size
181  throw exception::ZeroSize("qpp::logdet()");
182 
183  // check square matrix
185  throw exception::MatrixNotSquare("qpp::logdet()");
186  // END EXCEPTION CHECKS
187 
188  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
190  lu.matrixLU().template triangularView<Eigen::Upper>();
191  typename Derived::Scalar result = std::log(U(0, 0));
192 
193  for (idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
194  result += std::log(U(i, i));
195 
196  return result;
197 }
198 
206 template <typename Derived>
207 typename Derived::Scalar sum(const Eigen::MatrixBase<Derived>& A) {
208  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
209 
210  // EXCEPTION CHECKS
211 
212  // check zero-size
214  throw exception::ZeroSize("qpp::sum()");
215  // END EXCEPTION CHECKS
216 
217  return rA.sum();
218 }
219 
227 template <typename Derived>
228 typename Derived::Scalar prod(const Eigen::MatrixBase<Derived>& A) {
229  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
230 
231  // EXCEPTION CHECKS
232 
233  // check zero-size
235  throw exception::ZeroSize("qpp::prod()");
236  // END EXCEPTION CHECKS
237 
238  return rA.prod();
239 }
240 
247 template <typename Derived>
248 double norm(const Eigen::MatrixBase<Derived>& A) {
249  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
250 
251  // EXCEPTION CHECKS
252 
253  // check zero-size
255  throw exception::ZeroSize("qpp::norm()");
256  // END EXCEPTION CHECKS
257 
258  // convert matrix to complex then return its norm
259  return (rA.template cast<cplx>()).norm();
260 }
261 
268 template <typename Derived>
269 dyn_mat<typename Derived::Scalar>
270 normalize(const Eigen::MatrixBase<Derived>& A) {
271  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
272 
273  // EXCEPTION CHECKS
274 
275  // check zero size
277  throw exception::ZeroSize("qpp::normalize()");
278  // END EXCEPTION CHECKS
279 
281 
283  double normA = norm(rA);
284  try {
285  if (normA == 0)
286  throw std::overflow_error("Division by zero!");
287  } catch (...) {
288  std::cerr << "In qpp::normalize()\n";
289  throw;
290  }
291  result = rA / normA;
292  } else if (internal::check_square_mat(rA)) {
293  typename Derived::Scalar traceA = trace(rA);
294  try {
295  if (std::abs(traceA) == 0)
296  throw std::overflow_error("Division by zero!");
297  } catch (...) {
298  std::cerr << "In qpp::normalize()\n";
299  throw;
300  }
301  result = rA / trace(rA);
302  } else
303  throw exception::MatrixNotSquareNorVector("qpp::normalize()");
304 
305  return result;
306 }
307 
316 template <typename Derived>
317 std::pair<dyn_col_vect<cplx>, cmat> eig(const Eigen::MatrixBase<Derived>& A) {
318  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
319 
320  // EXCEPTION CHECKS
321 
322  // check zero-size
324  throw exception::ZeroSize("qpp::eig()");
325 
326  // check square matrix
328  throw exception::MatrixNotSquare("qpp::eig()");
329  // END EXCEPTION CHECKS
330 
331  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
332 
333  return std::make_pair(es.eigenvalues(), es.eigenvectors());
334 }
335 
343 template <typename Derived>
344 dyn_col_vect<cplx> evals(const Eigen::MatrixBase<Derived>& A) {
345  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
346 
347  // EXCEPTION CHECKS
348 
349  // check zero-size
351  throw exception::ZeroSize("qpp::evals()");
352 
353  // check square matrix
355  throw exception::MatrixNotSquare("qpp::evals()");
356  // END EXCEPTION CHECKS
357 
358  return eig(rA).first;
359 }
360 
368 template <typename Derived>
369 cmat evects(const Eigen::MatrixBase<Derived>& A) {
370  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
371 
372  // EXCEPTION CHECKS
373 
374  // check zero-size
376  throw exception::ZeroSize("qpp::evects()");
377 
378  // check square matrix
380  throw exception::MatrixNotSquare("qpp::evects()");
381  // END EXCEPTION CHECKS
382 
383  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
384 
385  return eig(rA).second;
386 }
387 
396 template <typename Derived>
397 std::pair<dyn_col_vect<double>, cmat>
398 heig(const Eigen::MatrixBase<Derived>& A) {
399  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
400 
401  // EXCEPTION CHECKS
402 
403  // check zero-size
405  throw exception::ZeroSize("qpp::heig()");
406 
407  // check square matrix
409  throw exception::MatrixNotSquare("qpp::heig()");
410  // END EXCEPTION CHECKS
411 
412  Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
413 
414  return std::make_pair(es.eigenvalues(), es.eigenvectors());
415 }
416 
424 template <typename Derived>
425 dyn_col_vect<double> hevals(const Eigen::MatrixBase<Derived>& A) {
426  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
427 
428  // EXCEPTION CHECKS
429 
430  // check zero-size
432  throw exception::ZeroSize("qpp::hevals()");
433 
434  // check square matrix
436  throw exception::MatrixNotSquare("qpp::hevals()");
437  // END EXCEPTION CHECKS
438 
439  return heig(rA).first;
440 }
441 
449 template <typename Derived>
450 cmat hevects(const Eigen::MatrixBase<Derived>& A) {
451  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
452 
453  // EXCEPTION CHECKS
454 
455  // check zero-size
457  throw exception::ZeroSize("qpp::hevects()");
458 
459  // check square matrix
461  throw exception::MatrixNotSquare("qpp::hevects()");
462  // END EXCEPTION CHECKS
463 
464  return heig(rA).second;
465 }
466 
476 template <typename Derived>
477 std::tuple<cmat, dyn_col_vect<double>, cmat>
478 
479 svd(const Eigen::MatrixBase<Derived>& A) {
480  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
481 
482  // EXCEPTION CHECKS
483 
484  // check zero-size
486  throw exception::ZeroSize("qpp::svd()");
487  // END EXCEPTION CHECKS
488 
489  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
490  rA, Eigen::DecompositionOptions::ComputeFullU |
491  Eigen::DecompositionOptions::ComputeFullV);
492 
493  return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
494 }
495 
503 template <typename Derived>
504 dyn_col_vect<double> svals(const Eigen::MatrixBase<Derived>& A) {
505  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
506 
507  // EXCEPTION CHECKS
508 
509  // check zero-size
511  throw exception::ZeroSize("qpp::svals()");
512  // END EXCEPTION CHECKS
513 
514  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
515 
516  return sv.singularValues();
517 }
518 
526 template <typename Derived>
527 cmat svdU(const Eigen::MatrixBase<Derived>& A) {
528  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
529 
530  // EXCEPTION CHECKS
531 
532  // check zero-size
534  throw exception::ZeroSize("qpp::svdU()");
535  // END EXCEPTION CHECKS
536 
537  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
538  rA, Eigen::DecompositionOptions::ComputeFullU);
539 
540  return sv.matrixU();
541 }
542 
550 template <typename Derived>
551 cmat svdV(const Eigen::MatrixBase<Derived>& A) {
552  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
553 
554  // EXCEPTION CHECKS
555 
556  // check zero-size
558  throw exception::ZeroSize("qpp::svdV()");
559  // END EXCEPTION CHECKS
560 
561  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
562  rA, Eigen::DecompositionOptions::ComputeFullV);
563 
564  return sv.matrixV();
565 }
566 
567 // Matrix functional calculus
568 
576 template <typename Derived>
577 cmat funm(const Eigen::MatrixBase<Derived>& A, cplx (*f)(const cplx&)) {
578  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
579 
580  // EXCEPTION CHECKS
581 
582  // check zero-size
584  throw exception::ZeroSize("qpp::funm()");
585 
586  // check square matrix
588  throw exception::MatrixNotSquare("qpp::funm()");
589  // END EXCEPTION CHECKS
590 
591  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
592  cmat evects = es.eigenvectors();
593  cmat evals = es.eigenvalues();
594  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
595  evals(i) = (*f)(evals(i)); // apply f(x) to each eigenvalue
596 
597  cmat evalsdiag = evals.asDiagonal();
598 
599  return evects * evalsdiag * evects.inverse();
600 }
601 
608 template <typename Derived>
609 cmat sqrtm(const Eigen::MatrixBase<Derived>& A) {
610  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
611 
612  // EXCEPTION CHECKS
613 
614  // check zero-size
616  throw exception::ZeroSize("qpp::sqrtm()");
617 
618  // check square matrix
620  throw exception::MatrixNotSquare("qpp::sqrtm()");
621  // END EXCEPTION CHECKS
622 
623  return funm(rA, &std::sqrt);
624 }
625 
632 template <typename Derived>
633 cmat absm(const Eigen::MatrixBase<Derived>& A) {
634  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
635 
636  // EXCEPTION CHECKS
637 
638  // check zero-size
640  throw exception::ZeroSize("qpp::absm()");
641 
642  // check square matrix
644  throw exception::MatrixNotSquare("qpp::absm()");
645  // END EXCEPTION CHECKS
646 
647  return sqrtm(adjoint(rA) * rA);
648 }
649 
656 template <typename Derived>
657 cmat expm(const Eigen::MatrixBase<Derived>& A) {
658  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
659 
660  // EXCEPTION CHECKS
661 
662  // check zero-size
664  throw exception::ZeroSize("qpp::expm()");
665 
666  // check square matrix
668  throw exception::MatrixNotSquare("qpp::expm()");
669  // END EXCEPTION CHECKS
670 
671  return funm(rA, &std::exp);
672 }
673 
680 template <typename Derived>
681 cmat logm(const Eigen::MatrixBase<Derived>& A) {
682  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
683 
684  // EXCEPTION CHECKS
685 
686  // check zero-size
688  throw exception::ZeroSize("qpp::logm()");
689 
690  // check square matrix
692  throw exception::MatrixNotSquare("qpp::logm()");
693  // END EXCEPTION CHECKS
694 
695  return funm(rA, &std::log);
696 }
697 
704 template <typename Derived>
705 cmat sinm(const Eigen::MatrixBase<Derived>& A) {
706  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
707 
708  // EXCEPTION CHECKS
709 
710  // check zero-size
712  throw exception::ZeroSize("qpp::sinm()");
713 
714  // check square matrix
716  throw exception::MatrixNotSquare("qpp::sinm()");
717  // END EXCEPTION CHECKS
718 
719  return funm(rA, &std::sin);
720 }
721 
728 template <typename Derived>
729 cmat cosm(const Eigen::MatrixBase<Derived>& A) {
730  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
731 
732  // EXCEPTION CHECKS
733 
734  // check zero-size
736  throw exception::ZeroSize("qpp::cosm()");
737 
738  // check square matrix
740  throw exception::MatrixNotSquare("qpp::cosm()");
741  // END EXCEPTION CHECKS
742 
743  return funm(rA, &std::cos);
744 }
745 
757 template <typename Derived>
758 cmat spectralpowm(const Eigen::MatrixBase<Derived>& A, const cplx z) {
759  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
760 
761  // EXCEPTION CHECKS
762 
763  // check zero-size
765  throw exception::ZeroSize("qpp::spectralpowm()");
766 
767  // check square matrix
769  throw exception::MatrixNotSquare("qpp::spectralpowm()");
770  // END EXCEPTION CHECKS
771 
772  // Define A^0 = Id, for z IDENTICALLY zero
773  if (real(z) == 0 && imag(z) == 0)
774  return cmat::Identity(rA.rows(), rA.rows());
775 
776  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
777  cmat evects = es.eigenvectors();
778  cmat evals = es.eigenvalues();
779  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
780  evals(i) = std::pow(evals(i), z);
781 
782  cmat evalsdiag = evals.asDiagonal();
783 
784  return evects * evalsdiag * evects.inverse();
785 }
786 
799 template <typename Derived>
800 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
801  idx n) {
802  // EXCEPTION CHECKS
803 
804  // check zero-size
806  throw exception::ZeroSize("qpp::powm()");
807 
808  // check square matrix
810  throw exception::MatrixNotSquare("qpp::powm()");
811  // END EXCEPTION CHECKS
812 
813  // if n = 1, return the matrix unchanged
814  if (n == 1)
815  return A;
816 
819 
820  // if n = 0, return the identity (as just prepared in result)
821  if (n == 0)
822  return result;
823 
824  dyn_mat<typename Derived::Scalar> cA = A.derived(); // copy
825 
826  // fast matrix power
827  for (; n > 0; n /= 2) {
828  if (n % 2)
829  result = (result * cA).eval();
830  cA = (cA * cA).eval();
831  }
832 
833  return result;
834 }
835 
844 template <typename Derived>
845 double schatten(const Eigen::MatrixBase<Derived>& A, double p) {
846  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
847 
848  // EXCEPTION CHECKS
849 
850  // check zero-size
852  throw exception::ZeroSize("qpp::schatten()");
853  if (p < 1)
854  throw exception::OutOfRange("qpp::schatten()");
855  // END EXCEPTION CHECKS
856 
857  if (p == infty) // infinity norm (largest singular value)
858  return svals(rA)(0);
859 
860  const dyn_col_vect<double> sv = svals(rA);
861  double result = 0;
862  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
863  result += std::pow(sv[i], p);
864 
865  return std::pow(result, 1. / p);
866 }
867 
868 // other functions
869 
878 template <typename OutputScalar, typename Derived>
879 dyn_mat<OutputScalar>
880 cwise(const Eigen::MatrixBase<Derived>& A,
881  OutputScalar (*f)(const typename Derived::Scalar&)) {
882  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
883 
884  // EXCEPTION CHECKS
885 
886  // check zero-size
888  throw exception::ZeroSize("qpp::cwise()");
889  // END EXCEPTION CHECKS
890 
891  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
892 
893 #ifdef WITH_OPENMP_
894 #pragma omp parallel for collapse(2)
895 #endif // WITH_OPENMP_
896  // column major order for speed
897  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
898  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
899  result(i, j) = (*f)(rA(i, j));
900 
901  return result;
902 }
903 
904 // Kronecker product of multiple matrices, preserve return type
905 // variadic template
916 template <typename T>
918  return head;
919 }
920 
931 template <typename T, typename... Args>
932 dyn_mat<typename T::Scalar> kron(const T& head, const Args&... tail) {
933  return internal::kron2(head, kron(tail...));
934 }
935 
945 template <typename Derived>
946 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As) {
947  // EXCEPTION CHECKS
948 
949  if (As.size() == 0)
950  throw exception::ZeroSize("qpp::kron()");
951 
952  for (auto&& elem : As)
953  if (!internal::check_nonzero_size(elem))
954  throw exception::ZeroSize("qpp::kron()");
955  // END EXCEPTION CHECKS
956 
957  dyn_mat<typename Derived::Scalar> result = As[0].derived();
958  for (idx i = 1; i < As.size(); ++i) {
959  result = kron(result, As[i]);
960  }
961 
962  return result;
963 }
964 
965 // Kronecker product of a list of matrices, preserve return type
966 // deduce the template parameters from initializer_list
977 template <typename Derived>
978 dyn_mat<typename Derived::Scalar>
979 kron(const std::initializer_list<Derived>& As) {
980  return kron(std::vector<Derived>(As));
981 }
982 
992 template <typename Derived>
993 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
994  idx n) {
995  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
996 
997  // EXCEPTION CHECKS
998 
999  // check zero-size
1001  throw exception::ZeroSize("qpp::kronpow()");
1002 
1003  // check out of range
1004  if (n == 0)
1005  throw exception::OutOfRange("qpp::kronpow()");
1006  // END EXCEPTION CHECKS
1007 
1008  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1009 
1010  return kron(As);
1011 }
1012 
1013 // Direct sum of multiple matrices, preserve return type
1014 // variadic template
1025 template <typename T>
1027  return head;
1028 }
1029 
1040 template <typename T, typename... Args>
1041 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args&... tail) {
1042  return internal::dirsum2(head, dirsum(tail...));
1043 }
1044 
1054 template <typename Derived>
1055 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As) {
1056  // EXCEPTION CHECKS
1057 
1058  if (As.size() == 0)
1059  throw exception::ZeroSize("qpp::dirsum()");
1060 
1061  for (auto&& elem : As)
1062  if (!internal::check_nonzero_size(elem))
1063  throw exception::ZeroSize("qpp::dirsum()");
1064  // END EXCEPTION CHECKS
1065 
1066  idx total_rows = 0, total_cols = 0;
1067  for (idx i = 0; i < As.size(); ++i) {
1068  total_rows += static_cast<idx>(As[i].rows());
1069  total_cols += static_cast<idx>(As[i].cols());
1070  }
1072  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1073 
1074  idx cur_row = 0, cur_col = 0;
1075  for (idx i = 0; i < As.size(); ++i) {
1076  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1077  cur_row += static_cast<idx>(As[i].rows());
1078  cur_col += static_cast<idx>(As[i].cols());
1079  }
1080 
1081  return result;
1082 }
1083 
1084 // Direct sum of a list of matrices, preserve return type
1085 // deduce the template parameters from initializer_list
1096 template <typename Derived>
1097 dyn_mat<typename Derived::Scalar>
1098 dirsum(const std::initializer_list<Derived>& As) {
1099  return dirsum(std::vector<Derived>(As));
1100 }
1101 
1111 template <typename Derived>
1112 dyn_mat<typename Derived::Scalar> dirsumpow(const Eigen::MatrixBase<Derived>& A,
1113  idx n) {
1114  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1115 
1116  // EXCEPTION CHECKS
1117 
1118  // check zero-size
1120  throw exception::ZeroSize("qpp::dirsumpow()");
1121 
1122  // check out of range
1123  if (n == 0)
1124  throw exception::OutOfRange("qpp::dirsumpow()");
1125  // END EXCEPTION CHECKS
1126 
1127  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1128 
1129  return dirsum(As);
1130 }
1131 
1143 template <typename Derived>
1144 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1145  idx rows, idx cols) {
1146  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1147 
1148  idx Arows = static_cast<idx>(rA.rows());
1149  idx Acols = static_cast<idx>(rA.cols());
1150 
1151  // EXCEPTION CHECKS
1152 
1153  // check zero-size
1155  throw exception::ZeroSize("qpp::reshape()");
1156 
1157  if (Arows * Acols != rows * cols)
1158  throw exception::DimsMismatchMatrix("qpp::reshape()");
1159  // END EXCEPTION CHECKS
1160 
1161  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1162  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1163 }
1164 
1177 template <typename Derived1, typename Derived2>
1178 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1179  const Eigen::MatrixBase<Derived2>& B) {
1180  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1181  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1182 
1183  // EXCEPTION CHECKS
1184 
1185  // check types
1186  if (!std::is_same<typename Derived1::Scalar,
1187  typename Derived2::Scalar>::value)
1188  throw exception::TypeMismatch("qpp::comm()");
1189 
1190  // check zero-size
1192  throw exception::ZeroSize("qpp::comm()");
1193 
1194  // check square matrices
1196  throw exception::MatrixNotSquare("qpp::comm()");
1197 
1198  // check equal dimensions
1199  if (rA.rows() != rB.rows())
1200  throw exception::DimsNotEqual("qpp::comm()");
1201  // END EXCEPTION CHECKS
1202 
1203  return rA * rB - rB * rA;
1204 }
1205 
1218 template <typename Derived1, typename Derived2>
1219 dyn_mat<typename Derived1::Scalar>
1220 anticomm(const Eigen::MatrixBase<Derived1>& A,
1221  const Eigen::MatrixBase<Derived2>& B) {
1222  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1223  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1224 
1225  // EXCEPTION CHECKS
1226 
1227  // check types
1228  if (!std::is_same<typename Derived1::Scalar,
1229  typename Derived2::Scalar>::value)
1230  throw exception::TypeMismatch("qpp::anticomm()");
1231 
1232  // check zero-size
1234  throw exception::ZeroSize("qpp::anticomm()");
1235 
1236  // check square matrices
1238  throw exception::MatrixNotSquare("qpp::anticomm()");
1239 
1240  // check equal dimensions
1241  if (rA.rows() != rB.rows())
1242  throw exception::DimsNotEqual("qpp::anticomm()");
1243  // END EXCEPTION CHECKS
1244 
1245  return rA * rB + rB * rA;
1246 }
1247 
1257 template <typename Derived>
1258 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& A) {
1259  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1260 
1261  // EXCEPTION CHECKS
1262 
1263  // check zero-size
1265  throw exception::ZeroSize("qpp::prj()");
1266 
1267  // check column vector
1268  if (!internal::check_cvector(rA))
1269  throw exception::MatrixNotCvector("qpp::prj()");
1270  // END EXCEPTION CHECKS
1271 
1272  double normA = norm(rA);
1273  if (normA > 0)
1274  return rA * adjoint(rA) / (normA * normA);
1275  else
1276  return dyn_mat<typename Derived::Scalar>::Zero(rA.rows(), rA.rows());
1277 }
1278 
1286 template <typename Derived>
1287 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& As) {
1288  // EXCEPTION CHECKS
1289 
1290  // check empty list
1292  throw exception::ZeroSize("qpp::grams()");
1293 
1294  for (auto&& elem : As)
1295  if (!internal::check_nonzero_size(elem))
1296  throw exception::ZeroSize("qpp::grams()");
1297 
1298  // check that As[0] is a column vector
1299  if (!internal::check_cvector(As[0]))
1300  throw exception::MatrixNotCvector("qpp::grams()");
1301 
1302  // now check that all the rest match the size of the first vector
1303  for (auto&& elem : As)
1304  if (elem.rows() != As[0].rows() || elem.cols() != 1)
1305  throw exception::DimsNotEqual("qpp::grams()");
1306  // END EXCEPTION CHECKS
1307 
1309  dyn_mat<typename Derived::Scalar>::Identity(As[0].rows(), As[0].rows());
1310 
1312  dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1313 
1314  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1315  // find the first non-zero vector in the list
1316  idx pos = 0;
1317  for (pos = 0; pos < As.size(); ++pos) {
1318  if (norm(As[pos]) > 0) // add it as the first element
1319  {
1320  outvecs.push_back(As[pos]);
1321  break;
1322  }
1323  }
1324 
1325  // start the process
1326  for (idx i = pos + 1; i < As.size(); ++i) {
1327  cut -= prj(outvecs[i - 1 - pos]);
1328  vi = cut * As[i];
1329  outvecs.push_back(vi);
1330  }
1331 
1332  dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1333 
1334  idx cnt = 0;
1335  for (auto&& elem : outvecs) {
1336  double normA = norm(elem);
1337  if (normA > 0) // we add only the non-zero vectors
1338  {
1339  result.col(cnt) = elem / normA;
1340  ++cnt;
1341  }
1342  }
1343 
1344  return result.block(0, 0, As[0].rows(), cnt);
1345 }
1346 
1347 // deduce the template parameters from initializer_list
1355 template <typename Derived>
1356 dyn_mat<typename Derived::Scalar>
1357 grams(const std::initializer_list<Derived>& As) {
1358  return grams(std::vector<Derived>(As));
1359 }
1360 
1368 template <typename Derived>
1369 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A) {
1370  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1371 
1372  // EXCEPTION CHECKS
1373 
1375  throw exception::ZeroSize("qpp::grams()");
1376  // END EXCEPTION CHECKS
1377 
1378  std::vector<dyn_mat<typename Derived::Scalar>> input;
1379 
1380  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1381  input.push_back(rA.col(i));
1382 
1383  return grams<dyn_mat<typename Derived::Scalar>>(input);
1384 }
1385 
1396 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims) {
1397  // EXCEPTION CHECKS
1398 
1399  if (!internal::check_dims(dims))
1400  throw exception::DimsInvalid("qpp::n2multiidx()");
1401 
1402  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1403  static_cast<idx>(1), std::multiplies<idx>()))
1404  throw exception::OutOfRange("qpp::n2multiidx()");
1405  // END EXCEPTION CHECKS
1406 
1407  // double the size for matrices reshaped as vectors
1408  idx result[2 * maxn];
1409  internal::n2multiidx(n, dims.size(), dims.data(), result);
1410 
1411  return std::vector<idx>(result, result + dims.size());
1412 }
1413 
1424 inline idx multiidx2n(const std::vector<idx>& midx,
1425  const std::vector<idx>& dims) {
1426  // EXCEPTION CHECKS
1427 
1428  if (!internal::check_dims(dims))
1429  throw exception::DimsInvalid("qpp::multiidx2n()");
1430 
1431  for (idx i = 0; i < dims.size(); ++i)
1432  if (midx[i] >= dims[i])
1433  throw exception::OutOfRange("qpp::multiidx2n()");
1434  // END EXCEPTION CHECKS
1435 
1436  return internal::multiidx2n(midx.data(), dims.size(), dims.data());
1437 }
1438 
1453 inline ket mket(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1454  idx n = mask.size();
1455 
1456  idx D = std::accumulate(std::begin(dims), std::end(dims),
1457  static_cast<idx>(1), std::multiplies<idx>());
1458 
1459  // EXCEPTION CHECKS
1460 
1461  // check zero size
1462  if (n == 0)
1463  throw exception::ZeroSize("qpp::mket()");
1464  // check valid dims
1465  if (!internal::check_dims(dims))
1466  throw exception::DimsInvalid("qpp::mket()");
1467  // check mask and dims have the same size
1468  if (mask.size() != dims.size())
1469  throw exception::SubsysMismatchDims("qpp::mket()");
1470  // check mask is a valid vector
1471  for (idx i = 0; i < n; ++i)
1472  if (mask[i] >= dims[i])
1473  throw exception::SubsysMismatchDims("qpp::mket()");
1474  // END EXCEPTION CHECKS
1475 
1476  ket result = ket::Zero(D);
1477  idx pos = multiidx2n(mask, dims);
1478  result(pos) = 1;
1479 
1480  return result;
1481 }
1482 
1496 inline ket mket(const std::vector<idx>& mask, idx d = 2) {
1497  idx n = mask.size();
1498  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1499 
1500  // EXCEPTION CHECKS
1501 
1502  // check zero size
1503  if (n == 0)
1504  throw exception::ZeroSize("qpp::mket()");
1505 
1506  // check valid dims
1507  if (d == 0)
1508  throw exception::DimsInvalid("qpp::mket()");
1509 
1510  // check mask is a valid vector
1511  for (idx i = 0; i < n; ++i)
1512  if (mask[i] >= d)
1513  throw exception::SubsysMismatchDims("qpp::mket()");
1514  // END EXCEPTION CHECKS
1515 
1516  ket result = ket::Zero(D);
1517  std::vector<idx> dims(n, d);
1518  idx pos = multiidx2n(mask, dims);
1519  result(pos) = 1;
1520 
1521  return result;
1522 }
1523 
1539 inline cmat mprj(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1540  idx n = mask.size();
1541 
1542  idx D = std::accumulate(std::begin(dims), std::end(dims),
1543  static_cast<idx>(1), std::multiplies<idx>());
1544 
1545  // EXCEPTION CHECKS
1546 
1547  // check zero size
1548  if (n == 0)
1549  throw exception::ZeroSize("qpp::mprj()");
1550  // check valid dims
1551  if (!internal::check_dims(dims))
1552  throw exception::DimsInvalid("qpp::mprj()");
1553  // check mask and dims have the same size
1554  if (mask.size() != dims.size())
1555  throw exception::SubsysMismatchDims("qpp::mprj()");
1556  // check mask is a valid vector
1557  for (idx i = 0; i < n; ++i)
1558  if (mask[i] >= dims[i])
1559  throw exception::SubsysMismatchDims("qpp::mprj()");
1560  // END EXCEPTION CHECKS
1561 
1562  cmat result = cmat::Zero(D, D);
1563  idx pos = multiidx2n(mask, dims);
1564  result(pos, pos) = 1;
1565 
1566  return result;
1567 }
1568 
1584 inline cmat mprj(const std::vector<idx>& mask, idx d = 2) {
1585  idx n = mask.size();
1586  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1587 
1588  // EXCEPTION CHECKS
1589 
1590  // check zero size
1591  if (n == 0)
1592  throw exception::ZeroSize("qpp::mprj()");
1593 
1594  // check valid dims
1595  if (d == 0)
1596  throw exception::DimsInvalid("qpp::mprj()");
1597 
1598  // check mask is a valid vector
1599  for (idx i = 0; i < n; ++i)
1600  if (mask[i] >= d)
1601  throw exception::SubsysMismatchDims("qpp::mprj()");
1602  // END EXCEPTION CHECKS
1603 
1604  cmat result = cmat::Zero(D, D);
1605  std::vector<idx> dims(n, d);
1606  idx pos = multiidx2n(mask, dims);
1607  result(pos, pos) = 1;
1608 
1609  return result;
1610 }
1611 
1620 template <typename InputIterator>
1621 std::vector<double> abssq(InputIterator first, InputIterator last) {
1622  std::vector<double> weights(std::distance(first, last));
1623  std::transform(first, last, std::begin(weights),
1624  [](cplx z) -> double { return std::norm(z); });
1625 
1626  return weights;
1627 }
1628 
1635 template <typename Container>
1636 std::vector<double>
1637 abssq(const Container& c,
1638  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1639 // we need the std::enable_if to SFINAE out Eigen expressions
1640 // that will otherwise match, instead of matching
1641 // the overload below:
1642 
1643 // template<typename Derived>
1644 // abssq(const Eigen::MatrixBase<Derived>& A)
1645 {
1646  return abssq(std::begin(c), std::end(c));
1647 }
1648 
1655 template <typename Derived>
1656 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A) {
1657  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1658 
1659  // EXCEPTION CHECKS
1660 
1661  // check zero-size
1663  throw exception::ZeroSize("qpp::abssq()");
1664  // END EXCEPTION CHECKS
1665 
1666  return abssq(rA.data(), rA.data() + rA.size());
1667 }
1668 
1677 template <typename InputIterator>
1678 typename std::iterator_traits<InputIterator>::value_type
1679 sum(InputIterator first, InputIterator last) {
1680  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1681 
1682  return std::accumulate(first, last, static_cast<value_type>(0));
1683 }
1684 
1692 template <typename Container>
1693 typename Container::value_type
1694 sum(const Container& c,
1695  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1696  return sum(std::begin(c), std::end(c));
1697 }
1698 
1707 template <typename InputIterator>
1708 typename std::iterator_traits<InputIterator>::value_type
1709 prod(InputIterator first, InputIterator last) {
1710  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1711 
1712  return std::accumulate(first, last, static_cast<value_type>(1),
1713  std::multiplies<value_type>());
1714 }
1715 
1723 template <typename Container>
1724 typename Container::value_type
1725 prod(const Container& c,
1726  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1727  return prod(std::begin(c), std::end(c));
1728 }
1729 
1742 template <typename Derived>
1743 dyn_col_vect<typename Derived::Scalar>
1744 rho2pure(const Eigen::MatrixBase<Derived>& A) {
1745  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1746 
1747  // EXCEPTION CHECKS
1748  // check zero-size
1750  throw exception::ZeroSize("qpp::rho2pure()");
1751  // check square matrix
1752  if (!internal::check_square_mat(rA))
1753  throw exception::MatrixNotSquare("qpp::rho2pure()");
1754  // END EXCEPTION CHECKS
1755 
1756  dyn_col_vect<double> tmp_evals = hevals(rA);
1757  cmat tmp_evects = hevects(rA);
1760  // find the non-zero eigenvector
1761  // there is only one, assuming the state is pure
1762  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k) {
1763  if (std::abs(tmp_evals(k)) > 0) {
1764  result = tmp_evects.col(k);
1765  break;
1766  }
1767  }
1768 
1769  return result;
1770 }
1771 
1780 inline std::vector<idx> complement(std::vector<idx> subsys, idx n) {
1781  // EXCEPTION CHECKS
1782 
1783  if (n < subsys.size())
1784  throw exception::OutOfRange("qpp::complement()");
1785  for (idx i = 0; i < subsys.size(); ++i)
1786  if (subsys[i] >= n)
1787  throw exception::SubsysMismatchDims("qpp::complement()");
1788  // END EXCEPTION CHECKS
1789 
1790  std::vector<idx> all(n);
1791  std::vector<idx> subsys_bar(n - subsys.size());
1792 
1793  std::iota(std::begin(all), std::end(all), 0);
1794  std::sort(std::begin(subsys), std::end(subsys));
1795  std::set_difference(std::begin(all), std::end(all), std::begin(subsys),
1796  std::end(subsys), std::begin(subsys_bar));
1797 
1798  return subsys_bar;
1799 }
1800 
1811 template <typename Derived>
1812 std::vector<double> rho2bloch(const Eigen::MatrixBase<Derived>& A) {
1813  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1814 
1815  // EXCEPTION CHECKS
1816 
1817  // check qubit matrix
1819  throw exception::NotQubitMatrix("qpp::rho2bloch()");
1820  // END EXCEPTION CHECKS
1821 
1822  std::vector<double> result(3);
1823  cmat X(2, 2), Y(2, 2), Z(2, 2);
1824  X << 0, 1, 1, 0;
1825  Y << 0, -1_i, 1_i, 0;
1826  Z << 1, 0, 0, -1;
1827  result[0] = std::real(trace(rA * X));
1828  result[1] = std::real(trace(rA * Y));
1829  result[2] = std::real(trace(rA * Z));
1830 
1831  return result;
1832 }
1833 
1842 inline cmat bloch2rho(const std::vector<double>& r) {
1843  // EXCEPTION CHECKS
1844 
1845  // check 3-dimensional vector
1846  if (r.size() != 3)
1847  throw exception::CustomException("qpp::bloch2rho",
1848  "r is not a 3-dimensional vector!");
1849  // END EXCEPTION CHECKS
1850 
1851  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1852  X << 0, 1, 1, 0;
1853  Y << 0, -1_i, 1_i, 0;
1854  Z << 1, 0, 0, -1;
1855  Id2 << 1, 0, 0, 1;
1856 
1857  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1858 }
1859 
1860 inline namespace literals {
1861 // Idea taken from http://techblog.altplus.co.jp/entry/2017/11/08/130921
1871 template <char... Bits>
1872 ket operator"" _ket() {
1873  constexpr idx n = sizeof...(Bits);
1874  constexpr char bits[n + 1] = {Bits..., '\0'};
1875  qpp::ket q = qpp::ket::Zero(std::pow(2, n));
1876  idx pos = 0;
1877 
1878  // EXCEPTION CHECKS
1879 
1880  // check valid multi-partite qubit state
1881  for (idx i = 0; i < n; ++i) {
1882  if (bits[i] != '0' && bits[i] != '1')
1883  throw exception::OutOfRange(R"xxx(qpp::operator "" _ket())xxx");
1884  }
1885  // END EXCEPTION CHECKS
1886 
1887  pos = std::stoi(bits, nullptr, 2);
1888  q(pos) = 1;
1889 
1890  return q;
1891 }
1892 
1902 template <char... Bits>
1903 bra operator"" _bra() {
1904  constexpr idx n = sizeof...(Bits);
1905  constexpr char bits[n + 1] = {Bits..., '\0'};
1906  qpp::bra q = qpp::ket::Zero(std::pow(2, n));
1907  idx pos = 0;
1908 
1909  // EXCEPTION CHECKS
1910 
1911  // check valid multi-partite qubit state
1912  for (idx i = 0; i < n; ++i) {
1913  if (bits[i] != '0' && bits[i] != '1')
1914  throw exception::OutOfRange(R"xxx(qpp::operator "" _bra())xxx");
1915  }
1916  // END EXCEPTION CHECKS
1917 
1918  pos = std::stoi(bits, nullptr, 2);
1919  q(pos) = 1;
1920 
1921  return q;
1922 }
1923 
1935 template <char... Bits>
1936 cmat operator"" _prj() {
1937  constexpr idx n = sizeof...(Bits);
1938  constexpr char bits[n + 1] = {Bits..., '\0'};
1939 
1940  // EXCEPTION CHECKS
1941 
1942  // check valid multi-partite qubit state
1943  for (idx i = 0; i < n; ++i) {
1944  if (bits[i] != '0' && bits[i] != '1')
1945  throw exception::OutOfRange(R"xxx(qpp::operator "" _prj())xxx");
1946  }
1947  // END EXCEPTION CHECKS
1948 
1949  return kron(operator""_ket<Bits...>(), operator""_bra<Bits...>());
1950 }
1951 } /* namespace literals */
1952 
1953 namespace internal {
1954 // hash combine, code taken from boost::hash_combine(), see
1955 // https://www.boost.org/doc/libs/1_69_0/doc/html/hash/reference.html#boost.hash_combine
1956 template <class T>
1957 void hash_combine(std::size_t& seed, const T& v) {
1958  std::hash<T> hasher;
1959  seed ^= hasher(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
1960 }
1961 } // namespace internal
1962 
1972 template <typename Derived>
1973 std::size_t hash_eigen(const Eigen::MatrixBase<Derived>& A,
1974  std::size_t seed = 0) {
1975  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1976 
1977  // EXCEPTION CHECKS
1978 
1979  // check zero-size
1981  throw exception::ZeroSize("qpp::hash_eigen()");
1982  // END EXCEPTION CHECKS
1983 
1984  auto* p = rA.data();
1985  idx sizeA = static_cast<idx>(rA.size());
1986  for (idx i = 0; i < sizeA; ++i) {
1987  internal::hash_combine(seed, std::real(p[i]));
1988  internal::hash_combine(seed, std::imag(p[i]));
1989  }
1990 
1991  return seed;
1992 }
1993 
1994 namespace internal {
1999 struct HashEigen {
2000  template <typename Derived>
2001  std::size_t operator()(const Eigen::MatrixBase<Derived>& A) const {
2002  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
2003  return hash_eigen(rA);
2004  }
2005 };
2006 
2013 struct EqualEigen {
2014  template <typename Derived>
2015  bool operator()(const Eigen::MatrixBase<Derived>& A,
2016  const Eigen::MatrixBase<Derived>& B) const {
2017  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
2018  const dyn_mat<typename Derived::Scalar>& rB = B.derived();
2019  if (rA.rows() == rB.rows() && rA.cols() == rB.cols())
2020  return rA == rB ? true : false;
2021  else
2022  return false;
2023  }
2024 };
2025 
2026 } /* namespace internal */
2027 } /* namespace qpp */
2028 
2029 #endif /* FUNCTIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:398
Dimensions not equal exception.
Definition: exception.h:284
Dimension(s) mismatch matrix size exception.
Definition: exception.h:300
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:174
void hash_combine(std::size_t &seed, const T &v)
Definition: functions.h:1957
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:67
Matrix is not square nor vector exception.
Definition: exception.h:239
cmat absm(const Eigen::MatrixBase< Derived > &A)
Matrix absolute value.
Definition: functions.h:633
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:277
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
Custom exception.
Definition: exception.h:600
cmat expm(const Eigen::MatrixBase< Derived > &A)
Matrix exponential.
Definition: functions.h:657
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:228
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:59
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:1026
Functor for hashing Eigen expressions.
Definition: functions.h:1999
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
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:43
std::vector< double > rho2bloch(const Eigen::MatrixBase< Derived > &A)
Computes the 3-dimensional real Bloch vector corresponding to the qubit density matrix A...
Definition: functions.h:1812
Subsystems mismatch dimensions exception.
Definition: exception.h:365
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1220
Matrix is not a column vector exception.
Definition: exception.h:164
Quantum++ main namespace.
Definition: circuits.h:35
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:317
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:248
bool check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:111
cmat spectralpowm(const Eigen::MatrixBase< Derived > &A, const cplx z)
Matrix power.
Definition: functions.h:758
std::size_t operator()(const Eigen::MatrixBase< Derived > &A) const
Definition: functions.h:2001
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:425
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:880
Invalid dimension(s) exception.
Definition: exception.h:269
Checks whether T is compatible with an STL-like iterable container.
Definition: traits.h:71
dyn_col_vect< typename Derived::Scalar > rho2pure(const Eigen::MatrixBase< Derived > &A)
Finds the pure state representation of a matrix proportional to a projector onto a pure state...
Definition: functions.h:1744
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:1112
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:369
cmat funm(const Eigen::MatrixBase< Derived > &A, cplx(*f)(const cplx &))
Functional calculus f(A)
Definition: functions.h:577
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:845
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:527
cmat bloch2rho(const std::vector< double > &r)
Computes the density matrix corresponding to the 3-dimensional real Bloch vector r.
Definition: functions.h:1842
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1396
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:89
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1424
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:504
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:917
Functor for comparing Eigen expressions for equality.
Definition: functions.h:2013
std::size_t hash_eigen(const Eigen::MatrixBase< Derived > &A, std::size_t seed=0)
Computes the hash of en Eigen matrix/vector/expression.
Definition: functions.h:1973
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1539
Derived::Scalar det(const Eigen::MatrixBase< Derived > &A)
Determinant.
Definition: functions.h:151
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:46
dyn_mat< typename Derived::Scalar > kronpow(const Eigen::MatrixBase< Derived > &A, idx n)
Kronecker power.
Definition: functions.h:993
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &A)
Projector.
Definition: functions.h:1258
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
Type mismatch exception.
Definition: exception.h:530
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:344
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors of Hermitian matrix.
Definition: functions.h:450
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:68
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:81
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:130
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:551
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:321
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:729
Argument out of range exception.
Definition: exception.h:515
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
cmat sinm(const Eigen::MatrixBase< Derived > &A)
Matrix sin.
Definition: functions.h:705
std::vector< double > abssq(InputIterator first, InputIterator last)
Computes the absolute values squared of an STL-like range of complex numbers.
Definition: functions.h:1621
Matrix is not 2 x 2 exception.
Definition: exception.h:411
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &As)
Gram-Schmidt orthogonalization.
Definition: functions.h:1287
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1178
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1453
dyn_mat< typename Derived::Scalar > inverse(const Eigen::MatrixBase< Derived > &A)
Inverse.
Definition: functions.h:110
cmat logm(const Eigen::MatrixBase< Derived > &A)
Matrix logarithm.
Definition: functions.h:681
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
bool operator()(const Eigen::MatrixBase< Derived > &A, const Eigen::MatrixBase< Derived > &B) const
Definition: functions.h:2015
cmat sqrtm(const Eigen::MatrixBase< Derived > &A)
Matrix square root.
Definition: functions.h:609
std::size_t idx
Non-negative integer index, make sure you use an unsigned type.
Definition: types.h:39
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
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:207
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1144
dyn_mat< typename Derived::Scalar > normalize(const Eigen::MatrixBase< Derived > &A)
Normalizes state vector (column or row vector) or density matrix.
Definition: functions.h:270
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:479
Object has zero size exception.
Definition: exception.h:134
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:240