53 template <
typename Derived1,
typename Derived2>
54 dyn_mat<typename Derived1::Scalar>
56 const Eigen::MatrixBase<Derived2>& A,
const std::vector<idx>& ctrl,
57 const std::vector<idx>& target,
const std::vector<idx>& dims) {
58 const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
65 if (!std::is_same<
typename Derived1::Scalar,
66 typename Derived2::Scalar>::value)
86 idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
87 for (
idx i = 1; i < ctrl.size(); ++i)
88 if (dims[ctrl[i]] != d)
100 std::vector<idx> target_dims(target.size());
101 for (
idx i = 0; i < target.size(); ++i)
102 target_dims[i] = dims[target[i]];
106 std::vector<idx> ctrlgate = ctrl;
107 ctrlgate.insert(std::end(ctrlgate), std::begin(target), std::end(target));
108 std::sort(std::begin(ctrlgate), std::end(ctrlgate));
117 std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
118 std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
119 for (
idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i) {
120 Ai.push_back(
powm(rA, i));
124 idx D =
static_cast<idx>(rstate.rows());
126 idx ctrlsize = ctrl.size();
127 idx ctrlgatesize = ctrlgate.size();
128 idx targetsize = target.size();
130 idx Dctrl =
static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
131 idx DA =
static_cast<idx>(rA.rows());
139 std::vector<idx> ctrlgate_bar =
complement(ctrlgate, n);
141 idx ctrlgate_barsize = ctrlgate_bar.size();
144 for (
idx i = 0; i < ctrlgate_barsize; ++i)
145 DCTRLA_bar *= dims[ctrlgate_bar[i]];
147 for (
idx k = 0; k < n; ++k)
149 for (
idx k = 0; k < targetsize; ++k)
150 CdimsA[k] = dims[target[k]];
151 for (
idx k = 0; k < ctrlsize; ++k)
153 for (
idx k = 0; k < ctrlgate_barsize; ++k)
154 CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
158 auto coeff_idx_ket = [&](
idx i_,
idx m_,
idx r_) noexcept
159 ->std::pair<typename Derived1::Scalar, idx> {
161 typename Derived1::Scalar coeff = 0;
170 for (
idx k = 0; k < ctrlsize; ++k) {
177 for (
idx k = 0; k < n - ctrlgatesize; ++k) {
178 Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
183 for (
idx k = 0; k < targetsize; ++k) {
184 Cmidx[target[k]] = CmidxA[k];
191 for (
idx n_ = 0; n_ < DA; ++n_) {
193 for (
idx k = 0; k < targetsize; ++k) {
194 Cmidx[target[k]] = CmidxA[k];
200 return std::make_pair(coeff, indx);
208 ->std::tuple<typename Derived1::Scalar, idx, idx> {
211 typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
228 for (
idx k = 0; k < ctrlsize; ++k) {
229 Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
230 Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
238 for (
idx k = 0; k < n - ctrlgatesize; ++k) {
239 Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
240 Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
246 for (
idx k = 0; k < target.size(); ++k) {
247 Cmidxrow[target[k]] = CmidxArow[k];
248 Cmidxcol[target[k]] = CmidxAcol[k];
256 bool all_ctrl_rows_equal =
true;
257 bool all_ctrl_cols_equal =
true;
259 idx first_ctrl_row, first_ctrl_col;
261 first_ctrl_row = CmidxCTRLrow[0];
262 first_ctrl_col = CmidxCTRLcol[0];
264 first_ctrl_row = first_ctrl_col = 1;
267 for (
idx k = 1; k < ctrlsize; ++k) {
268 if (CmidxCTRLrow[k] != first_ctrl_row) {
269 all_ctrl_rows_equal =
false;
273 for (
idx k = 1; k < ctrlsize; ++k) {
274 if (CmidxCTRLcol[k] != first_ctrl_col) {
275 all_ctrl_cols_equal =
false;
281 for (
idx n1_ = 0; n1_ < DA; ++n1_) {
283 for (
idx k = 0; k < targetsize; ++k) {
284 Cmidxrow[target[k]] = CmidxArow[k];
288 if (all_ctrl_rows_equal) {
289 lhs = Ai[first_ctrl_row](m1_, n1_);
291 lhs = (m1_ == n1_) ? 1 : 0;
294 for (
idx n2_ = 0; n2_ < DA; ++n2_) {
296 for (
idx k = 0; k < targetsize; ++k) {
297 Cmidxcol[target[k]] = CmidxAcol[k];
300 if (all_ctrl_cols_equal) {
301 rhs = Aidagger[first_ctrl_col](n2_, m2_);
303 rhs = (n2_ == m2_) ? 1 : 0;
308 coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
312 return std::make_tuple(coeff, idxrow, idxcol);
327 #pragma omp parallel for collapse(2) 328 #endif // WITH_OPENMP_ 329 for (
idx m = 0; m < DA; ++m)
330 for (
idx r = 0; r < DCTRLA_bar; ++r) {
333 result(coeff_idx_ket(1, m, r).second) =
334 coeff_idx_ket(1, m, r).first;
336 for (
idx i = 0; i < d; ++i) {
337 result(coeff_idx_ket(i, m, r).second) =
338 coeff_idx_ket(i, m, r).first;
357 #pragma omp parallel for collapse(4) 358 #endif // WITH_OPENMP_ 359 for (
idx m1 = 0; m1 < DA; ++m1)
360 for (
idx r1 = 0; r1 < DCTRLA_bar; ++r1)
361 for (
idx m2 = 0; m2 < DA; ++m2)
362 for (
idx r2 = 0; r2 < DCTRLA_bar; ++r2)
366 coeff_idx_rho(1, m1, r1, 1, m2, r2);
367 result(std::get<1>(coeff_idxes),
368 std::get<2>(coeff_idxes)) =
369 std::get<0>(coeff_idxes);
371 for (
idx i1 = 0; i1 < Dctrl; ++i1)
372 for (
idx i2 = 0; i2 < Dctrl; ++i2) {
374 coeff_idx_rho(i1, m1, r1, i2, m2, r2);
375 result(std::get<1>(coeff_idxes),
376 std::get<2>(coeff_idxes)) =
377 std::get<0>(coeff_idxes);
403 template <
typename Derived1,
typename Derived2>
404 dyn_mat<typename Derived1::Scalar>
406 const Eigen::MatrixBase<Derived2>& A,
const std::vector<idx>& ctrl,
407 const std::vector<idx>& target,
idx d = 2) {
408 const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
424 std::vector<idx> dims(n, d);
426 return applyCTRL(rstate, rA, ctrl, target, dims);
442 template <
typename Derived1,
typename Derived2>
443 dyn_mat<typename Derived1::Scalar>
444 apply(
const Eigen::MatrixBase<Derived1>& state,
445 const Eigen::MatrixBase<Derived2>& A,
const std::vector<idx>& target,
446 const std::vector<idx>& dims) {
447 const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
454 if (!std::is_same<
typename Derived1::Scalar,
455 typename Derived2::Scalar>::value)
483 std::vector<idx> subsys_dims(target.size());
484 for (
idx i = 0; i < target.size(); ++i)
485 subsys_dims[i] = dims[target[i]];
497 return applyCTRL(rstate, rA, {}, target, dims);
506 return applyCTRL(rstate, rA, {}, target, dims);
526 template <
typename Derived1,
typename Derived2>
527 dyn_mat<typename Derived1::Scalar>
528 apply(
const Eigen::MatrixBase<Derived1>& state,
529 const Eigen::MatrixBase<Derived2>& A,
const std::vector<idx>& target,
531 const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
547 std::vector<idx> dims(n, d);
549 return apply(rstate, rA, target, dims);
560 template <
typename Derived>
561 cmat apply(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks) {
562 const cmat& rA = A.derived();
574 if (Ks[0].rows() != rA.rows())
576 for (
auto&& elem : Ks)
577 if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
581 cmat result = cmat::Zero(rA.rows(), rA.rows());
584 #pragma omp parallel for 585 #endif // WITH_OPENMP_ 586 for (
idx i = 0; i < Ks.size(); ++i) {
589 #endif // WITH_OPENMP_ 590 { result += Ks[i] * rA *
adjoint(Ks[i]); }
606 template <
typename Derived>
607 cmat apply(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks,
608 const std::vector<idx>& target,
const std::vector<idx>& dims) {
609 const cmat& rA = A.derived();
637 std::vector<idx> subsys_dims(target.size());
638 for (
idx i = 0; i < target.size(); ++i)
639 subsys_dims[i] = dims[target[i]];
648 for (
auto&& elem : Ks)
649 if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
653 cmat result = cmat::Zero(rA.rows(), rA.rows());
655 for (
idx i = 0; i < Ks.size(); ++i)
656 result +=
apply(rA, Ks[i], target, dims);
671 template <
typename Derived>
672 cmat apply(
const Eigen::MatrixBase<Derived>& A,
const std::vector<cmat>& Ks,
673 const std::vector<idx>& target,
idx d = 2) {
674 const cmat& rA = A.derived();
688 std::vector<idx> dims(n, d);
690 return apply(rA, Ks, target, dims);
713 for (
auto&& elem : Ks)
714 if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
718 idx D =
static_cast<idx>(Ks[0].rows());
720 cmat result(D * D, D * D);
721 cmat MN = cmat::Zero(D, D);
722 bra A = bra::Zero(D);
723 ket B = ket::Zero(D);
724 cmat EMN = cmat::Zero(D, D);
727 #pragma omp parallel for collapse(2) 728 #endif // WITH_OPENMP_ 729 for (
idx m = 0; m < D; ++m) {
730 for (
idx n = 0; n < D; ++n) {
733 #endif // WITH_OPENMP_ 737 for (
idx i = 0; i < Ks.size(); ++i)
738 EMN += Ks[i] * MN *
adjoint(Ks[i]);
741 for (
idx a = 0; a < D; ++a) {
743 for (
idx b = 0; b < D; ++b) {
746 result(a * D + b, m * D + n) =
747 static_cast<cmat>(A * EMN * B).value();
752 EMN = cmat::Zero(D, D);
784 for (
auto&& elem : Ks)
785 if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
789 idx D =
static_cast<idx>(Ks[0].rows());
793 cmat MES = cmat::Zero(D * D, 1);
794 for (
idx a = 0; a < D; ++a)
799 cmat result = cmat::Zero(D * D, D * D);
802 #pragma omp parallel for 803 #endif // WITH_OPENMP_ 804 for (
idx i = 0; i < Ks.size(); ++i) {
807 #endif // WITH_OPENMP_ 809 result +=
kron(cmat::Identity(D, D), Ks[i]) * Omega *
839 if (D * D != static_cast<idx>(A.rows()))
845 std::vector<cmat> result;
847 for (
idx i = 0; i < D * D; ++i) {
848 if (std::abs(ev(i)) > 0)
849 result.push_back(std::sqrt(std::abs(ev(i))) *
872 if (D * D != static_cast<idx>(A.rows()))
876 cmat result(D * D, D * D);
879 #pragma omp parallel for collapse(4) 880 #endif // WITH_OPENMP_ 881 for (
idx a = 0; a < D; ++a)
882 for (
idx b = 0; b < D; ++b)
883 for (
idx m = 0; m < D; ++m)
884 for (
idx n = 0; n < D; ++n)
885 result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
906 if (D * D != static_cast<idx>(A.rows()))
910 cmat result(D * D, D * D);
913 #pragma omp parallel for collapse(4) 914 #endif // WITH_OPENMP_ 915 for (
idx a = 0; a < D; ++a)
916 for (
idx b = 0; b < D; ++b)
917 for (
idx m = 0; m < D; ++m)
918 for (
idx n = 0; n < D; ++n)
919 result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
937 template <
typename Derived>
939 const std::vector<idx>& dims) {
940 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
953 if (dims.size() != 2)
970 auto worker = [&](
idx i,
idx j) noexcept->typename Derived::Scalar {
971 typename Derived::Scalar
sum = 0;
972 for (
idx m = 0; m < DA; ++m)
973 sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
979 #pragma omp parallel for collapse(2) 980 #endif // WITH_OPENMP_ 982 for (
idx j = 0; j < DB; ++j)
983 for (
idx i = 0; i < DB; ++i)
984 result(i, j) = worker(i, j);
993 throw exception::DimsMismatchMatrix(
"qpp::ptrace1()");
995 auto worker = [&](
idx i,
idx j) noexcept->typename Derived::Scalar {
996 typename Derived::Scalar
sum = 0;
997 for (
idx m = 0; m < DA; ++m)
998 sum += rA(m * DB + i, m * DB + j);
1004 #pragma omp parallel for collapse(2) 1005 #endif // WITH_OPENMP_ 1007 for (
idx j = 0; j < DB; ++j)
1008 for (
idx i = 0; i < DB; ++i)
1009 result(i, j) = worker(i, j);
1015 throw exception::MatrixNotSquareNorCvector(
"qpp::ptrace1()");
1031 template <
typename Derived>
1047 std::vector<idx> dims(2, d);
1065 template <
typename Derived>
1067 const std::vector<idx>& dims) {
1068 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1081 if (dims.size() != 2)
1098 auto worker = [&](
idx i,
idx j) noexcept->typename Derived::Scalar {
1099 typename Derived::Scalar
sum = 0;
1100 for (
idx m = 0; m < DB; ++m)
1101 sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1107 #pragma omp parallel for collapse(2) 1108 #endif // WITH_OPENMP_ 1110 for (
idx j = 0; j < DA; ++j)
1111 for (
idx i = 0; i < DA; ++i)
1112 result(i, j) = worker(i, j);
1121 throw exception::DimsMismatchMatrix(
"qpp::ptrace2()");
1124 #pragma omp parallel for collapse(2) 1125 #endif // WITH_OPENMP_ 1127 for (
idx j = 0; j < DA; ++j)
1128 for (
idx i = 0; i < DA; ++i)
1129 result(i, j) =
trace(rA.block(i * DB, j * DB, DB, DB));
1135 throw exception::MatrixNotSquareNorCvector(
"qpp::ptrace1()");
1151 template <
typename Derived>
1167 std::vector<idx> dims(2, d);
1186 template <
typename Derived>
1188 const std::vector<idx>& target,
1189 const std::vector<idx>& dims) {
1190 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1207 idx D =
static_cast<idx>(rA.rows());
1208 idx n = dims.size();
1209 idx n_subsys = target.size();
1210 idx n_subsys_bar = n - n_subsys;
1212 for (
idx i = 0; i < n_subsys; ++i)
1213 Dsubsys *= dims[target[i]];
1214 idx Dsubsys_bar = D / Dsubsys;
1224 std::vector<idx> subsys_bar =
complement(target, n);
1225 std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1226 std::begin(Csubsys_bar));
1228 for (
idx i = 0; i < n; ++i) {
1231 for (
idx i = 0; i < n_subsys; ++i) {
1232 Csubsys[i] = target[i];
1233 Cdimssubsys[i] = dims[target[i]];
1235 for (
idx i = 0; i < n_subsys_bar; ++i) {
1236 Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1249 if (target.size() == dims.size()) {
1250 result(0, 0) = (
adjoint(rA) * rA).value();
1254 if (target.size() == 0)
1257 auto worker = [&](
idx i) noexcept->typename Derived::Scalar {
1267 Cmidxrowsubsys_bar);
1269 for (
idx k = 0; k < n_subsys_bar; ++k) {
1270 Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1271 Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1273 typename Derived::Scalar sm = 0;
1274 for (
idx a = 0; a < Dsubsys; ++a) {
1278 for (
idx k = 0; k < n_subsys; ++k)
1279 Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1290 for (
idx j = 0; j < Dsubsys_bar; ++j)
1294 Cmidxcolsubsys_bar);
1296 #pragma omp parallel for 1297 #endif // WITH_OPENMP_ 1298 for (
idx i = 0; i < Dsubsys_bar; ++i) {
1299 result(i, j) = worker(i);
1310 throw exception::DimsMismatchMatrix(
"qpp::ptrace()");
1312 if (target.size() == dims.size()) {
1313 result(0, 0) = rA.trace();
1317 if (target.size() == 0)
1320 auto worker = [&](
idx i) noexcept->typename Derived::Scalar {
1330 Cmidxrowsubsys_bar);
1332 for (
idx k = 0; k < n_subsys_bar; ++k) {
1333 Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1334 Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1336 typename Derived::Scalar sm = 0;
1337 for (
idx a = 0; a < Dsubsys; ++a) {
1341 for (
idx k = 0; k < n_subsys; ++k)
1342 Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1353 for (
idx j = 0; j < Dsubsys_bar; ++j)
1357 Cmidxcolsubsys_bar);
1359 #pragma omp parallel for 1360 #endif // WITH_OPENMP_ 1361 for (
idx i = 0; i < Dsubsys_bar; ++i) {
1362 result(i, j) = worker(i);
1370 throw exception::MatrixNotSquareNorCvector(
"qpp::ptrace()");
1387 template <
typename Derived>
1389 const std::vector<idx>& target,
1391 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1405 std::vector<idx> dims(n, d);
1407 return ptrace(rA, target, dims);
1423 template <
typename Derived>
1424 dyn_mat<typename Derived::Scalar>
1425 ptranspose(
const Eigen::MatrixBase<Derived>& A,
const std::vector<idx>& target,
1426 const std::vector<idx>& dims) {
1427 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1444 idx D =
static_cast<idx>(rA.rows());
1445 idx n = dims.size();
1446 idx n_subsys = target.size();
1452 for (
idx i = 0; i < n; ++i)
1454 for (
idx i = 0; i < n_subsys; ++i)
1455 Csubsys[i] = target[i];
1466 if (target.size() == dims.size())
1469 if (target.size() == 0)
1472 auto worker = [&](
idx i) noexcept->typename Derived::Scalar {
1477 for (
idx k = 0; k < n; ++k)
1478 midxcoltmp[k] = Cmidxcol[k];
1483 for (
idx k = 0; k < n_subsys; ++k)
1484 std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1491 for (
idx j = 0; j < D; ++j) {
1496 #pragma omp parallel for 1497 #endif // WITH_OPENMP_ 1498 for (
idx i = 0; i < D; ++i)
1499 result(i, j) = worker(i);
1509 throw exception::DimsMismatchMatrix(
"qpp::ptranspose()");
1511 if (target.size() == dims.size())
1512 return rA.transpose();
1514 if (target.size() == 0)
1517 auto worker = [&](
idx i) noexcept->typename Derived::Scalar {
1522 for (
idx k = 0; k < n; ++k)
1523 midxcoltmp[k] = Cmidxcol[k];
1528 for (
idx k = 0; k < n_subsys; ++k)
1529 std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1536 for (
idx j = 0; j < D; ++j) {
1541 #pragma omp parallel for 1542 #endif // WITH_OPENMP_ 1543 for (
idx i = 0; i < D; ++i)
1544 result(i, j) = worker(i);
1551 throw exception::MatrixNotSquareNorCvector(
"qpp::ptranspose()");
1567 template <
typename Derived>
1568 dyn_mat<typename Derived::Scalar>
1569 ptranspose(
const Eigen::MatrixBase<Derived>& A,
const std::vector<idx>& target,
1571 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1585 std::vector<idx> dims(n, d);
1602 template <
typename Derived>
1603 dyn_mat<typename Derived::Scalar>
1604 syspermute(
const Eigen::MatrixBase<Derived>& A,
const std::vector<idx>& perm,
1605 const std::vector<idx>& dims) {
1606 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1623 if (perm.size() != dims.size())
1627 idx D =
static_cast<idx>(rA.rows());
1628 idx n = dims.size();
1643 for (
idx i = 0; i < n; ++i) {
1647 result.resize(D, 1);
1649 auto worker = [&Cdims, &Cperm, n ](
idx i) noexcept->idx {
1659 for (
idx k = 0; k < n; ++k) {
1660 permdims[k] = Cdims[Cperm[k]];
1661 midxtmp[k] = midx[Cperm[k]];
1668 #pragma omp parallel for 1669 #endif // WITH_OPENMP_ 1670 for (
idx i = 0; i < D; ++i)
1671 result(worker(i)) = rA(i);
1683 throw exception::DimsMismatchMatrix(
"qpp::syspermute()");
1686 for (
idx i = 0; i < n; ++i) {
1688 Cdims[i + n] = dims[i];
1690 Cperm[i + n] = perm[i] + n;
1692 result.resize(D * D, 1);
1694 dyn_mat<typename Derived::Scalar> vectA =
1695 Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1696 const_cast<typename Derived::Scalar*
>(rA.data()), D * D, 1);
1698 auto worker = [&Cdims, &Cperm, n ](
idx i) noexcept->idx {
1708 for (
idx k = 0; k < 2 * n; ++k) {
1709 permdims[k] = Cdims[Cperm[k]];
1710 midxtmp[k] = midx[Cperm[k]];
1717 #pragma omp parallel for 1718 #endif // WITH_OPENMP_ 1719 for (
idx i = 0; i < D * D; ++i)
1720 result(worker(i)) = rA(i);
1726 throw exception::MatrixNotSquareNorCvector(
"qpp::syspermute()");
1741 template <
typename Derived>
1742 dyn_mat<typename Derived::Scalar>
1743 syspermute(
const Eigen::MatrixBase<Derived>& A,
const std::vector<idx>& perm,
1745 const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1759 std::vector<idx> dims(n, d);
1776 template <
typename Derived>
1778 const std::vector<idx>& target,
1779 idx d = 2,
bool swap =
true) {
1795 std::vector<idx> dims(n, d);
1821 idx n_subsys = target.size();
1825 for (
idx i = 0; i < n_subsys; ++i) {
1827 result =
apply(result, Gates::get_instance().H, {target[i]});
1829 for (
idx j = 2; j <= n_subsys - i; ++j) {
1832 Rj << 1, 0, 0, exp(2.0 *
pi * 1_i / std::pow(2, j));
1834 applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]});
1839 for (
idx i = 0; i < n_subsys / 2; ++i) {
1840 result =
apply(result, Gates::get_instance().SWAP,
1841 {target[i], target[n_subsys - i - 1]});
1846 for (
idx i = 0; i < n_subsys; ++i) {
1848 result =
apply(result, Gates::get_instance().Fd(d), {target[i]}, d);
1850 for (
idx j = 2; j <= n_subsys - i; ++j) {
1852 cmat Rj = cmat::Zero(d, d);
1853 for (
idx m = 0; m < d; ++m) {
1854 Rj(m, m) = exp(2.0 *
pi * m * 1_i / std::pow(d, j));
1857 applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]}, d);
1862 for (
idx i = 0; i < n_subsys / 2; ++i) {
1863 result =
apply(result, Gates::get_instance().SWAPd(d),
1864 {target[i], target[n_subsys - i - 1]}, d);
1884 template <
typename Derived>
1886 const std::vector<idx>& target,
1887 idx d = 2,
bool swap =
true) {
1903 std::vector<idx> dims(n, d);
1929 idx n_subsys = target.size();
1935 for (
idx i = 0; i < n_subsys / 2; ++i) {
1936 result =
apply(result, Gates::get_instance().SWAP,
1937 {target[i], target[n_subsys - i - 1]});
1940 for (
idx i = n_subsys; i-- > 0;) {
1942 for (
idx j = n_subsys - i + 1; j-- > 2;) {
1945 Rj << 1, 0, 0, exp(-2.0 *
pi * 1_i / std::pow(2, j));
1947 applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]});
1950 result =
apply(result, Gates::get_instance().H, {target[i]});
1955 for (
idx i = 0; i < n_subsys / 2; ++i) {
1956 result =
apply(result, Gates::get_instance().SWAPd(d),
1957 {target[i], target[n_subsys - i - 1]}, d);
1960 for (
idx i = n_subsys; i-- > 0;) {
1962 for (
idx j = n_subsys - i + 1; j-- > 2;) {
1964 cmat Rj = cmat::Zero(d, d);
1965 for (
idx m = 0; m < d; ++m) {
1966 Rj(m, m) = exp(-2.0 *
pi * m * 1_i / std::pow(d, j));
1969 applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]}, d);
1972 result =
apply(result,
adjoint(Gates::get_instance().Fd(d)),
1989 template <
typename Derived>
1991 idx d = 2,
bool swap =
true) {
2007 std::vector<idx> dims(n, d);
2028 std::vector<idx> subsys(n);
2029 std::iota(std::begin(subsys), std::end(subsys), 0);
2044 template <
typename Derived>
2046 idx d = 2,
bool swap =
true) {
2062 std::vector<idx> dims(n, d);
2083 std::vector<idx> subsys(n);
2084 std::iota(std::begin(subsys), std::end(subsys), 0);
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
Dimensions not equal exception.
Definition: exception.h:284
Dimension(s) mismatch matrix size exception.
Definition: exception.h:300
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:67
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:219
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:830
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:69
dyn_mat< typename Derived::Scalar > applyTFQ(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2, bool swap=true)
Applies the inverse (adjoint) qudit quantum Fourier transform to the part target of the multi-partite...
Definition: operations.h:1885
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:59
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:150
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
dyn_mat< typename Derived::Scalar > syspermute(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &perm, idx d=2)
Subsystem permutation.
Definition: operations.h:1743
Subsystems mismatch dimensions exception.
Definition: exception.h:365
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Quantum++ main namespace.
Definition: circuits.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:263
Matrix is not square nor column vector exception.
Definition: exception.h:209
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:425
Invalid dimension(s) exception.
Definition: exception.h:269
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:863
idx get_num_subsys(idx D, idx d)
Definition: util.h:370
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:897
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:165
Not bi-partite exception.
Definition: exception.h:486
dyn_mat< typename Derived1::Scalar > apply(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &target, const std::vector< idx > &dims)
Applies the gate A to the part target of the multi-partite state vector or density matrix state...
Definition: operations.h:444
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:938
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1396
Invalid permutation exception.
Definition: exception.h:380
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:89
dyn_mat< typename Derived1::Scalar > applyCTRL(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &ctrl, const std::vector< idx > &target, const std::vector< idx > &dims)
Applies the controlled-gate A to the part target of the multi-partite state vector or density matrix ...
Definition: operations.h:55
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1424
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:917
dyn_mat< typename Derived1::Scalar > applyCTRL(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &ctrl, const std::vector< idx > &target, idx d=2)
Applies the controlled-gate A to the part target of the multi-partite state vector or density matrix ...
Definition: operations.h:405
Matrix mismatch subsystems exception.
Definition: exception.h:254
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, idx d=2)
Partial trace.
Definition: operations.h:1152
Dimension(s) mismatch column vector size exception.
Definition: exception.h:316
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:46
dyn_col_vect< typename Derived::Scalar > TFQ(const Eigen::MatrixBase< Derived > &A, idx d=2, bool swap=true)
Inverse (adjoint) qudit quantum Fourier transform.
Definition: operations.h:1990
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
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:382
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors of Hermitian matrix.
Definition: functions.h:450
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:775
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
dyn_mat< typename Derived::Scalar > ptrace(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2)
Partial trace.
Definition: operations.h:1388
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:130
constexpr double pi
Definition: constants.h:72
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:704
dyn_mat< typename Derived::Scalar > applyQFT(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2, bool swap=true)
Applies the qudit quantum Fourier transform to the part target of the multi-partite state vector or d...
Definition: operations.h:1777
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
cmat apply(const Eigen::MatrixBase< Derived > &A, const std::vector< cmat > &Ks, const std::vector< idx > &target, idx d=2)
Applies the channel specified by the set of Kraus operators Ks to the part target of the multi-partit...
Definition: operations.h:672
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
Permutation mismatch dimensions exception.
Definition: exception.h:396
dyn_mat< typename Derived::Scalar > ptranspose(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2)
Partial transpose.
Definition: operations.h:1569
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, idx d=2)
Partial trace.
Definition: operations.h:1032
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
dyn_col_vect< typename Derived::Scalar > QFT(const Eigen::MatrixBase< Derived > &A, idx d=2, bool swap=true)
Qudit quantum Fourier transform.
Definition: operations.h:2045
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
Object has zero size exception.
Definition: exception.h:134