32 #ifndef INTERNAL_UTIL_H_ 33 #define INTERNAL_UTIL_H_ 43 inline void n2multiidx(idx n, idx numdims,
const idx*
const dims,
44 idx* result) noexcept {
50 for (idx i = 0; i < numdims; ++i)
56 for (idx i = 0; i < numdims; ++i) {
57 result[numdims - i - 1] = n % (dims[numdims - i - 1]);
58 n /= (dims[numdims - i - 1]);
64 #if (__GNUC__ && !__clang__) 65 #pragma GCC diagnostic push 66 #pragma GCC diagnostic ignored "-Warray-bounds" 67 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" 71 inline idx multiidx2n(
const idx*
const midx, idx numdims,
72 const idx*
const dims) noexcept {
82 idx part_prod[2 * maxn];
85 part_prod[numdims - 1] = 1;
86 for (idx i = 1; i < numdims; ++i) {
87 part_prod[numdims - i - 1] = part_prod[numdims - i] * dims[numdims - i];
88 result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
91 return result + midx[numdims - 1];
93 #if (__GNUC__ && !__clang__) 94 #pragma GCC diagnostic pop 98 template <
typename Derived>
99 bool check_square_mat(
const Eigen::MatrixBase<Derived>& A) {
100 return A.rows() == A.cols();
104 template <
typename Derived>
105 bool check_vector(
const Eigen::MatrixBase<Derived>& A) {
106 return A.rows() == 1 || A.cols() == 1;
110 template <
typename Derived>
111 bool check_rvector(
const Eigen::MatrixBase<Derived>& A) {
112 return A.rows() == 1;
116 template <
typename Derived>
117 bool check_cvector(
const Eigen::MatrixBase<Derived>& A) {
118 return A.cols() == 1;
122 template <
typename T>
123 bool check_nonzero_size(
const T& x) noexcept {
124 return x.size() != 0;
128 template <
typename T1,
typename T2>
129 bool check_matching_sizes(
const T1& lhs,
const T2& rhs) noexcept {
130 return lhs.size() == rhs.size();
134 inline bool check_dims(
const std::vector<idx>& dims) {
138 return std::find_if(std::begin(dims), std::end(dims),
139 [dims](idx i) ->
bool {
144 }) == std::end(dims);
149 template <
typename Derived>
150 bool check_dims_match_mat(
const std::vector<idx>& dims,
151 const Eigen::MatrixBase<Derived>& A) {
154 assert(dims.size() > 0);
155 assert(A.rows() == A.cols());
157 idx proddim = std::accumulate(std::begin(dims), std::end(dims),
158 static_cast<idx>(1), std::multiplies<idx>());
160 return proddim ==
static_cast<idx
>(A.rows());
164 template <
typename Derived>
165 bool check_dims_match_cvect(
const std::vector<idx>& dims,
166 const Eigen::MatrixBase<Derived>& A) {
169 assert(dims.size() > 0);
170 assert(A.rows() > 0);
171 assert(A.cols() == 1);
173 idx proddim = std::accumulate(std::begin(dims), std::end(dims),
174 static_cast<idx>(1), std::multiplies<idx>());
176 return proddim ==
static_cast<idx
>(A.rows());
180 template <
typename Derived>
181 bool check_dims_match_rvect(
const std::vector<idx>& dims,
182 const Eigen::MatrixBase<Derived>& A) {
185 assert(dims.size() > 0);
186 assert(A.cols() > 0);
187 assert(A.rows() == 1);
189 idx proddim = std::accumulate(std::begin(dims), std::end(dims),
190 static_cast<idx>(1), std::multiplies<idx>());
193 return proddim ==
static_cast<idx
>(A.cols());
197 inline bool check_eq_dims(
const std::vector<idx>& dims, idx dim) noexcept {
200 assert(dims.size() > 0);
210 inline bool check_no_duplicates(std::vector<idx> v) {
211 std::sort(std::begin(v), std::end(v));
212 if (std::unique(std::begin(v), std::end(v)) != std::end(v))
219 inline bool check_subsys_match_dims(
const std::vector<idx>& subsys,
220 const std::vector<idx>& dims) {
224 if (subsys.size() > dims.size())
228 if (!check_no_duplicates(subsys))
232 return std::find_if(std::begin(subsys), std::end(subsys),
233 [dims](idx i) ->
bool {
234 return i + 1 > dims.size();
235 }) == std::end(subsys);
239 template <
typename Derived>
240 bool check_qubit_matrix(
const Eigen::MatrixBase<Derived>& A) noexcept {
241 return A.rows() == 2 && A.cols() == 2;
245 template <
typename Derived>
246 bool check_qubit_cvector(
const Eigen::MatrixBase<Derived>& A) noexcept {
247 return A.rows() == 2 && A.cols() == 1;
251 template <
typename Derived>
252 bool check_qubit_rvector(
const Eigen::MatrixBase<Derived>& A) noexcept {
253 return A.rows() == 1 && A.cols() == 2;
257 template <
typename Derived>
258 bool check_qubit_vector(
const Eigen::MatrixBase<Derived>& A) noexcept {
259 return (A.rows() == 1 && A.cols() == 2) || (A.rows() == 2 && A.cols() == 1);
263 inline bool check_perm(
const std::vector<idx>& perm) {
267 std::vector<idx> ordered(perm.size());
268 std::iota(std::begin(ordered), std::end(ordered), 0);
270 return std::is_permutation(std::begin(ordered), std::end(ordered),
276 template <
typename Derived1,
typename Derived2>
277 dyn_mat<typename Derived1::Scalar> kron2(
const Eigen::MatrixBase<Derived1>& A,
278 const Eigen::MatrixBase<Derived2>& B) {
279 const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
280 const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
285 if (!std::is_same<
typename Derived1::Scalar,
286 typename Derived2::Scalar>::value)
287 throw exception::TypeMismatch(
"qpp::kron()");
290 if (!internal::check_nonzero_size(rA))
291 throw exception::ZeroSize(
"qpp::kron()");
294 if (!internal::check_nonzero_size(rB))
295 throw exception::ZeroSize(
"qpp::kron()");
298 idx Acols =
static_cast<idx
>(rA.cols());
299 idx Arows =
static_cast<idx
>(rA.rows());
300 idx Bcols =
static_cast<idx
>(rB.cols());
301 idx Brows =
static_cast<idx
>(rB.rows());
303 dyn_mat<typename Derived1::Scalar> result;
304 result.resize(Arows * Brows, Acols * Bcols);
307 #pragma omp parallel for collapse(2) 308 #endif // WITH_OPENMP_ 310 for (idx j = 0; j < Acols; ++j)
311 for (idx i = 0; i < Arows; ++i)
312 result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
319 template <
typename Derived1,
typename Derived2>
320 dyn_mat<typename Derived1::Scalar>
321 dirsum2(
const Eigen::MatrixBase<Derived1>& A,
322 const Eigen::MatrixBase<Derived2>& B) {
323 const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
324 const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
329 if (!std::is_same<
typename Derived1::Scalar,
330 typename Derived2::Scalar>::value)
331 throw exception::TypeMismatch(
"qpp::dirsum()");
334 if (!internal::check_nonzero_size(rA))
335 throw exception::ZeroSize(
"qpp::dirsum()");
338 if (!internal::check_nonzero_size(rB))
339 throw exception::ZeroSize(
"qpp::dirsum()");
342 idx Acols =
static_cast<idx
>(rA.cols());
343 idx Arows =
static_cast<idx
>(rA.rows());
344 idx Bcols =
static_cast<idx
>(rB.cols());
345 idx Brows =
static_cast<idx
>(rB.rows());
347 dyn_mat<typename Derived1::Scalar> result =
348 dyn_mat<typename Derived1::Scalar>::Zero(Arows + Brows, Acols + Bcols);
350 result.block(0, 0, Arows, Acols) = rA;
351 result.block(Arows, Acols, Brows, Bcols) = rB;
357 template <
typename T>
359 void variadic_vector_emplace(std::vector<T>&) {}
362 template <
typename T,
typename First,
typename... Args>
363 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&&... args) {
364 v.emplace_back(std::forward<First>(first));
365 variadic_vector_emplace(v, std::forward<Args>(args)...);
370 inline idx get_num_subsys(idx D, idx d) {
376 return static_cast<idx
>(std::llround(std::log2(D) / std::log2(d)));
382 inline idx get_dim_subsys(idx sz, idx N) {
389 return static_cast<idx
>(std::llround(std::sqrt(sz)));
391 return static_cast<idx
>(std::llround(std::pow(sz, 1. / N)));
395 template <
typename T,
396 typename std::enable_if<std::numeric_limits<T>::is_iec559 ||
397 is_complex<T>::value>::type* =
nullptr>
398 T abs_chop(
const T& x,
double chop = qpp::chop) {
399 if (std::abs(x) < chop)
406 template <
typename T,
407 typename std::enable_if<!(std::numeric_limits<T>::is_iec559 ||
408 is_complex<T>::value)>::type* =
nullptr>
409 T abs_chop(
const T& x,
double QPP_UNUSED_ chop = qpp::chop) {
414 struct Display_Impl_ {
415 template <
typename T>
417 std::ostream& display_impl_(
const T& A, std::ostream& os,
418 double chop = qpp::chop)
const {
419 std::ostringstream ostr;
422 std::vector<std::string> vstr;
425 for (idx i = 0; i < static_cast<idx>(A.rows()); ++i) {
426 for (idx j = 0; j < static_cast<idx>(A.cols()); ++j) {
429 ostr.str(std::string{});
432 double re =
static_cast<cplx
>(A(i, j)).real();
433 double im =
static_cast<cplx
>(A(i, j)).imag();
435 if (std::abs(re) < chop && std::abs(im) < chop) {
439 vstr.emplace_back(ostr.str());
440 }
else if (std::abs(re) < chop) {
442 vstr.emplace_back(ostr.str() +
"i");
443 }
else if (std::abs(im) < chop) {
445 vstr.emplace_back(ostr.str() +
" ");
450 strA += (im > 0 ?
" + " :
" - ");
452 ostr.str(std::string());
453 ostr << std::abs(im);
456 vstr.emplace_back(strA);
462 std::vector<idx> maxlengthcols(A.cols(), 0);
464 for (idx i = 0; i < static_cast<idx>(A.rows()); ++i)
465 for (idx j = 0; j < static_cast<idx>(A.cols()); ++j)
466 if (vstr[i * A.cols() + j].size() > maxlengthcols[j])
467 maxlengthcols[j] = vstr[i * A.cols() + j].size();
470 for (idx i = 0; i < static_cast<idx>(A.rows()); ++i) {
471 os << std::setw(static_cast<int>(maxlengthcols[0])) << std::right
472 << vstr[i * A.cols()];
474 for (idx j = 1; j < static_cast<idx>(A.cols()); ++j)
475 os << std::setw(static_cast<int>(maxlengthcols[j] + 2))
476 << std::right << vstr[i * A.cols() + j];
478 if (i < static_cast<idx>(A.rows()) - 1)
Quantum++ main namespace.
Definition: circuits.h:35