XACC
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 - 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 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
52  if (!internal::check_nonzero_size(rA))
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
74  if (!internal::check_nonzero_size(rA))
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
95  if (!internal::check_nonzero_size(rA))
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
116  if (!internal::check_nonzero_size(rA))
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
136  if (!internal::check_nonzero_size(rA))
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
157  if (!internal::check_nonzero_size(rA))
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
180  if (!internal::check_nonzero_size(rA))
181  throw exception::ZeroSize("qpp::logdet()");
182 
183  // check square matrix
184  if (!internal::check_square_mat(rA))
185  throw exception::MatrixNotSquare("qpp::logdet()");
186  // END EXCEPTION CHECKS
187 
188  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
189  dyn_mat<typename Derived::Scalar> U =
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
213  if (!internal::check_nonzero_size(rA))
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
234  if (!internal::check_nonzero_size(rA))
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
254  if (!internal::check_nonzero_size(rA))
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
276  if (!internal::check_nonzero_size(rA))
277  throw exception::ZeroSize("qpp::normalize()");
278  // END EXCEPTION CHECKS
279 
280  dyn_mat<typename Derived::Scalar> result;
281 
282  if (internal::check_cvector(rA) || internal::check_rvector(rA)) {
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
323  if (!internal::check_nonzero_size(rA))
324  throw exception::ZeroSize("qpp::eig()");
325 
326  // check square matrix
327  if (!internal::check_square_mat(rA))
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
350  if (!internal::check_nonzero_size(rA))
351  throw exception::ZeroSize("qpp::evals()");
352 
353  // check square matrix
354  if (!internal::check_square_mat(rA))
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
375  if (!internal::check_nonzero_size(rA))
376  throw exception::ZeroSize("qpp::evects()");
377 
378  // check square matrix
379  if (!internal::check_square_mat(rA))
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
404  if (!internal::check_nonzero_size(rA))
405  throw exception::ZeroSize("qpp::heig()");
406 
407  // check square matrix
408  if (!internal::check_square_mat(rA))
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
431  if (!internal::check_nonzero_size(rA))
432  throw exception::ZeroSize("qpp::hevals()");
433 
434  // check square matrix
435  if (!internal::check_square_mat(rA))
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
456  if (!internal::check_nonzero_size(rA))
457  throw exception::ZeroSize("qpp::hevects()");
458 
459  // check square matrix
460  if (!internal::check_square_mat(rA))
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
485  if (!internal::check_nonzero_size(rA))
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
510  if (!internal::check_nonzero_size(rA))
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
533  if (!internal::check_nonzero_size(rA))
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
557  if (!internal::check_nonzero_size(rA))
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
583  if (!internal::check_nonzero_size(rA))
584  throw exception::ZeroSize("qpp::funm()");
585 
586  // check square matrix
587  if (!internal::check_square_mat(rA))
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
615  if (!internal::check_nonzero_size(rA))
616  throw exception::ZeroSize("qpp::sqrtm()");
617 
618  // check square matrix
619  if (!internal::check_square_mat(rA))
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
639  if (!internal::check_nonzero_size(rA))
640  throw exception::ZeroSize("qpp::absm()");
641 
642  // check square matrix
643  if (!internal::check_square_mat(rA))
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
663  if (!internal::check_nonzero_size(rA))
664  throw exception::ZeroSize("qpp::expm()");
665 
666  // check square matrix
667  if (!internal::check_square_mat(rA))
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
687  if (!internal::check_nonzero_size(rA))
688  throw exception::ZeroSize("qpp::logm()");
689 
690  // check square matrix
691  if (!internal::check_square_mat(rA))
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
711  if (!internal::check_nonzero_size(rA))
712  throw exception::ZeroSize("qpp::sinm()");
713 
714  // check square matrix
715  if (!internal::check_square_mat(rA))
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
735  if (!internal::check_nonzero_size(rA))
736  throw exception::ZeroSize("qpp::cosm()");
737 
738  // check square matrix
739  if (!internal::check_square_mat(rA))
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
764  if (!internal::check_nonzero_size(rA))
765  throw exception::ZeroSize("qpp::spectralpowm()");
766 
767  // check square matrix
768  if (!internal::check_square_mat(rA))
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
805  if (!internal::check_nonzero_size(A))
806  throw exception::ZeroSize("qpp::powm()");
807 
808  // check square matrix
809  if (!internal::check_square_mat(A))
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 
817  dyn_mat<typename Derived::Scalar> result =
818  dyn_mat<typename Derived::Scalar>::Identity(A.rows(), A.rows());
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
851  if (!internal::check_nonzero_size(rA))
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
887  if (!internal::check_nonzero_size(rA))
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
915 template <typename T>
916 dyn_mat<typename T::Scalar> kron(const T& head) {
917  return head;
918 }
919 
929 template <typename T, typename... Args>
930 dyn_mat<typename T::Scalar> kron(const T& head, const Args&... tail) {
931  return internal::kron2(head, kron(tail...));
932 }
933 
942 template <typename Derived>
943 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As) {
944  // EXCEPTION CHECKS
945 
946  if (As.empty())
947  throw exception::ZeroSize("qpp::kron()");
948 
949  for (auto&& elem : As)
950  if (!internal::check_nonzero_size(elem))
951  throw exception::ZeroSize("qpp::kron()");
952  // END EXCEPTION CHECKS
953 
954  dyn_mat<typename Derived::Scalar> result = As[0].derived();
955  for (idx i = 1; i < As.size(); ++i) {
956  result = kron(result, As[i]);
957  }
958 
959  return result;
960 }
961 
962 // Kronecker product of a list of matrices, preserve return type
963 // deduce the template parameters from initializer_list
973 template <typename Derived>
974 dyn_mat<typename Derived::Scalar>
975 kron(const std::initializer_list<Derived>& As) {
976  return kron(std::vector<Derived>(As));
977 }
978 
988 template <typename Derived>
989 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
990  idx n) {
991  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
992 
993  // EXCEPTION CHECKS
994 
995  // check zero-size
996  if (!internal::check_nonzero_size(rA))
997  throw exception::ZeroSize("qpp::kronpow()");
998 
999  // check out of range
1000  if (n == 0)
1001  throw exception::OutOfRange("qpp::kronpow()");
1002  // END EXCEPTION CHECKS
1003 
1004  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1005 
1006  return kron(As);
1007 }
1008 
1009 // Direct sum of multiple matrices, preserve return type
1010 // variadic template
1020 template <typename T>
1021 dyn_mat<typename T::Scalar> dirsum(const T& head) {
1022  return head;
1023 }
1024 
1034 template <typename T, typename... Args>
1035 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args&... tail) {
1036  return internal::dirsum2(head, dirsum(tail...));
1037 }
1038 
1047 template <typename Derived>
1048 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As) {
1049  // EXCEPTION CHECKS
1050 
1051  if (As.empty())
1052  throw exception::ZeroSize("qpp::dirsum()");
1053 
1054  for (auto&& elem : As)
1055  if (!internal::check_nonzero_size(elem))
1056  throw exception::ZeroSize("qpp::dirsum()");
1057  // END EXCEPTION CHECKS
1058 
1059  idx total_rows = 0, total_cols = 0;
1060  for (idx i = 0; i < As.size(); ++i) {
1061  total_rows += static_cast<idx>(As[i].rows());
1062  total_cols += static_cast<idx>(As[i].cols());
1063  }
1064  dyn_mat<typename Derived::Scalar> result =
1065  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1066 
1067  idx cur_row = 0, cur_col = 0;
1068  for (idx i = 0; i < As.size(); ++i) {
1069  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1070  cur_row += static_cast<idx>(As[i].rows());
1071  cur_col += static_cast<idx>(As[i].cols());
1072  }
1073 
1074  return result;
1075 }
1076 
1077 // Direct sum of a list of matrices, preserve return type
1078 // deduce the template parameters from initializer_list
1088 template <typename Derived>
1089 dyn_mat<typename Derived::Scalar>
1090 dirsum(const std::initializer_list<Derived>& As) {
1091  return dirsum(std::vector<Derived>(As));
1092 }
1093 
1103 template <typename Derived>
1104 dyn_mat<typename Derived::Scalar> dirsumpow(const Eigen::MatrixBase<Derived>& A,
1105  idx n) {
1106  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1107 
1108  // EXCEPTION CHECKS
1109 
1110  // check zero-size
1111  if (!internal::check_nonzero_size(rA))
1112  throw exception::ZeroSize("qpp::dirsumpow()");
1113 
1114  // check out of range
1115  if (n == 0)
1116  throw exception::OutOfRange("qpp::dirsumpow()");
1117  // END EXCEPTION CHECKS
1118 
1119  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1120 
1121  return dirsum(As);
1122 }
1123 
1135 template <typename Derived>
1136 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1137  idx rows, idx cols) {
1138  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1139 
1140  idx Arows = static_cast<idx>(rA.rows());
1141  idx Acols = static_cast<idx>(rA.cols());
1142 
1143  // EXCEPTION CHECKS
1144 
1145  // check zero-size
1146  if (!internal::check_nonzero_size(rA))
1147  throw exception::ZeroSize("qpp::reshape()");
1148 
1149  if (Arows * Acols != rows * cols)
1150  throw exception::DimsMismatchMatrix("qpp::reshape()");
1151  // END EXCEPTION CHECKS
1152 
1153  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1154  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1155 }
1156 
1169 template <typename Derived1, typename Derived2>
1170 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1171  const Eigen::MatrixBase<Derived2>& B) {
1172  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1173  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1174 
1175  // EXCEPTION CHECKS
1176 
1177  // check types
1178  if (!std::is_same<typename Derived1::Scalar,
1179  typename Derived2::Scalar>::value)
1180  throw exception::TypeMismatch("qpp::comm()");
1181 
1182  // check zero-size
1183  if (!internal::check_nonzero_size(rA) || !internal::check_nonzero_size(rB))
1184  throw exception::ZeroSize("qpp::comm()");
1185 
1186  // check square matrices
1187  if (!internal::check_square_mat(rA) || !internal::check_square_mat(rB))
1188  throw exception::MatrixNotSquare("qpp::comm()");
1189 
1190  // check equal dimensions
1191  if (rA.rows() != rB.rows())
1192  throw exception::DimsNotEqual("qpp::comm()");
1193  // END EXCEPTION CHECKS
1194 
1195  return rA * rB - rB * rA;
1196 }
1197 
1210 template <typename Derived1, typename Derived2>
1211 dyn_mat<typename Derived1::Scalar>
1212 anticomm(const Eigen::MatrixBase<Derived1>& A,
1213  const Eigen::MatrixBase<Derived2>& B) {
1214  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1215  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1216 
1217  // EXCEPTION CHECKS
1218 
1219  // check types
1220  if (!std::is_same<typename Derived1::Scalar,
1221  typename Derived2::Scalar>::value)
1222  throw exception::TypeMismatch("qpp::anticomm()");
1223 
1224  // check zero-size
1225  if (!internal::check_nonzero_size(rA) || !internal::check_nonzero_size(rB))
1226  throw exception::ZeroSize("qpp::anticomm()");
1227 
1228  // check square matrices
1229  if (!internal::check_square_mat(rA) || !internal::check_square_mat(rB))
1230  throw exception::MatrixNotSquare("qpp::anticomm()");
1231 
1232  // check equal dimensions
1233  if (rA.rows() != rB.rows())
1234  throw exception::DimsNotEqual("qpp::anticomm()");
1235  // END EXCEPTION CHECKS
1236 
1237  return rA * rB + rB * rA;
1238 }
1239 
1249 template <typename Derived>
1250 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& A) {
1251  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1252 
1253  // EXCEPTION CHECKS
1254 
1255  // check zero-size
1256  if (!internal::check_nonzero_size(rA))
1257  throw exception::ZeroSize("qpp::prj()");
1258 
1259  // check column vector
1260  if (!internal::check_cvector(rA))
1261  throw exception::MatrixNotCvector("qpp::prj()");
1262  // END EXCEPTION CHECKS
1263 
1264  double normA = norm(rA);
1265  if (normA > 0)
1266  return rA * adjoint(rA) / (normA * normA);
1267  else
1268  return dyn_mat<typename Derived::Scalar>::Zero(rA.rows(), rA.rows());
1269 }
1270 
1278 template <typename Derived>
1279 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& As) {
1280  // EXCEPTION CHECKS
1281 
1282  // check empty list
1283  if (!internal::check_nonzero_size(As))
1284  throw exception::ZeroSize("qpp::grams()");
1285 
1286  for (auto&& elem : As)
1287  if (!internal::check_nonzero_size(elem))
1288  throw exception::ZeroSize("qpp::grams()");
1289 
1290  // check that As[0] is a column vector
1291  if (!internal::check_cvector(As[0]))
1292  throw exception::MatrixNotCvector("qpp::grams()");
1293 
1294  // now check that all the rest match the size of the first vector
1295  for (auto&& elem : As)
1296  if (elem.rows() != As[0].rows() || elem.cols() != 1)
1297  throw exception::DimsNotEqual("qpp::grams()");
1298  // END EXCEPTION CHECKS
1299 
1300  dyn_mat<typename Derived::Scalar> cut =
1301  dyn_mat<typename Derived::Scalar>::Identity(As[0].rows(), As[0].rows());
1302 
1303  dyn_mat<typename Derived::Scalar> vi =
1304  dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1305 
1306  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1307  // find the first non-zero vector in the list
1308  idx pos = 0;
1309  for (pos = 0; pos < As.size(); ++pos) {
1310  if (norm(As[pos]) > 0) // add it as the first element
1311  {
1312  outvecs.emplace_back(As[pos]);
1313  break;
1314  }
1315  }
1316 
1317  // start the process
1318  for (idx i = pos + 1; i < As.size(); ++i) {
1319  cut -= prj(outvecs[i - 1 - pos]);
1320  vi = cut * As[i];
1321  outvecs.emplace_back(vi);
1322  }
1323 
1324  dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1325 
1326  idx cnt = 0;
1327  for (auto&& elem : outvecs) {
1328  double normA = norm(elem);
1329  if (normA > 0) // we add only the non-zero vectors
1330  {
1331  result.col(cnt) = elem / normA;
1332  ++cnt;
1333  }
1334  }
1335 
1336  return result.block(0, 0, As[0].rows(), cnt);
1337 }
1338 
1339 // deduce the template parameters from initializer_list
1347 template <typename Derived>
1348 dyn_mat<typename Derived::Scalar>
1349 grams(const std::initializer_list<Derived>& As) {
1350  return grams(std::vector<Derived>(As));
1351 }
1352 
1360 template <typename Derived>
1361 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A) {
1362  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1363 
1364  // EXCEPTION CHECKS
1365 
1366  if (!internal::check_nonzero_size(rA))
1367  throw exception::ZeroSize("qpp::grams()");
1368  // END EXCEPTION CHECKS
1369 
1370  std::vector<dyn_mat<typename Derived::Scalar>> input;
1371 
1372  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1373  input.emplace_back(rA.col(i));
1374 
1375  return grams<dyn_mat<typename Derived::Scalar>>(input);
1376 }
1377 
1388 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims) {
1389  // EXCEPTION CHECKS
1390 
1391  if (!internal::check_dims(dims))
1392  throw exception::DimsInvalid("qpp::n2multiidx()");
1393 
1394  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1395  static_cast<idx>(1), std::multiplies<idx>()))
1396  throw exception::OutOfRange("qpp::n2multiidx()");
1397  // END EXCEPTION CHECKS
1398 
1399  // double the size for matrices reshaped as vectors
1400  idx result[2 * maxn];
1401  internal::n2multiidx(n, dims.size(), dims.data(), result);
1402 
1403  return std::vector<idx>(result, result + dims.size());
1404 }
1405 
1416 inline idx multiidx2n(const std::vector<idx>& midx,
1417  const std::vector<idx>& dims) {
1418  // EXCEPTION CHECKS
1419 
1420  if (!internal::check_dims(dims))
1421  throw exception::DimsInvalid("qpp::multiidx2n()");
1422 
1423  for (idx i = 0; i < dims.size(); ++i)
1424  if (midx[i] >= dims[i]) {
1425  throw exception::OutOfRange("qpp::multiidx2n()");
1426  }
1427  // END EXCEPTION CHECKS
1428 
1429  return internal::multiidx2n(midx.data(), dims.size(), dims.data());
1430 }
1431 
1445 inline ket mket(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1446  idx n = mask.size();
1447 
1448  idx D = std::accumulate(std::begin(dims), std::end(dims),
1449  static_cast<idx>(1), std::multiplies<idx>());
1450 
1451  // EXCEPTION CHECKS
1452 
1453  // check zero size
1454  if (n == 0)
1455  throw exception::ZeroSize("qpp::mket()");
1456  // check valid dims
1457  if (!internal::check_dims(dims))
1458  throw exception::DimsInvalid("qpp::mket()");
1459  // check mask and dims have the same size
1460  if (mask.size() != dims.size())
1461  throw exception::SubsysMismatchDims("qpp::mket()");
1462  // check mask is a valid vector
1463  for (idx i = 0; i < n; ++i)
1464  if (mask[i] >= dims[i])
1465  throw exception::SubsysMismatchDims("qpp::mket()");
1466  // END EXCEPTION CHECKS
1467 
1468  ket result = ket::Zero(D);
1469  idx pos = multiidx2n(mask, dims);
1470  result(pos) = 1;
1471 
1472  return result;
1473 }
1474 
1488 inline ket mket(const std::vector<idx>& mask, idx d = 2) {
1489  idx n = mask.size();
1490  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1491 
1492  // EXCEPTION CHECKS
1493 
1494  // check zero size
1495  if (n == 0)
1496  throw exception::ZeroSize("qpp::mket()");
1497 
1498  // check valid dims
1499  if (d == 0)
1500  throw exception::DimsInvalid("qpp::mket()");
1501 
1502  // check mask is a valid vector
1503  for (idx i = 0; i < n; ++i)
1504  if (mask[i] >= d)
1505  throw exception::SubsysMismatchDims("qpp::mket()");
1506  // END EXCEPTION CHECKS
1507 
1508  ket result = ket::Zero(D);
1509  std::vector<idx> dims(n, d);
1510  idx pos = multiidx2n(mask, dims);
1511  result(pos) = 1;
1512 
1513  return result;
1514 }
1515 
1530 inline cmat mprj(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1531  idx n = mask.size();
1532 
1533  idx D = std::accumulate(std::begin(dims), std::end(dims),
1534  static_cast<idx>(1), std::multiplies<idx>());
1535 
1536  // EXCEPTION CHECKS
1537 
1538  // check zero size
1539  if (n == 0)
1540  throw exception::ZeroSize("qpp::mprj()");
1541  // check valid dims
1542  if (!internal::check_dims(dims))
1543  throw exception::DimsInvalid("qpp::mprj()");
1544  // check mask and dims have the same size
1545  if (mask.size() != dims.size())
1546  throw exception::SubsysMismatchDims("qpp::mprj()");
1547  // check mask is a valid vector
1548  for (idx i = 0; i < n; ++i)
1549  if (mask[i] >= dims[i])
1550  throw exception::SubsysMismatchDims("qpp::mprj()");
1551  // END EXCEPTION CHECKS
1552 
1553  cmat result = cmat::Zero(D, D);
1554  idx pos = multiidx2n(mask, dims);
1555  result(pos, pos) = 1;
1556 
1557  return result;
1558 }
1559 
1574 inline cmat mprj(const std::vector<idx>& mask, idx d = 2) {
1575  idx n = mask.size();
1576  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1577 
1578  // EXCEPTION CHECKS
1579 
1580  // check zero size
1581  if (n == 0)
1582  throw exception::ZeroSize("qpp::mprj()");
1583 
1584  // check valid dims
1585  if (d == 0)
1586  throw exception::DimsInvalid("qpp::mprj()");
1587 
1588  // check mask is a valid vector
1589  for (idx i = 0; i < n; ++i)
1590  if (mask[i] >= d)
1591  throw exception::SubsysMismatchDims("qpp::mprj()");
1592  // END EXCEPTION CHECKS
1593 
1594  cmat result = cmat::Zero(D, D);
1595  std::vector<idx> dims(n, d);
1596  idx pos = multiidx2n(mask, dims);
1597  result(pos, pos) = 1;
1598 
1599  return result;
1600 }
1601 
1610 template <typename InputIterator>
1611 std::vector<double> abssq(InputIterator first, InputIterator last) {
1612  std::vector<double> weights(std::distance(first, last));
1613  std::transform(first, last, std::begin(weights),
1614  [](cplx z) -> double { return std::norm(z); });
1615 
1616  return weights;
1617 }
1618 
1625 template <typename Container>
1626 std::vector<double>
1627 abssq(const Container& c,
1628  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1629 // we need the std::enable_if to SFINAE out Eigen expressions
1630 // that will otherwise match, instead of matching
1631 // the overload below:
1632 
1633 // template<typename Derived>
1634 // abssq(const Eigen::MatrixBase<Derived>& A)
1635 {
1636  return abssq(std::begin(c), std::end(c));
1637 }
1638 
1645 template <typename Derived>
1646 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A) {
1647  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1648 
1649  // EXCEPTION CHECKS
1650 
1651  // check zero-size
1652  if (!internal::check_nonzero_size(rA))
1653  throw exception::ZeroSize("qpp::abssq()");
1654  // END EXCEPTION CHECKS
1655 
1656  return abssq(rA.data(), rA.data() + rA.size());
1657 }
1658 
1667 template <typename InputIterator>
1668 typename std::iterator_traits<InputIterator>::value_type
1669 sum(InputIterator first, InputIterator last) {
1670  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1671 
1672  return std::accumulate(first, last, static_cast<value_type>(0));
1673 }
1674 
1682 template <typename Container>
1683 typename Container::value_type
1684 sum(const Container& c,
1685  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1686  return sum(std::begin(c), std::end(c));
1687 }
1688 
1697 template <typename InputIterator>
1698 typename std::iterator_traits<InputIterator>::value_type
1699 prod(InputIterator first, InputIterator last) {
1700  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1701 
1702  return std::accumulate(first, last, static_cast<value_type>(1),
1703  std::multiplies<value_type>());
1704 }
1705 
1713 template <typename Container>
1714 typename Container::value_type
1715 prod(const Container& c,
1716  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1717  return prod(std::begin(c), std::end(c));
1718 }
1719 
1732 template <typename Derived>
1733 dyn_col_vect<typename Derived::Scalar>
1734 rho2pure(const Eigen::MatrixBase<Derived>& A) {
1735  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1736 
1737  // EXCEPTION CHECKS
1738  // check zero-size
1739  if (!internal::check_nonzero_size(rA))
1740  throw exception::ZeroSize("qpp::rho2pure()");
1741  // check square matrix
1742  if (!internal::check_square_mat(rA))
1743  throw exception::MatrixNotSquare("qpp::rho2pure()");
1744  // END EXCEPTION CHECKS
1745 
1746  dyn_col_vect<double> tmp_evals = hevals(rA);
1747  cmat tmp_evects = hevects(rA);
1748  dyn_col_vect<typename Derived::Scalar> result =
1749  dyn_col_vect<typename Derived::Scalar>::Zero(rA.rows());
1750  // find the non-zero eigenvector
1751  // there is only one, assuming the state is pure
1752  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k) {
1753  if (std::abs(tmp_evals(k)) > 0) {
1754  result = tmp_evects.col(k);
1755  break;
1756  }
1757  }
1758 
1759  return result;
1760 }
1761 
1770 inline std::vector<idx> complement(std::vector<idx> subsys, idx n) {
1771  // EXCEPTION CHECKS
1772 
1773  if (n < subsys.size())
1774  throw exception::OutOfRange("qpp::complement()");
1775  for (idx i = 0; i < subsys.size(); ++i)
1776  if (subsys[i] >= n)
1777  throw exception::SubsysMismatchDims("qpp::complement()");
1778  // END EXCEPTION CHECKS
1779 
1780  std::vector<idx> all(n);
1781  std::vector<idx> subsys_bar(n - subsys.size());
1782 
1783  std::iota(std::begin(all), std::end(all), 0);
1784  std::sort(std::begin(subsys), std::end(subsys));
1785  std::set_difference(std::begin(all), std::end(all), std::begin(subsys),
1786  std::end(subsys), std::begin(subsys_bar));
1787 
1788  return subsys_bar;
1789 }
1790 
1801 template <typename Derived>
1802 std::vector<double> rho2bloch(const Eigen::MatrixBase<Derived>& A) {
1803  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1804 
1805  // EXCEPTION CHECKS
1806 
1807  // check qubit matrix
1808  if (!internal::check_qubit_matrix(rA))
1809  throw exception::NotQubitMatrix("qpp::rho2bloch()");
1810  // END EXCEPTION CHECKS
1811 
1812  std::vector<double> result(3);
1813  cmat X(2, 2), Y(2, 2), Z(2, 2);
1814  X << 0, 1, 1, 0;
1815  Y << 0, -1_i, 1_i, 0;
1816  Z << 1, 0, 0, -1;
1817  result[0] = std::real(trace(rA * X));
1818  result[1] = std::real(trace(rA * Y));
1819  result[2] = std::real(trace(rA * Z));
1820 
1821  return result;
1822 }
1823 
1832 inline cmat bloch2rho(const std::vector<double>& r) {
1833  // EXCEPTION CHECKS
1834 
1835  // check 3-dimensional vector
1836  if (r.size() != 3)
1837  throw exception::CustomException("qpp::bloch2rho",
1838  "r is not a 3-dimensional vector!");
1839  // END EXCEPTION CHECKS
1840 
1841  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1842  X << 0, 1, 1, 0;
1843  Y << 0, -1_i, 1_i, 0;
1844  Z << 1, 0, 0, -1;
1845  Id2 << 1, 0, 0, 1;
1846 
1847  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1848 }
1849 
1850 inline namespace literals {
1851 // Idea taken from http://techblog.altplus.co.jp/entry/2017/11/08/130921
1861 template <char... Bits>
1862 ket operator"" _ket() {
1863  constexpr idx n = sizeof...(Bits);
1864  constexpr char bits[n + 1] = {Bits..., '\0'};
1865  qpp::ket q = qpp::ket::Zero(static_cast<idx>(std::llround(std::pow(2, n))));
1866  idx pos = 0;
1867 
1868  // EXCEPTION CHECKS
1869 
1870  // check valid multi-partite qubit state
1871  for (idx i = 0; i < n; ++i) {
1872  if (bits[i] != '0' && bits[i] != '1')
1873  throw exception::OutOfRange(R"xxx(qpp::operator "" _ket())xxx");
1874  }
1875  // END EXCEPTION CHECKS
1876 
1877  pos = std::stoi(bits, nullptr, 2);
1878  q(pos) = 1;
1879 
1880  return q;
1881 }
1882 
1892 template <char... Bits>
1893 bra operator"" _bra() {
1894  constexpr idx n = sizeof...(Bits);
1895  constexpr char bits[n + 1] = {Bits..., '\0'};
1896  qpp::bra q = qpp::ket::Zero(static_cast<idx>(std::llround(std::pow(2, n))));
1897  idx pos = 0;
1898 
1899  // EXCEPTION CHECKS
1900 
1901  // check valid multi-partite qubit state
1902  for (idx i = 0; i < n; ++i) {
1903  if (bits[i] != '0' && bits[i] != '1')
1904  throw exception::OutOfRange(R"xxx(qpp::operator "" _bra())xxx");
1905  }
1906  // END EXCEPTION CHECKS
1907 
1908  pos = std::stoi(bits, nullptr, 2);
1909  q(pos) = 1;
1910 
1911  return q;
1912 }
1913 
1925 template <char... Bits>
1926 cmat operator"" _prj() {
1927  constexpr idx n = sizeof...(Bits);
1928  constexpr char bits[n + 1] = {Bits..., '\0'};
1929 
1930  // EXCEPTION CHECKS
1931 
1932  // check valid multi-partite qubit state
1933  for (idx i = 0; i < n; ++i) {
1934  if (bits[i] != '0' && bits[i] != '1')
1935  throw exception::OutOfRange(R"xxx(qpp::operator "" _prj())xxx");
1936  }
1937  // END EXCEPTION CHECKS
1938 
1939  return kron(operator""_ket<Bits...>(), operator""_bra<Bits...>());
1940 }
1941 } /* namespace literals */
1942 
1943 namespace internal {
1944 // hash combine, code taken from boost::hash_combine(), see
1945 // https://www.boost.org/doc/libs/1_69_0/doc/html/hash/reference.html#boost.hash_combine
1953 template <class T>
1954 void hash_combine(std::size_t& seed, const T& v) {
1955  std::hash<T> hasher;
1956  seed ^= hasher(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
1957 }
1958 } /* namespace internal */
1959 
1969 template <typename Derived>
1970 std::size_t hash_eigen(const Eigen::MatrixBase<Derived>& A,
1971  std::size_t seed = 0) {
1972  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1973 
1974  // EXCEPTION CHECKS
1975 
1976  // check zero-size
1977  if (!internal::check_nonzero_size(rA))
1978  throw exception::ZeroSize("qpp::hash_eigen()");
1979  // END EXCEPTION CHECKS
1980 
1981  auto* p = rA.data();
1982  idx sizeA = static_cast<idx>(rA.size());
1983  for (idx i = 0; i < sizeA; ++i) {
1984  internal::hash_combine(seed, std::real(p[i]));
1985  internal::hash_combine(seed, std::imag(p[i]));
1986  }
1987 
1988  return seed;
1989 }
1990 
1991 namespace internal {
1996 struct HashEigen {
1997  template <typename Derived>
1998  std::size_t operator()(const Eigen::MatrixBase<Derived>& A) const {
1999  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
2000  return hash_eigen(rA);
2001  }
2002 };
2003 
2010 struct EqualEigen {
2011  template <typename Derived>
2012  bool operator()(const Eigen::MatrixBase<Derived>& A,
2013  const Eigen::MatrixBase<Derived>& B) const {
2014  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
2015  const dyn_mat<typename Derived::Scalar>& rB = B.derived();
2016  if (rA.rows() == rB.rows() && rA.cols() == rB.cols())
2017  return rA == rB ? true : false;
2018  else
2019  return false;
2020  }
2021 };
2022 
2029 struct EqualSameSizeStringDits {
2030  bool operator()(const std::string& s1, const std::string& s2) const {
2031  std::vector<std::string> tk1, tk2;
2032  std::string w1, w2;
2033  std::stringstream ss1{s1}, ss2{s2};
2034 
2035  // tokenize the string into words (assumes words are separated by space)
2036  while (ss1 >> w1)
2037  tk1.emplace_back(w1);
2038  while (ss2 >> w2)
2039  tk2.emplace_back(w2);
2040 
2041  // compare lexicographically
2042  auto it1 = std::begin(tk1);
2043  auto it2 = std::begin(tk2);
2044  while (it1 != std::end(tk1) && it2 != std::end(tk2)) {
2045  auto n1 = std::stoll(*it1++);
2046  auto n2 = std::stoll(*it2++);
2047  if (n1 < n2)
2048  return true;
2049  else if (n1 == n2)
2050  continue;
2051  else if (n1 > n2)
2052  return false;
2053  }
2054  return false;
2055  }
2056 };
2057 
2058 } /* namespace internal */
2059 } /* namespace qpp */
2060 
2061 #endif /* FUNCTIONS_H_ */
qpp
Quantum++ main namespace.
Definition: circuits.h:35
Generated by   doxygen 1.8.11