From 7f2c83db624f7e6c0ea629a8b12b6b63a550f8b1 Mon Sep 17 00:00:00 2001 From: boschmitt <7152025+boschmitt@users.noreply.github.com> Date: Wed, 9 Oct 2024 12:39:51 +0200 Subject: [PATCH 1/5] Fix multiplication phase correction issue Signed-off-by: boschmitt <7152025+boschmitt@users.noreply.github.com> --- runtime/cudaq/spin/spin_op.cpp | 75 ++++++++++++++++++---------------- 1 file changed, 39 insertions(+), 36 deletions(-) diff --git a/runtime/cudaq/spin/spin_op.cpp b/runtime/cudaq/spin/spin_op.cpp index f1acaf1d70..3313bde0c6 100644 --- a/runtime/cudaq/spin/spin_op.cpp +++ b/runtime/cudaq/spin/spin_op.cpp @@ -63,47 +63,50 @@ actionOnBra(spin_op &term, const std::string &bitConfiguration) { } std::pair, std::vector> -mult(std::vector row, std::vector other_row, +mult(const std::vector& p1, const std::vector p2, std::complex &rowCoeff, std::complex &otherCoeff) { - // This is term_i * otherTerm_j - std::vector tmp(row.size()), tmp2(row.size()); - std::size_t numQubits = row.size() / 2; - - for (std::size_t i = 0; i < 2 * numQubits; i++) - tmp[i] = row[i] ^ other_row[i]; + std::size_t numQubits = p1.size() / 2; + std::vector result(p1.size(), false); + int phase = 0; + + int p1Phase = 0; + int p2Phase = 0; + int resultPhase = 0; + int cPhase = 0; + + // Calculate the result of teh Pauli multiplication and the phase shift + for (std::size_t i = 0; i < numQubits; ++i) { + bool p1_x = p1[i]; + bool p1_z = p1[i + numQubits]; + bool p2_x = p2[i]; + bool p2_z = p2[i + numQubits]; + + // Compute the resulting Pauli operator + result[i] = p1_x ^ p2_x; + result[i + numQubits] = p1_z ^ p2_z; + + p1Phase += p1_x & p1_z; + p2Phase += p2_x & p2_z; + cPhase += p1_x & p2_z; + resultPhase += result[i] & result[i + numQubits]; + } + phase = 2*cPhase + p1Phase + p2Phase - resultPhase; - for (std::size_t i = 0; i < numQubits; i++) - tmp2[i] = (row[i] && other_row[numQubits + i]) || - (row[i + numQubits] && other_row[i]); + // Normalize the phase to a value in the range [0, 3] + phase %= 4; + if (phase < 0) + phase += 4; - int orig_phase = 0, other_phase = 0; - for (std::size_t i = 0; i < numQubits; i++) { - if (row[i] && row[i + numQubits]) - orig_phase++; + // Phase correction factors based on the total phase + using namespace std::complex_literals; + std::array, 4> phaseCoeffArr{1.0, -1i, -1.0, 1i}; + std::complex final_coeff = rowCoeff * phaseCoeffArr[phase] * otherCoeff; - if (other_row[i] && other_row[i + numQubits]) - other_phase++; - } + // Handle the "-0" issue + if (std::abs(final_coeff.real()) < 1e-12) + final_coeff.real(0); - int sum = 0; - for (auto a : tmp2) - if (a) - sum++; - - auto _phase = orig_phase + other_phase + 2 * sum; - // Based on the phase, figure out an extra coeff to apply - for (std::size_t i = 0; i < numQubits; i++) - if (tmp[i] && tmp[i + numQubits]) - _phase -= 1; - - _phase %= 4; - std::complex imaginary(0, 1); - std::array, 4> phaseCoeffArr{1.0, -1. * imaginary, -1.0, - imaginary}; - auto phase_coeff = phaseCoeffArr[_phase]; - auto coeff = rowCoeff; - coeff *= phase_coeff * otherCoeff; - return std::make_pair(coeff, tmp); + return std::make_pair(final_coeff, result); } } // namespace details From b315074dbce9276abed23f654e27e214c984cad7 Mon Sep 17 00:00:00 2001 From: boschmitt <7152025+boschmitt@users.noreply.github.com> Date: Thu, 10 Oct 2024 00:36:59 +0200 Subject: [PATCH 2/5] Simplify spin_op multiplication Signed-off-by: boschmitt <7152025+boschmitt@users.noreply.github.com> --- runtime/cudaq/spin/spin_op.cpp | 134 ++++++++++++++++----------------- 1 file changed, 65 insertions(+), 69 deletions(-) diff --git a/runtime/cudaq/spin/spin_op.cpp b/runtime/cudaq/spin/spin_op.cpp index 3313bde0c6..bee95426f8 100644 --- a/runtime/cudaq/spin/spin_op.cpp +++ b/runtime/cudaq/spin/spin_op.cpp @@ -9,6 +9,7 @@ #include "common/EigenDense.h" #include "common/EigenSparse.h" #include "common/FmtCore.h" +#include #include #include #include @@ -21,7 +22,6 @@ #include #include #include -#include #include #include #include @@ -62,51 +62,51 @@ actionOnBra(spin_op &term, const std::string &bitConfiguration) { return std::make_pair(newConfiguration, coeff); } -std::pair, std::vector> -mult(const std::vector& p1, const std::vector p2, - std::complex &rowCoeff, std::complex &otherCoeff) { - std::size_t numQubits = p1.size() / 2; - std::vector result(p1.size(), false); - int phase = 0; - - int p1Phase = 0; - int p2Phase = 0; - int resultPhase = 0; +std::pair, std::complex> +mult(const std::vector& p1, const std::vector& p2, + const std::complex &p1Coeff, const std::complex &p2Coeff) { + auto [minSize, maxSize] = std::minmax({p1.size(), p2.size()}); + std::size_t minNumQubits = minSize / 2; + std::size_t maxNumQubits = maxSize / 2; + std::vector result(maxSize, false); + int yCount = 0; int cPhase = 0; - // Calculate the result of teh Pauli multiplication and the phase shift - for (std::size_t i = 0; i < numQubits; ++i) { + for (std::size_t i = 0; i < minNumQubits; ++i) { bool p1_x = p1[i]; - bool p1_z = p1[i + numQubits]; + bool p1_z = p1[i + (p1.size() / 2)]; bool p2_x = p2[i]; - bool p2_z = p2[i + numQubits]; + bool p2_z = p2[i + (p2.size() / 2)]; // Compute the resulting Pauli operator result[i] = p1_x ^ p2_x; - result[i + numQubits] = p1_z ^ p2_z; + result[i + maxNumQubits] = p1_z ^ p2_z; - p1Phase += p1_x & p1_z; - p2Phase += p2_x & p2_z; + yCount += (p1_x & p1_z) + (p2_x & p2_z) - (result[i] & result[i + maxNumQubits]); cPhase += p1_x & p2_z; - resultPhase += result[i] & result[i + numQubits]; } - phase = 2*cPhase + p1Phase + p2Phase - resultPhase; + + const std::vector& big = p1.size() < p2.size() ? p2 : p1; + for (std::size_t i = minNumQubits; i < maxNumQubits; ++i) { + result[i] = big[i]; + result[i + maxNumQubits] = big[i + maxNumQubits]; + } // Normalize the phase to a value in the range [0, 3] - phase %= 4; - if (phase < 0) - phase += 4; + int phaseFactor = (2*cPhase + yCount) % 4; + if (phaseFactor < 0) + phaseFactor += 4; // Phase correction factors based on the total phase using namespace std::complex_literals; - std::array, 4> phaseCoeffArr{1.0, -1i, -1.0, 1i}; - std::complex final_coeff = rowCoeff * phaseCoeffArr[phase] * otherCoeff; + std::array, 4> phase{1.0, -1i, -1.0, 1i}; + std::complex resultCoeff = p1Coeff * phase[phaseFactor] * p2Coeff; // Handle the "-0" issue - if (std::abs(final_coeff.real()) < 1e-12) - final_coeff.real(0); + if (std::abs(resultCoeff.real()) < 1e-12) + resultCoeff.real(0); - return std::make_pair(final_coeff, result); + return {result, resultCoeff}; } } // namespace details @@ -428,56 +428,52 @@ spin_op &spin_op::operator-=(const spin_op &v) noexcept { } spin_op &spin_op::operator*=(const spin_op &v) noexcept { - spin_op copy = v; - if (v.num_qubits() > num_qubits()) - expandToNQubits(copy.num_qubits()); - else if (v.num_qubits() < num_qubits()) - copy.expandToNQubits(num_qubits()); - - std::unordered_map, std::complex> newTerms; - std::size_t ourRow = 0, theirRow = 0; - std::vector> composedCoeffs(num_terms() * - copy.num_terms()); - std::vector> composition(num_terms() * copy.num_terms()); - std::map> indexMap; - auto nElements = composition.size(); - for (std::size_t i = 0; i < nElements; i++) { - auto pair = std::make_pair(ourRow, theirRow); - indexMap.emplace(i, pair); - if (theirRow == copy.num_terms() - 1) { - theirRow = 0; - ourRow++; - } else - theirRow++; + using term_and_coeff = std::pair, std::complex>; + std::size_t numTerms = num_terms() * v.num_terms(); + std::vector result(numTerms); + + std::size_t min = std::min(num_terms(), v.num_terms()); + + // Ugh, I take the iterators from the unordered_map to minimize pointer + // chasing. This way we don't need to walk from the begining all the way to + // the element on every iteration. + // + // TODO: Name this things better. + using Iter = std::unordered_map>::const_iterator; + std::vector vec1; + std::vector vec2; + vec1.reserve(terms.size()); + vec2.reserve(v.terms.size()); + for (auto it = terms.begin(); it != terms.end(); ++it) { + vec1.push_back(it); } + for (auto it = v.terms.begin(); it != v.terms.end(); ++it) { + vec2.push_back(it); + } + #if defined(_OPENMP) // Threshold to start OpenMP parallelization. // 16 ~ 4-term * 4-term constexpr std::size_t spin_op_omp_threshold = 16; -#pragma omp parallel for shared(composition) if (nElements > \ - spin_op_omp_threshold) +#pragma omp parallel for shared(result) if (numTerms > spin_op_omp_threshold) #endif - for (std::size_t i = 0; i < nElements; i++) { - auto [j, k] = indexMap[i]; - auto s = terms.begin(); - auto t = copy.terms.begin(); - std::advance(s, j); - std::advance(t, k); - auto res = details::mult(s->first, t->first, s->second, t->second); - composition[i] = res.second; - composedCoeffs[i] = res.first; + for (std::size_t i = 0; i < numTerms; ++i) { + Iter s = vec1[i % min]; + Iter t = vec2[i / min]; + if (terms.size() > v.terms.size()) { + s = vec1[i / min]; + t = vec2[i % min]; + } + result[i] = details::mult(s->first, t->first, s->second, t->second); } - for (std::size_t i = 0; i < nElements; i++) { - auto iter = newTerms.find(composition[i]); - if (iter == newTerms.end()) - newTerms.emplace(composition[i], composedCoeffs[i]); - else - iter->second += composedCoeffs[i]; + terms.clear(); + terms.reserve(numTerms); + for (auto&& [term, coeff] : result) { + auto [it, created] = terms.emplace(term, coeff); + if (!created) + it->second += coeff; } - - terms = newTerms; - return *this; } From a7bcb90771d7756a1c9c803c3446881c7e1b846e Mon Sep 17 00:00:00 2001 From: boschmitt <7152025+boschmitt@users.noreply.github.com> Date: Thu, 17 Oct 2024 17:48:08 +0200 Subject: [PATCH 3/5] Cleanup and add more tests Signed-off-by: boschmitt <7152025+boschmitt@users.noreply.github.com> --- runtime/cudaq/spin/spin_op.cpp | 34 +++---- unittests/spin_op/SpinOpTester.cpp | 151 ++++++++++++++++++----------- 2 files changed, 110 insertions(+), 75 deletions(-) diff --git a/runtime/cudaq/spin/spin_op.cpp b/runtime/cudaq/spin/spin_op.cpp index bee95426f8..d76ce98f5a 100644 --- a/runtime/cudaq/spin/spin_op.cpp +++ b/runtime/cudaq/spin/spin_op.cpp @@ -431,25 +431,19 @@ spin_op &spin_op::operator*=(const spin_op &v) noexcept { using term_and_coeff = std::pair, std::complex>; std::size_t numTerms = num_terms() * v.num_terms(); std::vector result(numTerms); - std::size_t min = std::min(num_terms(), v.num_terms()); - // Ugh, I take the iterators from the unordered_map to minimize pointer - // chasing. This way we don't need to walk from the begining all the way to - // the element on every iteration. - // - // TODO: Name this things better. + // Take the `unordered_map` iterators to minimize pointer chasing when doing + // the cartesian product of the terms of these spin operators. using Iter = std::unordered_map>::const_iterator; - std::vector vec1; - std::vector vec2; - vec1.reserve(terms.size()); - vec2.reserve(v.terms.size()); - for (auto it = terms.begin(); it != terms.end(); ++it) { - vec1.push_back(it); - } - for (auto it = v.terms.begin(); it != v.terms.end(); ++it) { - vec2.push_back(it); - } + std::vector thisTermIt; + std::vector otherTermIt; + thisTermIt.reserve(terms.size()); + otherTermIt.reserve(v.terms.size()); + for (auto it = terms.begin(); it != terms.end(); ++it) + thisTermIt.push_back(it); + for (auto it = v.terms.begin(); it != v.terms.end(); ++it) + otherTermIt.push_back(it); #if defined(_OPENMP) // Threshold to start OpenMP parallelization. @@ -458,11 +452,11 @@ spin_op &spin_op::operator*=(const spin_op &v) noexcept { #pragma omp parallel for shared(result) if (numTerms > spin_op_omp_threshold) #endif for (std::size_t i = 0; i < numTerms; ++i) { - Iter s = vec1[i % min]; - Iter t = vec2[i / min]; + Iter s = thisTermIt[i % min]; + Iter t = otherTermIt[i / min]; if (terms.size() > v.terms.size()) { - s = vec1[i / min]; - t = vec2[i % min]; + s = thisTermIt[i / min]; + t = otherTermIt[i % min]; } result[i] = details::mult(s->first, t->first, s->second, t->second); } diff --git a/unittests/spin_op/SpinOpTester.cpp b/unittests/spin_op/SpinOpTester.cpp index d8176c89b3..945934a6aa 100644 --- a/unittests/spin_op/SpinOpTester.cpp +++ b/unittests/spin_op/SpinOpTester.cpp @@ -10,6 +10,82 @@ #include "cudaq/spin_op.h" +enum Pauli : int8_t { I = 0, X, Y, Z }; +constexpr Pauli paulis[4] = {Pauli::I, Pauli::X, Pauli::Y, Pauli::Z }; + +// Function to multiply two single-qubit Pauli operators +static std::pair, Pauli> multiply_paulis(Pauli a, Pauli b) { + using namespace std::complex_literals; + // I X Y Z + constexpr std::complex table[4][4] = {{ 1., 1., 1, 1}, // I + { 1., 1., 1i, -1i}, // X + { 1., -1i, 1, 1i}, // Y + { 1., 1i, -1i, 1} // Z + }; + if (a == b) return {1.0, Pauli::I}; + if (a == Pauli::I) return {1.0, b}; + if (b == Pauli::I) return {1.0, a}; + return {table[a][b], paulis[a ^ b]}; +} + +// Function to multiply two multi-qubit Pauli words +static std::pair, std::vector> +multiply_pauli_words(const std::vector& a, const std::vector& b, bool verbose = false) { + std::complex phase = 1.0; + std::string info; + std::vector result(a.size(), Pauli::I); + for (size_t i = 0; i < a.size(); ++i) { + auto [p, r] = multiply_paulis(a[i], b[i]); + phase *= p; + result[i] = r; + } + return {phase, result}; +} + +// Generates a pauli word out of a binary representation of it. +static std::vector generate_pauli_word(int64_t id, int64_t num_qubits) { + constexpr int64_t mask = 0x3; + std::vector word(num_qubits, Pauli::I); + for (int64_t i = 0; i < num_qubits; ++i) { + assert((id & mask) < 4); + word[i] = paulis[id & mask]; + id >>= 2; + } + return word; +} + +static std::string generate_pauli_string(const std::vector& word) { + constexpr char paulis_name[4] = {'I', 'X', 'Y', 'Z' }; + std::string result(word.size(), 'I'); + for (int64_t i = 0; i < word.size(); ++i) + result[i] = paulis_name[word[i]]; + return result; +} + +static cudaq::spin_op generate_cudaq_spin(int64_t id, int64_t num_qubits, bool addI = true) { + constexpr int64_t mask = 0x3; + cudaq::spin_op result; + for (int64_t i = 0; i < num_qubits; ++i) { + switch (paulis[id & mask]) { + case Pauli::I: + if (addI) + result *= cudaq::spin::i(i); + break; + case Pauli::X: + result *= cudaq::spin::x(i); + break; + case Pauli::Y: + result *= cudaq::spin::y(i); + break; + case Pauli::Z: + result *= cudaq::spin::z(i); + break; + } + id >>= 2; + } + return result; +} + using namespace cudaq::spin; TEST(SpinOpTester, checkConstruction) { @@ -82,61 +158,26 @@ TEST(SpinOpTester, checkBug178) { } TEST(SpinOpTester, checkMultiplication) { - auto mult = x(0) * y(1); - mult.dump(); - auto mult3 = y(0) * y(1); - mult3.dump(); - - auto tmp = 2 * y(1); - tmp.dump(); - - auto mult2 = x(3) * tmp; - mult2.dump(); - - std::cout << "X * Z: -iY\n"; - (x(3) * z(3)).dump(); - EXPECT_EQ(y(3), x(3) * z(3)); - EXPECT_EQ((x(3) * z(3)).get_coefficient(), std::complex(0, -1)); - - std::cout << "X * X: I\n"; - (x(2) * x(2)).dump(); - EXPECT_EQ(cudaq::spin_op(), x(2) * x(2)); - EXPECT_EQ((x(2) * x(2)).get_coefficient(), std::complex(1, 0)); - - std::cout << "Y * Y: I\n"; - (y(14) * y(14)).dump(); - EXPECT_EQ(cudaq::spin_op(), y(14) * y(14)); - EXPECT_EQ((y(14) * y(14)).get_coefficient(), std::complex(1, 0)); - - std::cout << "Z * Z: I\n"; - (z(0) * z(0)).dump(); - EXPECT_EQ(cudaq::spin_op(), z(0) * z(0)); - EXPECT_EQ((z(0) * z(0)).get_coefficient(), std::complex(1, 0)); - - std::cout << "X * Y: iZ\n"; - (x(3) * y(3)).dump(); - EXPECT_EQ(z(3), x(3) * y(3)); - EXPECT_EQ((x(3) * y(3)).get_coefficient(), std::complex(0, 1)); - - std::cout << "I * I: I\n"; - (i(2) * i(2)).dump(); - EXPECT_EQ(i(2), i(2) * i(2)); - EXPECT_EQ((i(2) * i(2)).get_coefficient(), std::complex(1, 0)); - - std::cout << "I * Z: Z\n"; - (i(3) * i(3)).dump(); - EXPECT_EQ(z(3), i(3) * z(3)); - EXPECT_EQ((i(3) * z(3)).get_coefficient(), std::complex(1, 0)); - - auto tmp2 = 2 * x(0) * x(1) * y(2) * y(3) + 3 * y(0) * y(1) * x(2) * x(3); - std::cout << "START\n"; - tmp2 = tmp2 * tmp2; - tmp2.dump(); - - EXPECT_EQ(2, tmp2.num_terms()); - auto expected = - 13 * i(0) * i(1) * i(2) * i(3) + 12 * z(0) * z(1) * z(2) * z(3); - EXPECT_EQ(expected, tmp2); + for (int num_qubits = 1; num_qubits <= 4; ++num_qubits) { + int64_t num_words = std::pow(4, num_qubits); + for (int64_t i = 0; i < num_words; ++i) { + for (int64_t j = 0; j < num_words; ++j) { + // Expected result: + std::vector a_word = generate_pauli_word(i, num_qubits); + std::vector b_word = generate_pauli_word(j, num_qubits); + auto [phase, result] = multiply_pauli_words(a_word, b_word); + + // Result: + cudaq::spin_op a_spin = generate_cudaq_spin(i, num_qubits); + cudaq::spin_op b_spin = generate_cudaq_spin(j, num_qubits, false); + cudaq::spin_op result_spin = a_spin * b_spin; + + // Check result + EXPECT_EQ(generate_pauli_string(result), result_spin.to_string(false)); + EXPECT_EQ(phase, result_spin.get_coefficient()); + } + } + } } TEST(SpinOpTester, canBuildDeuteron) { From 83c914f309b7a5c3f52ecd8ac3e81da5e7bdde4e Mon Sep 17 00:00:00 2001 From: boschmitt <7152025+boschmitt@users.noreply.github.com> Date: Thu, 17 Oct 2024 17:58:27 +0200 Subject: [PATCH 4/5] Fix formatting Signed-off-by: boschmitt <7152025+boschmitt@users.noreply.github.com> --- runtime/cudaq/spin/spin_op.cpp | 14 ++++--- unittests/spin_op/SpinOpTester.cpp | 61 +++++++++++++++++------------- 2 files changed, 42 insertions(+), 33 deletions(-) diff --git a/runtime/cudaq/spin/spin_op.cpp b/runtime/cudaq/spin/spin_op.cpp index d76ce98f5a..5a0aaa0592 100644 --- a/runtime/cudaq/spin/spin_op.cpp +++ b/runtime/cudaq/spin/spin_op.cpp @@ -63,7 +63,7 @@ actionOnBra(spin_op &term, const std::string &bitConfiguration) { } std::pair, std::complex> -mult(const std::vector& p1, const std::vector& p2, +mult(const std::vector &p1, const std::vector &p2, const std::complex &p1Coeff, const std::complex &p2Coeff) { auto [minSize, maxSize] = std::minmax({p1.size(), p2.size()}); std::size_t minNumQubits = minSize / 2; @@ -82,18 +82,19 @@ mult(const std::vector& p1, const std::vector& p2, result[i] = p1_x ^ p2_x; result[i + maxNumQubits] = p1_z ^ p2_z; - yCount += (p1_x & p1_z) + (p2_x & p2_z) - (result[i] & result[i + maxNumQubits]); + yCount += + (p1_x & p1_z) + (p2_x & p2_z) - (result[i] & result[i + maxNumQubits]); cPhase += p1_x & p2_z; } - const std::vector& big = p1.size() < p2.size() ? p2 : p1; + const std::vector &big = p1.size() < p2.size() ? p2 : p1; for (std::size_t i = minNumQubits; i < maxNumQubits; ++i) { result[i] = big[i]; result[i + maxNumQubits] = big[i + maxNumQubits]; } // Normalize the phase to a value in the range [0, 3] - int phaseFactor = (2*cPhase + yCount) % 4; + int phaseFactor = (2 * cPhase + yCount) % 4; if (phaseFactor < 0) phaseFactor += 4; @@ -435,7 +436,8 @@ spin_op &spin_op::operator*=(const spin_op &v) noexcept { // Take the `unordered_map` iterators to minimize pointer chasing when doing // the cartesian product of the terms of these spin operators. - using Iter = std::unordered_map>::const_iterator; + using Iter = + std::unordered_map>::const_iterator; std::vector thisTermIt; std::vector otherTermIt; thisTermIt.reserve(terms.size()); @@ -463,7 +465,7 @@ spin_op &spin_op::operator*=(const spin_op &v) noexcept { terms.clear(); terms.reserve(numTerms); - for (auto&& [term, coeff] : result) { + for (auto &&[term, coeff] : result) { auto [it, created] = terms.emplace(term, coeff); if (!created) it->second += coeff; diff --git a/unittests/spin_op/SpinOpTester.cpp b/unittests/spin_op/SpinOpTester.cpp index 945934a6aa..51b24e3942 100644 --- a/unittests/spin_op/SpinOpTester.cpp +++ b/unittests/spin_op/SpinOpTester.cpp @@ -11,26 +11,32 @@ #include "cudaq/spin_op.h" enum Pauli : int8_t { I = 0, X, Y, Z }; -constexpr Pauli paulis[4] = {Pauli::I, Pauli::X, Pauli::Y, Pauli::Z }; +constexpr Pauli paulis[4] = {Pauli::I, Pauli::X, Pauli::Y, Pauli::Z}; // Function to multiply two single-qubit Pauli operators -static std::pair, Pauli> multiply_paulis(Pauli a, Pauli b) { +static std::pair, Pauli> multiply_paulis(Pauli a, + Pauli b) { using namespace std::complex_literals; - // I X Y Z - constexpr std::complex table[4][4] = {{ 1., 1., 1, 1}, // I - { 1., 1., 1i, -1i}, // X - { 1., -1i, 1, 1i}, // Y - { 1., 1i, -1i, 1} // Z + // I X Y Z + constexpr std::complex table[4][4] = { + {1., 1., 1, 1}, // I + {1., 1., 1i, -1i}, // X + {1., -1i, 1, 1i}, // Y + {1., 1i, -1i, 1} // Z }; - if (a == b) return {1.0, Pauli::I}; - if (a == Pauli::I) return {1.0, b}; - if (b == Pauli::I) return {1.0, a}; + if (a == b) + return {1.0, Pauli::I}; + if (a == Pauli::I) + return {1.0, b}; + if (b == Pauli::I) + return {1.0, a}; return {table[a][b], paulis[a ^ b]}; } // Function to multiply two multi-qubit Pauli words static std::pair, std::vector> -multiply_pauli_words(const std::vector& a, const std::vector& b, bool verbose = false) { +multiply_pauli_words(const std::vector &a, const std::vector &b, + bool verbose = false) { std::complex phase = 1.0; std::string info; std::vector result(a.size(), Pauli::I); @@ -54,32 +60,33 @@ static std::vector generate_pauli_word(int64_t id, int64_t num_qubits) { return word; } -static std::string generate_pauli_string(const std::vector& word) { - constexpr char paulis_name[4] = {'I', 'X', 'Y', 'Z' }; +static std::string generate_pauli_string(const std::vector &word) { + constexpr char paulis_name[4] = {'I', 'X', 'Y', 'Z'}; std::string result(word.size(), 'I'); for (int64_t i = 0; i < word.size(); ++i) result[i] = paulis_name[word[i]]; return result; } -static cudaq::spin_op generate_cudaq_spin(int64_t id, int64_t num_qubits, bool addI = true) { +static cudaq::spin_op generate_cudaq_spin(int64_t id, int64_t num_qubits, + bool addI = true) { constexpr int64_t mask = 0x3; cudaq::spin_op result; for (int64_t i = 0; i < num_qubits; ++i) { switch (paulis[id & mask]) { - case Pauli::I: - if (addI) - result *= cudaq::spin::i(i); - break; - case Pauli::X: - result *= cudaq::spin::x(i); - break; - case Pauli::Y: - result *= cudaq::spin::y(i); - break; - case Pauli::Z: - result *= cudaq::spin::z(i); - break; + case Pauli::I: + if (addI) + result *= cudaq::spin::i(i); + break; + case Pauli::X: + result *= cudaq::spin::x(i); + break; + case Pauli::Y: + result *= cudaq::spin::y(i); + break; + case Pauli::Z: + result *= cudaq::spin::z(i); + break; } id >>= 2; } From 428bf22b694a648e97897562c7a0a54fc8f1d922 Mon Sep 17 00:00:00 2001 From: boschmitt <7152025+boschmitt@users.noreply.github.com> Date: Thu, 17 Oct 2024 20:37:04 +0200 Subject: [PATCH 5/5] More cleanup Signed-off-by: boschmitt <7152025+boschmitt@users.noreply.github.com> --- runtime/cudaq/spin/spin_op.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/runtime/cudaq/spin/spin_op.cpp b/runtime/cudaq/spin/spin_op.cpp index 5a0aaa0592..a898efd0f1 100644 --- a/runtime/cudaq/spin/spin_op.cpp +++ b/runtime/cudaq/spin/spin_op.cpp @@ -9,7 +9,6 @@ #include "common/EigenDense.h" #include "common/EigenSparse.h" #include "common/FmtCore.h" -#include #include #include #include @@ -434,8 +433,8 @@ spin_op &spin_op::operator*=(const spin_op &v) noexcept { std::vector result(numTerms); std::size_t min = std::min(num_terms(), v.num_terms()); - // Take the `unordered_map` iterators to minimize pointer chasing when doing - // the cartesian product of the terms of these spin operators. + // Put the `unordered_map` iterators into vectors to minimize pointer chasing + // when doing the cartesian product of the spin operators' terms. using Iter = std::unordered_map>::const_iterator; std::vector thisTermIt;