From f7f0f21de420b747da86d06670169a5a30738a41 Mon Sep 17 00:00:00 2001 From: Jun Doi Date: Tue, 20 Feb 2024 13:10:40 +0900 Subject: [PATCH 1/5] Add support > 63 qubits for extended stabilizer --- src/framework/bitvector.hpp | 374 +++++++++ .../extended_stabilizer/ch_runner.hpp | 127 +-- .../chlib/chstabilizer.hpp | 729 ++++++++---------- .../extended_stabilizer/chlib/core.hpp | 127 +-- .../extended_stabilizer_state.hpp | 39 +- src/simulators/extended_stabilizer/gates.hpp | 4 + 6 files changed, 874 insertions(+), 526 deletions(-) create mode 100644 src/framework/bitvector.hpp diff --git a/src/framework/bitvector.hpp b/src/framework/bitvector.hpp new file mode 100644 index 0000000000..9851346106 --- /dev/null +++ b/src/framework/bitvector.hpp @@ -0,0 +1,374 @@ +/** + * This code is part of Qiskit. + * + * (C) Copyright IBM 2018, 2019, 2024. + * + * This code is licensed under the Apache License, Version 2.0. You may + * obtain a copy of this license in the LICENSE.txt file in the root directory + * of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. + * + * Any modifications or derivative works of this code must retain this + * copyright notice, and modified files need to carry a notice indicating + * that they have been altered from the originals. + */ + +#ifndef _aer_framework_bitvector_hpp_ +#define _aer_framework_bitvector_hpp_ + +#include "framework/rng.hpp" +#include "framework/types.hpp" +#include "framework/utils.hpp" + +namespace AER { + +//============================================================================ +// Bit vestor class +//============================================================================ + +class BitVector { +protected: + reg_t bits_; + uint_t num_bits_; + const static size_t REG_SIZE = 64; + const static size_t REG_BITS = 6; + const static size_t REG_MASK = (1ull << REG_BITS) - 1; + +public: + BitVector() { num_bits_ = 0; } + BitVector(uint_t nbits, bool init = false) { allocate(nbits, init); } + BitVector(const BitVector &src) { + bits_ = src.bits_; + num_bits_ = src.num_bits_; + } + + uint_t num_bits() { return num_bits_; } + uint_t length() { return bits_.size(); } + + void allocate(uint_t n, bool init = false) { + uint_t size = n >> REG_BITS; + if (size == 0) + size = 1; + if (init) + bits_.resize(size, 0ull); + else + bits_.resize(size); + num_bits_ = n; + } + + BitVector &operator=(const BitVector &src) { + bits_ = src.bits_; + num_bits_ = src.num_bits_; + return *this; + } + BitVector &operator=(const std::string &src) { + from_string(src); + return *this; + } + BitVector &operator=(const reg_t &src) { + from_vector(src); + return *this; + } + + // copy with swap + void map(const BitVector &src, const reg_t map); + + // bit access + inline bool get(const uint_t idx) const { + uint_t pos = idx >> REG_BITS; + uint_t bit = idx & REG_MASK; + return (((bits_[pos] >> bit) & 1ull) == 1ull); + } + inline bool operator[](const uint_t idx) const { return get(idx); } + inline uint_t &operator()(const uint_t pos) { return bits_[pos]; } + inline uint_t operator()(const uint_t pos) const { return bits_[pos]; } + + void set(const uint_t idx, const bool val) { + uint_t pos = idx >> REG_BITS; + uint_t bit = idx & REG_MASK; + uint_t mask = ~(1ull << bit); + bits_[pos] &= mask; + bits_[pos] |= (((uint_t)val) << bit); + } + + void zero(void); + void rand(RngEngine &rng); + + bool operator==(const BitVector &src) const; + bool operator!=(const BitVector &src) const; + bool is_zero(void); + + // logical operators + BitVector &operator|=(const BitVector &src); + BitVector &operator&=(const BitVector &src); + BitVector &operator^=(const BitVector &src); + + friend BitVector operator|(const BitVector &v0, const BitVector &v1); + friend BitVector operator&(const BitVector &v0, const BitVector &v1); + friend BitVector operator^(const BitVector &v0, const BitVector &v1); + friend BitVector operator~(const BitVector &v); + + void set_or(const uint_t idx, const bool b) { + uint_t pos = idx >> REG_BITS; + uint_t bit = idx & REG_MASK; + uint_t val = (uint_t)b << bit; + bits_[pos] |= val; + } + + void set_and(const uint_t idx, const bool b) { + uint_t pos = idx >> REG_BITS; + uint_t bit = idx & REG_MASK; + uint_t val = (~(1ull << bit)) | ((uint_t)b << bit); + bits_[pos] &= val; + } + + void set_xor(const uint_t idx, const bool b) { + uint_t pos = idx >> REG_BITS; + uint_t bit = idx & REG_MASK; + uint_t val = (uint_t)b << bit; + bits_[pos] ^= val; + } + + void set_not(const uint_t idx) { + uint_t pos = idx >> REG_BITS; + uint_t bit = idx & REG_MASK; + uint_t val = 1ull << bit; + bits_[pos] ^= val; + } + + // convert from other data + void from_uint(const uint_t nbits, const uint_t src); + void from_string(const std::string &src); + void from_vector(const reg_t &src); + void from_vector_with_map(const reg_t &src, const reg_t &map); + + // convert to other data types + std::string to_string(); + std::string to_hex_string(bool prefix = true); + reg_t to_vector(); + + bool parity(void); + uint_t popcount(void); +}; + +void BitVector::map(const BitVector &src, const reg_t map) { + allocate(map.size()); + + for (uint_t i = 0; i < map.size(); i++) { + set(i, src[map[i]]); + } +} + +void BitVector::from_uint(const uint_t nbits, const uint_t src) { + allocate(nbits); + bits_[0] = src; +} + +void BitVector::from_string(const std::string &src) { + allocate(src.size()); + + uint_t pos = 0; + for (uint_t i = 0; i < bits_.size(); i++) { + uint_t n = REG_SIZE; + uint_t val = 0; + if (n > num_bits_ - pos) + n = num_bits_ - pos; + for (uint_t j = 0; j < n; j++) { + val |= (((uint_t)(src[num_bits_ - 1 - pos] == '1')) << j); + pos++; + } + bits_[i] = val; + } +} + +void BitVector::from_vector(const reg_t &src) { + allocate(src.size()); + + uint_t pos = 0; + for (uint_t i = 0; i < bits_.size(); i++) { + uint_t n = REG_SIZE; + uint_t val = 0; + if (n > num_bits_ - pos) + n = num_bits_ - pos; + for (uint_t j = 0; j < n; j++) { + val |= ((src[pos++] & 1ull) << j); + } + bits_[i] = val; + } +} + +void BitVector::from_vector_with_map(const reg_t &src, const reg_t &map) { + allocate(src.size()); + + uint_t pos = 0; + for (uint_t i = 0; i < bits_.size(); i++) { + uint_t n = REG_SIZE; + uint_t val = 0; + if (n > num_bits_ - pos) + n = num_bits_ - pos; + for (uint_t j = 0; j < n; j++) { + val |= ((src[map[pos++]] & 1ull) << j); + } + bits_[i] = val; + } +} + +std::string BitVector::to_string(void) { + std::string str; + for (uint_t i = 0; i < num_bits_; i++) { + if (get(num_bits_ - 1 - i)) + str += '1'; + else + str += '0'; + } + return str; +} + +std::string BitVector::to_hex_string(bool prefix) { + // initialize output string + std::string hex = (prefix) ? "0x" : ""; + + for (uint_t i = 0; i < bits_.size(); i++) { + if (i == 0) { + uint_t n = num_bits_ & (REG_SIZE - 1); + uint_t val = bits_[bits_.size() - 1] & ((1ull << n) - 1); + + std::stringstream ss; + ss << std::hex << val; + hex += ss.str(); + } else { + std::stringstream ss; + ss << std::hex << bits_[bits_.size() - 1 - i]; + std::string part = ss.str(); + part.insert(0, (REG_SIZE / 4) - part.size(), '0'); + hex += part; + } + } + return hex; +} + +reg_t BitVector::to_vector(void) { + reg_t ret(num_bits_); + for (uint_t i = 0; i < num_bits_; i++) { + ret[i] = (uint_t)get(i); + } + return ret; +} + +void BitVector::zero(void) { + for (uint_t i = 0; i < bits_.size(); i++) { + bits_[i] = 0; + } +} + +void BitVector::rand(RngEngine &rng) { + if (bits_.size() == 1) { + bits_[0] = rng.rand_int(0ull, (1ull << num_bits_) - 1ull); + } else { + double r = rng.rand(); + double c = 0.5; + + for (int_t i = num_bits_ - 1; i >= 0; i--) { + if (r >= c) { + set(i, true); + r -= c; + } else + set(i, false); + c /= 2.0; + } + } +} + +bool BitVector::operator==(const BitVector &src) const { + if (num_bits_ != src.num_bits_) + return false; + for (uint_t i = 0; i < bits_.size(); i++) { + if (bits_[i] != src.bits_[i]) + return false; + } + return true; +} + +bool BitVector::operator!=(const BitVector &src) const { + if (num_bits_ != src.num_bits_) + return true; + for (uint_t i = 0; i < bits_.size(); i++) { + if (bits_[i] != src.bits_[i]) + return true; + } + return false; +} + +bool BitVector::is_zero(void) { + if (num_bits_ == 0) + return true; + for (uint_t i = 0; i < bits_.size(); i++) { + if (bits_[i] != 0) + return false; + } + return true; +} + +BitVector &BitVector::operator|=(const BitVector &src) { + uint_t size = std::min(src.bits_.size(), bits_.size()); + for (uint_t i = 0; i < size; i++) + bits_[i] |= src.bits_[i]; + return *this; +} + +BitVector &BitVector::operator&=(const BitVector &src) { + uint_t size = std::min(src.bits_.size(), bits_.size()); + for (uint_t i = 0; i < size; i++) + bits_[i] &= src.bits_[i]; + return *this; +} + +BitVector &BitVector::operator^=(const BitVector &src) { + uint_t size = std::min(src.bits_.size(), bits_.size()); + for (uint_t i = 0; i < size; i++) + bits_[i] ^= src.bits_[i]; + return *this; +} + +BitVector operator|(const BitVector &v0, const BitVector &v1) { + BitVector ret(v0); + return ret |= v1; +} + +BitVector operator&(const BitVector &v0, const BitVector &v1) { + BitVector ret(v0); + return ret &= v1; +} + +BitVector operator^(const BitVector &v0, const BitVector &v1) { + BitVector ret(v0); + return ret ^= v1; +} + +BitVector operator~(const BitVector &v) { + BitVector ret(v); + for (uint_t i = 0; i < ret.bits_.size(); i++) { + ret.bits_[i] = ~ret.bits_[i]; + } + return ret; +} + +bool BitVector::parity(void) { + uint_t ret = 0; + for (uint_t i = 0; i < bits_.size(); i++) { + ret ^= AER::Utils::hamming_parity(bits_[i]); + } + return ret != 0; +} + +uint_t BitVector::popcount(void) { + uint_t ret = 0; + for (uint_t i = 0; i < bits_.size(); i++) { + ret += AER::Utils::popcount(bits_[i]); + } + return ret; +} + +//------------------------------------------------------------------------------ +} // end namespace AER +//------------------------------------------------------------------------------ +#endif // _aer_framework_bitvector_hpp_ diff --git a/src/simulators/extended_stabilizer/ch_runner.hpp b/src/simulators/extended_stabilizer/ch_runner.hpp index 489d6b77ad..1887a30825 100644 --- a/src/simulators/extended_stabilizer/ch_runner.hpp +++ b/src/simulators/extended_stabilizer/ch_runner.hpp @@ -37,6 +37,7 @@ #include #endif +namespace AER { namespace CHSimulator { using chstabilizer_t = StabilizerState; @@ -62,7 +63,7 @@ class Runner { bool accept_; complex_t old_ampsum_; - uint_t x_string_; + BitVector x_string_; uint_t last_proposal_; void init_metropolis(AER::RngEngine &rng); @@ -135,23 +136,24 @@ class Runner { std::vector generators, AER::RngEngine &rng); // Metropolis Estimation for sampling from the output distribution - uint_t metropolis_estimation(uint_t n_steps, AER::RngEngine &rng); - std::vector metropolis_estimation(uint_t n_steps, uint_t n_shots, - AER::RngEngine &rng); + BitVector metropolis_estimation(uint_t n_steps, AER::RngEngine &rng); + std::vector metropolis_estimation(uint_t n_steps, uint_t n_shots, + AER::RngEngine &rng); // NE Sampling Routines - uint_t ne_single_sample(uint_t default_samples, uint_t repetitions, - bool preserve_states, const AER::reg_t &qubits, - AER::RngEngine &rng); + BitVector ne_single_sample(uint_t default_samples, uint_t repetitions, + bool preserve_states, const AER::reg_t &qubits, + AER::RngEngine &rng); std::vector ne_probabilities(uint_t default_samples, uint_t repetitions, const AER::reg_t &qubits, AER::RngEngine &rng); // Efficient Sampler for the output distribution of a stabilizer state - uint_t stabilizer_sampler(AER::RngEngine &rng); - std::vector stabilizer_sampler(uint_t n_shots, AER::RngEngine &rng); + BitVector stabilizer_sampler(AER::RngEngine &rng); + std::vector stabilizer_sampler(uint_t n_shots, + AER::RngEngine &rng); // Utilities for the state-vector snapshot. - complex_t amplitude(uint_t x_measure); + complex_t amplitude(BitVector &x_measure); AER::Vector statevector(); }; @@ -169,6 +171,7 @@ void Runner::initialize(uint_t num_qubits) { num_threads_ = 1; states_ = std::vector(1, chstabilizer_t(num_qubits)); coefficients_.push_back(complex_t(1., 0.)); + x_string_.allocate(n_qubits_, true); } void Runner::initialize_decomposition(uint_t n_states, double delta) { @@ -407,10 +410,11 @@ double Runner::norm_estimation(uint_t n_samples, uint_t repetitions, const int_t NQUBITS = n_qubits_; std::vector xi_samples(repetitions, 0.); for (uint_t m = 0; m < repetitions; m++) { - std::vector adiag_1(n_samples, 0ULL); - std::vector adiag_2(n_samples, 0ULL); - std::vector> a(n_samples, - std::vector(n_qubits_, 0ULL)); + std::vector adiag_1(n_samples, BitVector(n_qubits_, true)); + std::vector adiag_2(n_samples, BitVector(n_qubits_, true)); + std::vector> a( + n_samples, + std::vector(n_qubits_, BitVector(n_qubits_, true))); #pragma omp parallel if (num_threads_ > 1) num_threads(num_threads_) { #ifdef _WIN32 @@ -422,13 +426,13 @@ double Runner::norm_estimation(uint_t n_samples, uint_t repetitions, for (int_t i = 0; i < NQUBITS; i++) { for (int_t j = i; j < NQUBITS; j++) { if (rng.rand() < 0.5) { - a[l][i] |= (1ULL << j); - a[l][j] |= (1ULL << i); + a[l][i].set_or(j, true); + a[l][j].set_or(i, true); } } - adiag_1[l] |= (a[l][i] & (1ULL << i)); + adiag_1[l].set_or(i, a[l][i][i]); if (rng.rand() < 0.5) { - adiag_2[l] |= (1ULL << i); + adiag_2[l].set_or(i, true); } } } @@ -452,9 +456,10 @@ double Runner::norm_estimation(uint_t n_samples, uint_t repetitions, return norm_estimation(n_samples, repetitions, rng); } -uint_t Runner::ne_single_sample(uint_t default_samples, uint_t repetitions, - bool preserve_states, const AER::reg_t &qubits, - AER::RngEngine &rng) { +BitVector Runner::ne_single_sample(uint_t default_samples, uint_t repetitions, + bool preserve_states, + const AER::reg_t &qubits, + AER::RngEngine &rng) { uint_t n_samples = std::llrint(4 * std::pow(qubits.size(), 2)); if (default_samples > n_samples) { n_samples = default_samples; @@ -462,10 +467,10 @@ uint_t Runner::ne_single_sample(uint_t default_samples, uint_t repetitions, double denominator = norm_estimation(n_samples, repetitions, rng); std::vector generators; std::vector states_cache(states_); - uint_t out_string = ZERO; + BitVector out_string(n_qubits_, true); for (uint_t i = 0; i < qubits.size(); i++) { - pauli_t generator; - generator.Z = (1ULL << qubits[i]); + pauli_t generator(n_qubits_); + generator.Z.set(qubits[i], true); apply_pauli(generator); // Compute probability this bit is 0, given the previous assignemnts double numerator = norm_estimation(n_samples, repetitions, rng); @@ -480,7 +485,7 @@ uint_t Runner::ne_single_sample(uint_t default_samples, uint_t repetitions, generators.push_back(generator); states_ = states_cache; apply_pauli_projector(generators); - out_string ^= (1ULL << qubits[i]); + out_string.set_xor(qubits[i], one); denominator = (1 - p_zero_given_observations) * denominator; } } @@ -499,8 +504,8 @@ std::vector Runner::ne_probabilities(uint_t default_samples, std::vector states_cache(states_); std::vector generators; for (uint_t i = 0; i < qubits.size(); i++) { - pauli_t generator; - generator.Z = 1ULL << qubits[i]; + pauli_t generator(n_qubits_); + generator.Z.set(qubits[i], true); generators.push_back(generator); } double norm = norm_estimation(default_samples, repetitions, rng); @@ -521,7 +526,7 @@ std::vector Runner::ne_probabilities(uint_t default_samples, return probs; } -uint_t Runner::metropolis_estimation(uint_t n_steps, AER::RngEngine &rng) { +BitVector Runner::metropolis_estimation(uint_t n_steps, AER::RngEngine &rng) { init_metropolis(rng); for (uint_t i = 0; i < n_steps; i++) { metropolis_step(rng); @@ -529,10 +534,10 @@ uint_t Runner::metropolis_estimation(uint_t n_steps, AER::RngEngine &rng) { return x_string_; } -std::vector Runner::metropolis_estimation(uint_t n_steps, - uint_t n_shots, - AER::RngEngine &rng) { - std::vector shots(n_shots, zer); +std::vector Runner::metropolis_estimation(uint_t n_steps, + uint_t n_shots, + AER::RngEngine &rng) { + std::vector shots(n_shots, zer); shots[0] = metropolis_estimation(n_steps, rng); for (uint_t i = 1; i < n_shots; i++) { metropolis_step(rng); @@ -544,8 +549,7 @@ std::vector Runner::metropolis_estimation(uint_t n_steps, void Runner::init_metropolis(AER::RngEngine &rng) { accept_ = false; // Random initial x_string from RngEngine - uint_t max = (1ULL << n_qubits_) - 1; - x_string_ = rng.rand_int(ZERO, max); + x_string_.rand(rng); // = rng.rand_int(ZERO, max); last_proposal_ = 0; double local_real = 0., local_imag = 0.; const int_t END = num_states_; @@ -564,7 +568,7 @@ void Runner::init_metropolis(AER::RngEngine &rng) { void Runner::metropolis_step(AER::RngEngine &rng) { uint_t proposal = rng.rand(0ULL, n_qubits_); if (accept_) { - x_string_ ^= (one << last_proposal_); + x_string_.set_xor(last_proposal_, one); } double real_part = 0., imag_part = 0.; if (accept_ == 0) { @@ -610,19 +614,20 @@ void Runner::metropolis_step(AER::RngEngine &rng) { } } -uint_t Runner::stabilizer_sampler(AER::RngEngine &rng) { - uint_t max = (1ULL << n_qubits_) - 1; - return states_[0].Sample(rng.rand_int(ZERO, max)); +BitVector Runner::stabilizer_sampler(AER::RngEngine &rng) { + BitVector t(n_qubits_); + t.rand(rng); + return states_[0].Sample(t); } -std::vector Runner::stabilizer_sampler(uint_t n_shots, - AER::RngEngine &rng) { +std::vector Runner::stabilizer_sampler(uint_t n_shots, + AER::RngEngine &rng) { if (num_states_ > 1) { throw std::invalid_argument( "CH::Runner::stabilizer_sampler: This method can only be used for a " "single Stabilizer state.\n"); } - std::vector shots; + std::vector shots; shots.reserve(n_shots); for (uint_t i = 0; i < n_shots; i++) { shots.push_back(stabilizer_sampler(rng)); @@ -630,7 +635,7 @@ std::vector Runner::stabilizer_sampler(uint_t n_shots, return shots; } -complex_t Runner::amplitude(uint_t x_measure) { +complex_t Runner::amplitude(BitVector &x_measure) { double real_part = 0., imag_part = 0.; // Splitting the reduction guarantees support on more OMP versions. const int_t END = num_states_; @@ -645,13 +650,18 @@ complex_t Runner::amplitude(uint_t x_measure) { } AER::Vector Runner::statevector() { + if (n_qubits_ > 63) { + throw std::runtime_error("savestatevector is not supported > 63 qubits"); + } uint_t ceil = 1ULL << n_qubits_; AER::Vector svector(ceil, false); double norm_squared = 0; for (uint_t i = 0; i < ceil; i++) { - svector[i] = amplitude(i); + BitVector t(n_qubits_); + t(0) = i; + svector[i] = amplitude(t); norm_squared += svector[i].real() * svector[i].real() + svector[i].imag() * svector[i].imag(); } @@ -687,26 +697,39 @@ std::vector Runner::serialize_decomposition() const { json_t Runner::serialize_state(uint_t rank) const { json_t js = json_t::object(); std::vector gamma; - std::vector M; - std::vector F; - std::vector G; + std::vector M; + std::vector F; + std::vector G; gamma.reserve(n_qubits_); M = states_[rank].MMatrix(); F = states_[rank].FMatrix(); G = states_[rank].GMatrix(); - uint_t gamma1 = states_[rank].Gamma1(); - uint_t gamma2 = states_[rank].Gamma2(); + BitVector gamma1 = states_[rank].Gamma1(); + BitVector gamma2 = states_[rank].Gamma2(); for (uint_t i = 0; i < n_qubits_; i++) { - gamma.push_back(((gamma1 >> i) & 1ULL) + 2 * ((gamma2 >> i) & 1ULL)); + gamma.push_back((unsigned)gamma1[i] + 2 * (unsigned)gamma2[i]); } js["gamma"] = gamma; - js["M"] = M; - js["F"] = F; - js["G"] = G; + std::vector s(n_qubits_); + for (uint_t i = 0; i < n_qubits_; i++) { + s[i] = M[i].to_string(); + } + js["M"] = s; + for (uint_t i = 0; i < n_qubits_; i++) { + s[i] = F[i].to_string(); + } + js["F"] = s; + for (uint_t i = 0; i < n_qubits_; i++) { + s[i] = G[i].to_string(); + } + js["G"] = s; js["internal_cofficient"] = states_[rank].Omega().to_complex(); js["coefficient"] = coefficients_[rank]; return js; } } // namespace CHSimulator +//------------------------------------------------------------------------------ +} // end namespace AER +//------------------------------------------------------------------------------ #endif diff --git a/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp b/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp index 1d4e27fd39..adb4ecbb0f 100644 --- a/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp +++ b/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp @@ -26,6 +26,7 @@ #include "core.hpp" #include "framework/utils.hpp" +namespace AER { namespace CHSimulator { // Clifford simulator based on the CH-form class StabilizerState { @@ -37,13 +38,13 @@ class StabilizerState { // n = number of qubits // qubits are numbered as q=0,1,...,n-1 - uint_fast64_t NQubits() const { return n; } - scalar_t Omega() const { return omega; } - uint_fast64_t Gamma1() const { return gamma1; } - uint_fast64_t Gamma2() const { return gamma2; } - std::vector GMatrix() const { return G; } - std::vector FMatrix() const { return F; } - std::vector MMatrix() const { return M; } + uint_fast64_t NQubits() const { return n_qubits_; } + scalar_t Omega() const { return omega_; } + BitVector Gamma1() const { return gamma1_; } + BitVector Gamma2() const { return gamma2_; } + std::vector GMatrix() const { return G_; } + std::vector FMatrix() const { return F_; } + std::vector MMatrix() const { return M_; } // Clifford gates void CX(unsigned q, unsigned r); // q=control, r=target @@ -56,30 +57,29 @@ class StabilizerState { void Sdag(unsigned q); // S inverse // state preps - void CompBasisVector(uint_fast64_t x); // prepares |x> - void HadamardBasisVector(uint_fast64_t x); // prepares H(x)|00...0> + void CompBasisVector(BitVector &x); // prepares |x> + void HadamardBasisVector(BitVector &x); // prepares H(x)|00...0> // measurements - scalar_t Amplitude(uint_fast64_t x); // computes the amplitude - uint_fast64_t Sample(); // returns a sample from the distribution ||^2 - uint_fast64_t Sample(uint_fast64_t v_mask); + scalar_t Amplitude(BitVector &x); // computes the amplitude + BitVector Sample(); // returns a sample from the distribution ||^2 + BitVector Sample(BitVector &v_mask); void MeasurePauli(const pauli_t P); // applies a gate (I+P)/2 // where P is an arbitrary Pauli operator void MeasurePauliProjector(const std::vector &generators); - inline scalar_t ScalarPart() { return omega; } + inline scalar_t ScalarPart() { return omega_; } // InnerProduct & Norm Estimation - scalar_t InnerProduct(const uint_fast64_t &A_diag1, - const uint_fast64_t &A_diag2, - const std::vector &A); + scalar_t InnerProduct(const BitVector &A_diag1, const BitVector &A_diag2, + const std::vector &A); // Metropolis updates: // To initialize Metropolis by a state x call Amplitude(x) scalar_t ProposeFlip(unsigned flip_pos); // returns the amplitude // where x'=bitflip(x,q) // x = current Metropolis state - inline void AcceptFlip() { P = Q; } // accept the proposed bit flip + inline void AcceptFlip() { P_ = Q_; } // accept the proposed bit flip friend double NormEstimate(std::vector &states, @@ -98,21 +98,21 @@ class StabilizerState { #endif private: - unsigned n; + unsigned n_qubits_; // define CH-data (F,G,M,gamma,v,s,omega) see Section IIIA for notations // stabilizer tableaux of the C-layer - uint_fast64_t gamma1; // phase vector gamma = gamma1 + 2*gamma2 (mod 4) - uint_fast64_t gamma2; + BitVector gamma1_; // phase vector gamma = gamma1 + 2*gamma2 (mod 4) + BitVector gamma2_; // each column of F,G,M is represented by uint_fast64_t integer - std::vector F; // F[i] = i-th column of F - std::vector G; // G[i] = i-th column of G - std::vector M; // M[i] = i-th column of M + std::vector F_; // F[i] = i-th column of F + std::vector G_; // G[i] = i-th column of G + std::vector M_; // M[i] = i-th column of M - uint_fast64_t v; // H-layer U_H=H(v) - uint_fast64_t s; // initial state |s> - scalar_t omega; // scalar factor such that |phi>=omega*U_C*U_H|s> + BitVector v_; // H-layer U_H=H(v) + BitVector s_; // initial state |s> + scalar_t omega_; // scalar factor such that |phi>=omega*U_C*U_H|s> // multiplies the C-layer on the right by a C-type gate void RightCX(unsigned q, unsigned r); // q=control, r=target @@ -122,25 +122,16 @@ class StabilizerState { void RightSdag(unsigned q); // computes a Pauli operator U_C^{-1}X(x)U_C - pauli_t GetPauliX(uint_fast64_t x); + pauli_t GetPauliX(BitVector &x); // replace the initial state |s> in the CH-form by a superposition // (|t> + i^b |u>)*sqrt(1/2) as described in Proposition 3 // update the CH-form - void UpdateSvector(uint_fast64_t t, uint_fast64_t u, unsigned b); - - // F-transposed and M-transposed. Need this to compute amplitudes - std::vector FT; // FT[i] = i-th row of F - std::vector MT; // MT[i] = i-th row of M - // compute the transposed matrices - void TransposeF(); - void TransposeM(); - bool isReadyFT; // true if F-transposed is available - bool isReadyMT; // true if M-transposed is available + void UpdateSvector(BitVector t, BitVector u, unsigned b); // auxiliary Pauli operators - pauli_t P; - pauli_t Q; + pauli_t P_; + pauli_t Q_; }; using cdouble = std::complex; @@ -152,241 +143,211 @@ using cdouble = std::complex; // Clifford simulator based on the CH-form for for n<=64 qubits StabilizerState::StabilizerState(const unsigned n_qubits) - : n(n_qubits), // number of qubits - gamma1(zer), gamma2(zer), F(n, zer), G(n, zer), M(n, zer), v(zer), s(zer), - omega(), FT(n, zer), MT(n, zer), isReadyFT(false), isReadyMT(false) { + : n_qubits_(n_qubits), // number of qubits + gamma1_(n_qubits_, true), gamma2_(n_qubits_, true), + F_(n_qubits_, BitVector(n_qubits_, true)), + G_(n_qubits_, BitVector(n_qubits_, true)), + M_(n_qubits_, BitVector(n_qubits_, true)), v_(n_qubits_, true), + s_(n_qubits_, true), omega_(), P_(n_qubits_), Q_(n_qubits_) { // initialize by the basis vector 00...0 - if (n > 63) { - throw std::invalid_argument( - "The CH simulator only supports up to 63 qubits.\n"); - } // G and F are identity matrices - for (unsigned q = 0; q < n; q++) { - G[q] = (one << q); - F[q] = (one << q); + for (unsigned q = 0; q < n_qubits_; q++) { + G_[q].set(q, one); + F_[q].set(q, one); } - omega.makeOne(); + omega_.makeOne(); } StabilizerState::StabilizerState(const StabilizerState &rhs) - : n(rhs.n), // number of qubits - gamma1(rhs.gamma1), gamma2(rhs.gamma2), F(rhs.F), G(rhs.G), M(rhs.M), - v(rhs.v), s(rhs.s), omega(rhs.omega), FT(rhs.FT), MT(rhs.MT), - isReadyFT(rhs.isReadyFT), isReadyMT(rhs.isReadyMT) {} - -void StabilizerState::CompBasisVector(uint_fast64_t x) { - s = x; - v = zer; - gamma1 = zer; - gamma2 = zer; - omega.makeOne(); + : n_qubits_(rhs.n_qubits_), // number of qubits + gamma1_(rhs.gamma1_), gamma2_(rhs.gamma2_), F_(rhs.F_), G_(rhs.G_), + M_(rhs.M_), v_(rhs.v_), s_(rhs.s_), omega_(rhs.omega_), P_(rhs.P_), + Q_(rhs.Q_) {} + +void StabilizerState::CompBasisVector(BitVector &x) { + s_ = x; + v_.zero(); + gamma1_.zero(); + gamma2_.zero(); + omega_.makeOne(); // G and F are identity matrices, M is zero matrix - for (unsigned q = 0; q < n; q++) { - M[q] = zer; - G[q] = (one << q); - F[q] = (one << q); + for (unsigned q = 0; q < n_qubits_; q++) { + M_[q].zero(); + G_[q].set(q, one); + F_[q].set(q, one); } - isReadyFT = false; - isReadyMT = false; } -void StabilizerState::HadamardBasisVector(uint_fast64_t x) { - s = zer; - v = x; - gamma1 = zer; - gamma2 = zer; - omega.makeOne(); +void StabilizerState::HadamardBasisVector(BitVector &x) { + s_.zero(); + v_ = x; + gamma1_.zero(); + gamma2_.zero(); + omega_.makeOne(); // G and F are identity matrices, M is zero matrix - for (unsigned q = 0; q < n; q++) { - M[q] = zer; - G[q] = (one << q); - F[q] = (one << q); + for (unsigned q = 0; q < n_qubits_; q++) { + M_[q].zero(); + G_[q].set(q, one); + F_[q].set(q, one); } - isReadyFT = false; - isReadyMT = false; } void StabilizerState::RightS(unsigned q) { - isReadyMT = false; // we are going to change M - M[q] ^= F[q]; + M_[q] ^= F_[q]; // update phase vector: gamma[p] gets gamma[p] - F_{p,q} (mod 4) for all p - gamma2 ^= F[q] ^ (gamma1 & F[q]); - gamma1 ^= F[q]; + gamma2_ ^= F_[q] ^ (gamma1_ & F_[q]); + gamma1_ ^= F_[q]; } void StabilizerState::RightSdag(unsigned q) { - isReadyMT = false; // we are going to change M - M[q] ^= F[q]; + M_[q] ^= F_[q]; // update phase vector: gamma[p] gets gamma[p] + F_{p,q} (mod 4) for all p - gamma1 ^= F[q]; - gamma2 ^= F[q] ^ (gamma1 & F[q]); + gamma1_ ^= F_[q]; + gamma2_ ^= F_[q] ^ (gamma1_ & F_[q]); } void StabilizerState::RightZ(unsigned q) { // update phase vector: gamma[p] gets gamma[p] + 2F_{p,q} (mod 4) for all p - gamma2 ^= F[q]; + gamma2_ ^= F_[q]; } void StabilizerState::S(unsigned q) { - isReadyMT = false; // we are going to change M - uint_fast64_t C = (one << q); - for (unsigned p = 0; p < n; p++) { - M[p] ^= ((G[p] >> q) & one) * C; + for (unsigned p = 0; p < n_qubits_; p++) { + M_[p].set_xor(q, G_[p][q]); } // update phase vector: gamma[q] gets gamma[q] - 1 - gamma1 ^= C; - gamma2 ^= ((gamma1 >> q) & one) * C; + gamma1_.set_xor(q, one); + gamma2_.set_xor(q, gamma1_[q]); } void StabilizerState::Sdag(unsigned q) { - isReadyMT = false; // we are going to change M - uint_fast64_t C = (one << q); - for (unsigned p = 0; p < n; p++) { - M[p] ^= ((G[p] >> q) & one) * C; + for (unsigned p = 0; p < n_qubits_; p++) { + M_[p].set_xor(q, G_[p][q]); } // update phase vector: gamma[q] gets gamma[q] + 1 - gamma2 ^= ((gamma1 >> q) & one) * C; - gamma1 ^= C; + gamma2_.set_xor(q, gamma1_[q]); + gamma1_.set_xor(q, one); } void StabilizerState::Z(unsigned q) { // update phase vector: gamma[q] gets gamma[q] + 2 - gamma2 ^= (one << q); + gamma2_.set_xor(q, one); } void StabilizerState::X(unsigned q) { // Commute the Pauli through UC - if (!isReadyMT) { - TransposeM(); - } - if (!isReadyFT) { - TransposeF(); + BitVector x_string(n_qubits_); + BitVector z_string(n_qubits_); + + for (uint_t i = 0; i < n_qubits_; i++) { + x_string.set(i, F_[i][q]); + z_string.set(i, M_[i][q]); } - uint_fast64_t x_string = FT[q]; - uint_fast64_t z_string = MT[q]; + // Initial phase correction - int phase = 2 * ((gamma1 >> q) & one) + 4 * ((gamma2 >> q) & one); + uint_t phase = 2 * (uint_t)gamma1_[q] + 4 * (uint_t)gamma2_[q]; // Commute the z_string through the hadamard layer // Each z that hits a hadamard becomes a Pauli X - s ^= (z_string & v); + s_ ^= (z_string & v_); // Remaining Z gates add a global phase from their action on the s string - phase += 4 * (AER::Utils::hamming_parity((z_string & ~(v)) & s)); + z_string = (z_string & ~v_) & s_; + phase += 4 * z_string.parity(); // Commute the x_string through the hadamard layer // Any remaining X gates update s - s ^= (x_string & ~(v)); + s_ ^= (x_string & ~v_); // New z gates add a global phase from their action on the s string - phase += 4 * (AER::Utils::hamming_parity((x_string & v) & s)); + x_string = (x_string & v_) & s_; + phase += 4 * (uint_t)x_string.parity(); // Update the global phase - omega.e = (omega.e + phase) % 8; + omega_.e = (omega_.e + phase) % 8; } void StabilizerState::Y(unsigned q) { Z(q); X(q); // Add a global phase of -i - omega.e = (omega.e + 2) % 8; + omega_.e = (omega_.e + 2) % 8; } void StabilizerState::RightCX(unsigned q, unsigned r) { - G[q] ^= G[r]; - F[r] ^= F[q]; - M[q] ^= M[r]; + G_[q] ^= G_[r]; + F_[r] ^= F_[q]; + M_[q] ^= M_[r]; } void StabilizerState::CX(unsigned q, unsigned r) { - isReadyMT = false; // we are going to change M and F - isReadyFT = false; - uint_fast64_t C = (one << q); - uint_fast64_t T = (one << r); bool b = false; - for (unsigned p = 0; p < n; p++) { - b = (b != ((M[p] & C) && (F[p] & T))); - G[p] ^= ((G[p] >> q) & one) * T; - F[p] ^= ((F[p] >> r) & one) * C; - M[p] ^= ((M[p] >> r) & one) * C; + for (unsigned p = 0; p < n_qubits_; p++) { + b = (b != (M_[p][q] && F_[p][r])); + G_[p].set_xor(r, G_[p][q]); + F_[p].set_xor(q, F_[p][r]); + M_[p].set_xor(q, M_[p][r]); } // update phase vector as // gamma[q] gets gamma[q] + gamma[r] + 2*b (mod 4) if (b) - gamma2 ^= C; - b = ((gamma1 & C) && (gamma1 & T)); - gamma1 ^= ((gamma1 >> r) & one) * C; - gamma2 ^= ((gamma2 >> r) & one) * C; + gamma2_.set_xor(q, one); + b = gamma1_[q] && gamma1_[r]; + gamma1_.set_xor(q, gamma1_[r]); + gamma2_.set_xor(q, gamma2_[r]); if (b) - gamma2 ^= C; + gamma2_.set_xor(q, one); } void StabilizerState::RightCZ(unsigned q, unsigned r) { - isReadyMT = false; // we are going to change M - M[q] ^= F[r]; - M[r] ^= F[q]; - gamma2 ^= (F[q] & F[r]); + M_[q] ^= F_[r]; + M_[r] ^= F_[q]; + gamma2_ ^= (F_[q] & F_[r]); } void StabilizerState::CZ(unsigned q, unsigned r) { - isReadyMT = false; // we are going to change M - uint_fast64_t C = (one << q); - uint_fast64_t T = (one << r); - for (unsigned p = 0; p < n; p++) { - M[p] ^= ((G[p] >> r) & one) * C; - M[p] ^= ((G[p] >> q) & one) * T; + for (unsigned p = 0; p < n_qubits_; p++) { + M_[p].set_xor(q, G_[p][r]); + M_[p].set_xor(r, G_[p][q]); } } -pauli_t StabilizerState::GetPauliX(uint_fast64_t x) { +pauli_t StabilizerState::GetPauliX(BitVector &x) { // make sure that M-transposed and F-transposed have been already computed - if (!isReadyMT) { - TransposeM(); - } - if (!isReadyFT) { - TransposeF(); - } - pauli_t R; - - for (unsigned pos = 0; pos < n; pos++) { - if (x & (one << pos)) { - pauli_t P1; // make P1=U_C^{-1} X_{pos} U_C - P1.e = 1 * ((gamma1 >> pos) & one); - P1.e += 2 * ((gamma2 >> pos) & one); - P1.X = FT[pos]; - P1.Z = MT[pos]; + pauli_t R(n_qubits_); + + for (unsigned pos = 0; pos < n_qubits_; pos++) { + if (x[pos]) { + pauli_t P1(n_qubits_); // make P1=U_C^{-1} X_{pos} U_C + P1.e = 1 * (uint_t)gamma1_[pos]; + P1.e += 2 * (uint_t)gamma2_[pos]; + for (uint_t i = 0; i < n_qubits_; i++) { + P1.X.set(i, F_[i][pos]); + P1.Z.set(i, M_[i][pos]); + } R *= P1; } } return R; } -scalar_t StabilizerState::Amplitude(uint_fast64_t x) { - // compute transposed matrices if needed - if (!isReadyMT) { - TransposeM(); - } - if (!isReadyFT) { - TransposeF(); - } +scalar_t StabilizerState::Amplitude(BitVector &x) { // compute Pauli U_C^{-1} X(x) U_C - P = GetPauliX(x); + P_ = GetPauliX(x); - if (!omega.eps) - return omega; // the state is zero + if (!omega_.eps) + return omega_; // the state is zero // now the amplitude = complex conjugate of // Z-part of P is absorbed into 0^n scalar_t amp; - amp.e = 2 * P.e; - int p = (int)AER::Utils::popcount(v); + amp.e = 2 * P_.e; + int p = (int)v_.popcount(); amp.p = -1 * p; // each Hadamard gate contributes 1/sqrt(2) bool isNonZero = true; - for (unsigned q = 0; q < n; q++) { - uint_fast64_t pos = (one << q); - if (v & pos) { - amp.e += - 4 * ((s & pos) && (P.X & pos)); // minus sign that comes from <1|H|1> + for (unsigned q = 0; q < n_qubits_; q++) { + if (v_[q]) { + amp.e += 4 * (s_[q] && P_.X[q]); // minus sign that comes from <1|H|1> } else { - isNonZero = (((P.X & pos) == (s & pos))); + isNonZero = (P_.X[q] == s_[q]); } if (!isNonZero) break; @@ -401,21 +362,23 @@ scalar_t StabilizerState::Amplitude(uint_fast64_t x) { } // multiply amp by omega - amp.p += omega.p; - amp.e = (amp.e + omega.e) % 8; + amp.p += omega_.p; + amp.e = (amp.e + omega_.e) % 8; return amp; } scalar_t StabilizerState::ProposeFlip(unsigned flip_pos) { // Q gets Pauli operator U_C^{-1} X_{flip_pos} U_C - Q.e = 1 * ((gamma1 >> flip_pos) & one); - Q.e += 2 * ((gamma2 >> flip_pos) & one); - Q.X = FT[flip_pos]; - Q.Z = MT[flip_pos]; - Q *= P; + Q_.e = 1 * gamma1_[flip_pos]; + Q_.e += 2 * gamma2_[flip_pos]; + for (uint_t i = 0; i < n_qubits_; i++) { + Q_.X.set(i, F_[i][flip_pos]); + Q_.Z.set(i, M_[i][flip_pos]); + } + Q_ *= P_; - if (!omega.eps) - return omega; // the state is zero + if (!omega_.eps) + return omega_; // the state is zero // the rest is the same as Amplitude() except that P becomes Q @@ -425,18 +388,16 @@ scalar_t StabilizerState::ProposeFlip(unsigned flip_pos) { // the rest is the same as Amplitude() except that P is replaced by Q scalar_t amp; - amp.e = 2 * Q.e; + amp.e = 2 * Q_.e; // each Hadamard gate contributes 1/sqrt(2) - amp.p = -1 * (int)(AER::Utils::popcount(v)); + amp.p = -1 * (int)(v_.popcount()); bool isNonZero = true; - for (unsigned q = 0; q < n; q++) { - uint_fast64_t pos = (one << q); - if (v & pos) { - amp.e += - 4 * ((s & pos) && (Q.X & pos)); // minus sign that comes from <1|H|1> + for (unsigned q = 0; q < n_qubits_; q++) { + if (v_[q]) { + amp.e += 4 * (s_[q] && Q_.X[q]); // minus sign that comes from <1|H|1> } else { - isNonZero = (((Q.X & pos) == (s & pos))); + isNonZero = (Q_.X[q] == s_[q]); } if (!isNonZero) break; @@ -451,82 +412,55 @@ scalar_t StabilizerState::ProposeFlip(unsigned flip_pos) { } // multiply amp by omega - amp.p += omega.p; - amp.e = (amp.e + omega.e) % 8; + amp.p += omega_.p; + amp.e = (amp.e + omega_.e) % 8; return amp; } -uint_fast64_t StabilizerState::Sample() { - uint_fast64_t x = zer; - for (unsigned q = 0; q < n; q++) { - bool w = !!(s & (one << q)); - w ^= ((v & (one << q)) && (rand() % 2)); +BitVector StabilizerState::Sample() { + BitVector x(n_qubits_); + for (unsigned q = 0; q < n_qubits_; q++) { + bool w = !!(s_[q]); + w ^= (v_[q] && (rand() % 2)); if (w) - x ^= G[q]; + x ^= G_[q]; } return x; } -uint_fast64_t StabilizerState::Sample(uint_fast64_t v_mask) { +BitVector StabilizerState::Sample(BitVector &v_mask) { // v_mask is a uniform random binary string we use to sample the bits // of v in a single step. - uint_fast64_t x = zer; - uint_fast64_t masked_v = v & v_mask; - for (unsigned q = 0; q < n; q++) { - bool w = !!(s & (one << q)); - w ^= !!(masked_v & (one << q)); + BitVector x(n_qubits_); + BitVector masked_v = v_ & v_mask; + for (unsigned q = 0; q < n_qubits_; q++) { + bool w = !!(s_[q]); + w ^= !!(masked_v[q]); if (w) - x ^= G[q]; + x ^= G_[q]; } return x; } -void StabilizerState::TransposeF() { - for (unsigned p = 0; p < n; p++) // look at p-th row of F - { - uint_fast64_t row = zer; - uint_fast64_t mask = (one << p); - for (unsigned q = 0; q < n; q++) // look at q-th column of F - if (F[q] & mask) - row ^= (one << q); - FT[p] = row; - } - isReadyFT = true; -} - -void StabilizerState::TransposeM() { - for (unsigned p = 0; p < n; p++) // look at p-th row of M - { - uint_fast64_t row = zer; - uint_fast64_t mask = (one << p); - for (unsigned q = 0; q < n; q++) // look at q-th column of M - if (M[q] & mask) - row ^= (one << q); - MT[p] = row; - } - isReadyMT = true; -} - -void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, - unsigned b) { +void StabilizerState::UpdateSvector(BitVector t, BitVector u, unsigned b) { // take care of the trivial case: t=u if (t == u) // multiply omega by (1+i^b)/sqrt(2) switch (b) { case 0: - omega.p += 1; - s = t; + omega_.p += 1; + s_ = t; return; case 1: - s = t; - omega.e = (omega.e + 1) % 8; + s_ = t; + omega_.e = (omega_.e + 1) % 8; return; case 2: - s = t; - omega.eps = 0; + s_ = t; + omega_.eps = 0; return; case 3: - s = t; - omega.e = (omega.e + 7) % 8; + s_ = t; + omega_.e = (omega_.e + 7) % 8; return; default: // we should not get here @@ -535,42 +469,38 @@ void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, } // now t and u are distinct - isReadyFT = false; // below we are going to change F and M - isReadyMT = false; // naming of variables roughly follows Section IIIA - uint_fast64_t ut = u ^ t; - uint_fast64_t nu0 = (~v) & ut; - uint_fast64_t nu1 = v & ut; + BitVector ut = u ^ t; + BitVector nu0 = (~v_) & ut; + BitVector nu1 = v_ & ut; // b %= 4; unsigned q = 0; - uint_fast64_t qpos = zer; - if (nu0) { + if (!nu0.is_zero()) { // the subset nu0 is non-empty // find the first element of nu0 q = 0; - while (!(nu0 & (one << q))) + while (!nu0[q]) q++; - qpos = (one << q); - if (!(nu0 & qpos)) { + if (!nu0[q]) { throw std::logic_error( "Failed to find first bit of nu despite being non-empty."); } // if nu0 has size >1 then multiply U_C on the right by the first half of // the circuit VC - nu0 ^= qpos; // set q-th bit to zero - if (nu0) - for (unsigned q1 = q + 1; q1 < n; q1++) - if (nu0 & (one << q1)) + nu0.set_xor(q, one); // set q-th bit to zero + if (!nu0.is_zero()) + for (unsigned q1 = q + 1; q1 < n_qubits_; q1++) + if (nu0[q1]) RightCX(q, q1); // if nu1 has size >0 then apply the second half of the circuit VC - if (nu1) - for (unsigned q1 = 0; q1 < n; q1++) - if (nu1 & (one << q1)) + if (!nu1.is_zero()) + for (unsigned q1 = 0; q1 < n_qubits_; q1++) + if (nu1[q1]) RightCZ(q, q1); } // if (nu0) @@ -578,19 +508,18 @@ void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, // if we got here when nu0 is empty // find the first element of nu1 q = 0; - while (!(nu1 & (one << q))) + while (!nu1[q]) q++; - qpos = (one << q); - if (!(nu1 & qpos)) { + if (!nu1[q]) { throw std::logic_error("Expect at least one non-zero element in nu1.\n"); } // if nu1 has size >1 then apply the circuit VC - nu1 ^= qpos; - if (nu1) - for (unsigned q1 = q + 1; q1 < n; q1++) - if (nu1 & (one << q1)) + nu1.set_xor(q, one); + if (!nu1.is_zero()) + for (unsigned q1 = q + 1; q1 < n_qubits_; q1++) + if (nu1[q1]) RightCX(q1, q); } // if (nu0) else @@ -598,14 +527,14 @@ void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, // update the initial state // if t_q=1 then switch t_q and u_q // uint_fast64_t y,z; - if (t & qpos) { - s = u; - omega.e = (omega.e + 2 * b) % 8; + if (t[q]) { + s_ = u; + omega_.e = (omega_.e + 2 * b) % 8; b = (4 - b) % 4; - if (!!(u & qpos)) { + if (!!(u[q])) { } } else - s = t; + s_ = t; // change the order of H and S gates // H^{a} S^{b} |+> = eta^{e1} S^{e2} H^{e3} |e4> @@ -629,7 +558,7 @@ void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, // e3 = not(a) + a * (b mod 2) // e4 = not(a)*(b>=2) + a*( (b==1) || (b==2) ) - bool a = ((v & qpos) > 0); + bool a = v_[q]; unsigned e1 = a * (b % 2) * (3 * b - 2); unsigned e2 = b % 2; bool e3 = ((!a) != (a && ((b % 2) > 0))); @@ -637,15 +566,13 @@ void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, // update CH-form // set q-th bit of s to e4 - s &= ~qpos; - s ^= e4 * qpos; + s_.set(q, e4); // set q-th bit of v to e3 - v &= ~qpos; - v ^= e3 * qpos; + v_.set(q, e3); // update the scalar factor omega - omega.e = (omega.e + e1) % 8; + omega_.e = (omega_.e + e1) % 8; // multiply the C-layer on the right by S^{e2} on the q-th qubit if (e2) @@ -653,96 +580,102 @@ void StabilizerState::UpdateSvector(uint_fast64_t t, uint_fast64_t u, } void StabilizerState::H(unsigned q) { - isReadyMT = false; // we are going to change M and F - isReadyFT = false; // extract the q-th row of F,G,M - uint_fast64_t rowF = zer; - uint_fast64_t rowG = zer; - uint_fast64_t rowM = zer; - for (unsigned j = 0; j < n; j++) { - uint_fast64_t pos = (one << j); - rowF ^= pos * ((F[j] >> q) & one); - rowG ^= pos * ((G[j] >> q) & one); - rowM ^= pos * ((M[j] >> q) & one); + BitVector rowF(n_qubits_, true); + BitVector rowG(n_qubits_, true); + BitVector rowM(n_qubits_, true); + for (unsigned j = 0; j < n_qubits_; j++) { + rowF.set_xor(j, F_[j][q]); + rowG.set_xor(j, G_[j][q]); + rowM.set_xor(j, M_[j][q]); } // after commuting H through the C and H laters it maps |s> to a state // sqrt(0.5)*[ (-1)^alpha |t> + i^{gamma[p]} (-1)^beta |u> ] // // compute t,s,alpha,beta - uint_fast64_t t = s ^ (rowG & v); - uint_fast64_t u = s ^ (rowF & (~v)) ^ (rowM & v); + BitVector t = s_ ^ (rowG & v_); + BitVector u = s_ ^ (rowF & (~v_)) ^ (rowM & v_); - unsigned alpha = AER::Utils::popcount(rowG & (~v) & s); - unsigned beta = - AER::Utils::popcount((rowM & (~v) & s) ^ (rowF & v & (rowM ^ s))); + BitVector pc = rowG & (~v_) & s_; + unsigned alpha = pc.popcount(); + pc = (rowM & (~v_) & s_) ^ (rowF & v_ & (rowM ^ s_)); + unsigned beta = pc.popcount(); if (alpha % 2) - omega.e = (omega.e + 4) % 8; + omega_.e = (omega_.e + 4) % 8; // get the phase gamma[q] - unsigned phase = ((gamma1 >> q) & one) + 2 * ((gamma2 >> q) & one); + unsigned phase = (unsigned)gamma1_[q] + 2 * (unsigned)gamma2_[q]; unsigned b = (phase + 2 * alpha + 2 * beta) % 4; // now the initial state is sqrt(0.5)*(|t> + i^b |u>) // take care of the trivial case if (t == u) { - s = t; + s_ = t; if (!((b == 1) || (b == 3))) // otherwise the state is not normalized { throw std::logic_error( "State is not properly normalised, b should be 1 or 3.\n"); } if (b == 1) - omega.e = (omega.e + 1) % 8; + omega_.e = (omega_.e + 1) % 8; else - omega.e = (omega.e + 7) % 8; + omega_.e = (omega_.e + 7) % 8; } else UpdateSvector(t, u, b); } void StabilizerState::MeasurePauli(pauli_t PP) { // compute Pauli R = U_C^{-1} P U_C - pauli_t R; + pauli_t R(n_qubits_); R.e = PP.e; - for (unsigned j = 0; j < n; j++) - if ((PP.X >> j) & one) { + for (unsigned j = 0; j < n_qubits_; j++) + if (PP.X[j]) { // multiply R by U_C^{-1} X_j U_C // extract the j-th rows of F and M - uint_fast64_t rowF = zer; - uint_fast64_t rowM = zer; - for (unsigned i = 0; i < n; i++) { - rowF ^= (one << i) * ((F[i] >> j) & one); - rowM ^= (one << i) * ((M[i] >> j) & one); + BitVector rowF(n_qubits_, true); + BitVector rowM(n_qubits_, true); + for (unsigned i = 0; i < n_qubits_; i++) { + rowF.set_xor(i, F_[i][j]); + rowM.set_xor(i, M_[i][j]); } - R.e += 2 * AER::Utils::popcount( - R.Z & rowF); // extra sign from Pauli commutation + BitVector pc = R.Z; + pc &= rowF; + R.e += 2 * pc.popcount(); // extra sign from Pauli commutation R.Z ^= rowM; R.X ^= rowF; - R.e += ((gamma1 >> j) & one) + 2 * ((gamma2 >> j) & one); + R.e += (uint_t)gamma1_[j] + 2 * (uint_t)gamma2_[j]; } - for (unsigned q = 0; q < n; q++) - R.Z ^= (one << q) * (AER::Utils::popcount(PP.Z & G[q]) % 2); + for (unsigned q = 0; q < n_qubits_; q++) { + BitVector pc = PP.Z; + pc &= G_[q]; + R.Z.set_xor(q, (pc.popcount() & 1)); + } // now R=U_C^{-1} PP U_C // next conjugate R by U_H - uint_fast64_t tempX = ((~v) & R.X) ^ (v & R.Z); - uint_fast64_t tempZ = ((~v) & R.Z) ^ (v & R.X); + BitVector tempX = ((~v_) & R.X) ^ (v_ & R.Z); + BitVector tempZ = ((~v_) & R.Z) ^ (v_ & R.X); // the sign flips each time a Hadamard hits Y on some qubit - R.e = (R.e + 2 * AER::Utils::popcount(v & R.X & R.Z)) % 4; + BitVector pc = v_; + pc &= R.X; + pc &= R.Z; + R.e = (R.e + 2 * pc.popcount()) % 4; R.X = tempX; R.Z = tempZ; // now the initial state |s> becomes 0.5*(|s> + R |s>) = 0.5*(|s> + i^b |s ^ // R.X>) - unsigned b = (R.e + 2 * AER::Utils::popcount(R.Z & s)) % 4; - UpdateSvector(s, s ^ R.X, b); + BitVector sz = s_; + sz &= R.Z; + unsigned b = (R.e + 2 * (unsigned)sz.popcount()) % 4; + BitVector sr = s_; + sr ^= R.X; + UpdateSvector(s_, sr, b); // account for the extra factor sqrt(1/2) - omega.p -= 1; - - isReadyMT = false; // we have changed change M and F - isReadyFT = false; + omega_.p -= 1; } void StabilizerState::MeasurePauliProjector( @@ -751,95 +684,100 @@ void StabilizerState::MeasurePauliProjector( { for (uint_fast64_t i = 0; i < generators.size(); i++) { this->MeasurePauli(generators[i]); - if (omega.eps == 0) { + if (omega_.eps == 0) { break; } } } -scalar_t StabilizerState::InnerProduct(const uint_fast64_t &A_diag1, - const uint_fast64_t &A_diag2, - const std::vector &A) { - uint_fast64_t K_diag1 = zer, K_diag2 = zer, J_diag1 = gamma1, - J_diag2 = gamma2; - std::vector J(n, zer); - std::vector K(n, zer); - std::vector placeholder(n, zer); - if (!isReadyMT) { - TransposeM(); - } - if (!isReadyFT) { - TransposeF(); - } +scalar_t StabilizerState::InnerProduct(const BitVector &A_diag1, + const BitVector &A_diag2, + const std::vector &A) { + BitVector K_diag1(n_qubits_, true), K_diag2(n_qubits_, true); + BitVector J_diag1 = gamma1_; + BitVector J_diag2 = gamma2_; + std::vector J(n_qubits_, BitVector(n_qubits_, true)); + std::vector K(n_qubits_, BitVector(n_qubits_, true)); + std::vector placeholder(n_qubits_, BitVector(n_qubits_, true)); + // Setup the J matrix - for (size_t i = 0; i < n; i++) { - for (size_t j = i; j < n; j++) { - if (AER::Utils::hamming_parity(MT[i] & FT[j])) { - J[i] |= (one << j); - J[j] |= (one << i); + for (size_t i = 0; i < n_qubits_; i++) { + BitVector mt(n_qubits_, true); + for (uint_t k = 0; k < n_qubits_; k++) + mt.set(k, M_[k][i]); + + for (size_t j = i; j < n_qubits_; j++) { + BitVector ft(n_qubits_, true); + for (uint_t k = 0; k < n_qubits_; k++) + ft.set(k, F_[k][j] & mt[k]); + if (ft.parity()) { + J[i].set_or(j, true); + J[j].set_or(i, true); } } } + // Calculate the matrix J = A+J J_diag2 = (A_diag2 ^ J_diag2 ^ (A_diag1 & J_diag1)); J_diag1 = (A_diag1 ^ J_diag1); - for (size_t i = 0; i < n; i++) { + for (size_t i = 0; i < n_qubits_; i++) { J[i] ^= A[i]; } // placeholder = J*G) - for (size_t i = 0; i < n; i++) { + for (size_t i = 0; i < n_qubits_; i++) { // Grab column i of J, it's symmetric - uint_fast64_t col_i = J[i]; - for (size_t j = 0; j < n; j++) { - if (AER::Utils::hamming_parity(col_i & G[j])) { - placeholder[j] |= (one << i); + BitVector col_i = J[i]; + for (size_t j = 0; j < n_qubits_; j++) { + BitVector p = col_i; + p &= G_[j]; + if (p.parity()) { + placeholder[j].set_or(i, true); } } } // K = GT*placeholder - for (size_t i = 0; i < n; i++) { - uint_fast64_t col_i = G[i]; - uint_fast64_t shift = (one << i); - for (size_t j = i; j < n; j++) { - if (AER::Utils::hamming_parity(col_i & placeholder[j])) { - K[j] |= shift; - K[i] |= (one << j); + for (size_t i = 0; i < n_qubits_; i++) { + BitVector col_i(n_qubits_); + col_i = G_[i]; + for (size_t j = i; j < n_qubits_; j++) { + BitVector p = col_i & placeholder[j]; + if (p.parity()) { + K[j].set_or(i, true); + K[i].set_or(j, true); } } - for (size_t r = 0; r < n; r++) { - if ((col_i >> r) & one) { - uint_fast64_t one_bit = (J_diag1 >> r) & one; - if ((K_diag1 >> i) & one_bit) { - K_diag2 ^= shift; + for (size_t r = 0; r < n_qubits_; r++) { + if (col_i[r]) { + if (K_diag1[i] & J_diag1[r]) { + K_diag2.set_xor(i, true); } - K_diag1 ^= one_bit * shift; - K_diag2 ^= ((J_diag2 >> r) & one) * shift; + K_diag1.set_xor(i, J_diag1[r]); + K_diag2.set_xor(i, J_diag2[r]); } - for (size_t k = r + 1; k < n; k++) { - if ((J[k] >> r) & (col_i >> r) & (col_i >> k) & one) { - K_diag2 ^= shift; + for (size_t k = r + 1; k < n_qubits_; k++) { + if (J[k][r] & col_i[r] & col_i[k]) { + K_diag2.set_xor(i, true); } } } } unsigned col = 0; - QuadraticForm q(AER::Utils::popcount(v)); + QuadraticForm q(v_.popcount()); // We need to setup a quadratic form to evaluate the Exponential Sum - for (size_t i = 0; i < n; i++) { - if ((v >> i) & one) { - uint_fast64_t shift = (one << col); + for (size_t i = 0; i < n_qubits_; i++) { + if (v_[i]) { // D = Diag(K(1,1)) + 2*[s + s*K](1) // J = K(1,1); - q.D1 ^= ((K_diag1 >> i) & one) * shift; - q.D2 ^= ((((K_diag2 >> i) ^ (s >> i)) & one) ^ - AER::Utils::hamming_parity(K[i] & s)) * - shift; + q.D1.set_xor(col, K_diag1[i]); + BitVector p = K[i]; + p &= s_; + q.D2.set_xor(col, K_diag2[i] ^ s_[i] ^ p.parity()); // q.D2 ^= ((s >> i) & one) * shift; // q.D2 ^= AER::Utils::hamming_parity(K[i] & s) * shift; unsigned row = 0; - for (size_t j = 0; j < n; j++) { - if ((v >> j) & one) { - q.J[col] |= (((K[i] >> j) & one) << row); + for (size_t j = 0; j < n_qubits_; j++) { + if (v_[j]) { + q.J[col].set_or(row, K[i][j]); row++; } } @@ -847,13 +785,14 @@ scalar_t StabilizerState::InnerProduct(const uint_fast64_t &A_diag1, } } // Q = 4* (s.v) + sKs - q.Q = AER::Utils::hamming_parity(s & v) * 4; - for (size_t i = 0; i < n; i++) { - if ((s >> i) & one) { - q.Q = (q.Q + 4 * ((K_diag2 >> i) & one) + 2 * ((K_diag1 >> i) & one)) % + BitVector p = s_ & v_; + q.Q = (uint_t)p.parity() * 4; + for (size_t i = 0; i < n_qubits_; i++) { + if (s_[i]) { + q.Q = (q.Q + 4 * (uint_t)K_diag2[i] + 2 * (uint_t)K_diag1[i]) % 8; // + 2*((K_diag1 >> i) & one))%8; - for (size_t j = i + 1; j < n; j++) { - if ((s >> j) & (K[j] >> i) & one) { + for (size_t j = i + 1; j < n_qubits_; j++) { + if (s_[j] & K[j][i]) { q.Q ^= 4; } } @@ -861,9 +800,9 @@ scalar_t StabilizerState::InnerProduct(const uint_fast64_t &A_diag1, } scalar_t amp = q.ExponentialSum(); // Reweight by 2^{-(n+|v|)}/2 - amp.p -= (n + AER::Utils::popcount(v)); + amp.p -= (n_qubits_ + v_.popcount()); // We need to further multiply by omega* - scalar_t psi_amp(omega); + scalar_t psi_amp(omega_); psi_amp.conjugate(); amp *= psi_amp; return amp; @@ -871,9 +810,9 @@ scalar_t StabilizerState::InnerProduct(const uint_fast64_t &A_diag1, double NormEstimate(std::vector &states, const std::vector> &phases, - const std::vector &Samples_d1, - const std::vector &Samples_d2, - const std::vector> &Samples) { + const std::vector &Samples_d1, + const std::vector &Samples_d2, + const std::vector> &Samples) { // Norm estimate for a state |psi> = \sum_{i} c_{i}|phi_{i}> double xi = 0; unsigned L = Samples_d1.size(); @@ -903,13 +842,12 @@ double NormEstimate(std::vector &states, return std::pow(2, states[0].NQubits()) * (xi / L); } -double -ParallelNormEstimate(std::vector &states, - const std::vector> &phases, - const std::vector &Samples_d1, - const std::vector &Samples_d2, - const std::vector> &Samples, - int n_threads) { +double ParallelNormEstimate(std::vector &states, + const std::vector> &phases, + const std::vector &Samples_d1, + const std::vector &Samples_d2, + const std::vector> &Samples, + int n_threads) { double xi = 0; unsigned L = Samples_d1.size(); unsigned chi = states.size(); @@ -938,4 +876,7 @@ ParallelNormEstimate(std::vector &states, } } // namespace CHSimulator +//------------------------------------------------------------------------------ +} // end namespace AER +//------------------------------------------------------------------------------ #endif diff --git a/src/simulators/extended_stabilizer/chlib/core.hpp b/src/simulators/extended_stabilizer/chlib/core.hpp index e4e18dba7a..af6a07c575 100644 --- a/src/simulators/extended_stabilizer/chlib/core.hpp +++ b/src/simulators/extended_stabilizer/chlib/core.hpp @@ -22,8 +22,10 @@ #include #include +#include "framework/bitvector.hpp" #include "framework/utils.hpp" +namespace AER { namespace CHSimulator { static const std::array RE_PHASE = {1, 1, 0, -1, -1, -1, 0, 1}; @@ -125,12 +127,13 @@ struct scalar_t { // Below we represent n-bit strings by uint64_t integers struct pauli_t { // n-qubit Pauli operators: i^e * X(x) * Z(z) - uint_fast64_t X; // n-bit string - uint_fast64_t Z; // n-bit string - unsigned e = 0; // takes values 0,1,2,3 + BitVector X; // n-bit string + BitVector Z; // n-bit string + unsigned e = 0; // takes values 0,1,2,3 // constructor makes the identity Pauli operator pauli_t(); + pauli_t(uint_t n); pauli_t(const pauli_t &p) = default; // multiplication of Pauli operators @@ -140,9 +143,9 @@ struct pauli_t { struct QuadraticForm { unsigned int n; int Q; - uint_fast64_t D1; - uint_fast64_t D2; - std::vector J; + BitVector D1; + BitVector D2; + std::vector J; scalar_t ExponentialSum(); QuadraticForm(unsigned n_qubits); QuadraticForm(const QuadraticForm &rhs); @@ -179,11 +182,14 @@ scalar_t scalar_t::operator*(const scalar_t &rhs) const { return out; } -pauli_t::pauli_t() : X(zer), Z(zer) {} +pauli_t::pauli_t() {} +pauli_t::pauli_t(uint_t n) : X(n, true), Z(n, true) {} pauli_t &pauli_t::operator*=(const pauli_t &rhs) { - unsigned overlap = - AER::Utils::popcount(Z & rhs.X); // commute rhs.X to the left + BitVector pc = Z; + pc &= rhs.X; + unsigned overlap = pc.popcount(); + // AER::Utils::popcount(Z & rhs.X); // commute rhs.X to the left X ^= rhs.X; Z ^= rhs.Z; e = (e + rhs.e + 2 * overlap) % 4; @@ -198,11 +204,8 @@ class QubitException : public std::exception { } qubex; QuadraticForm::QuadraticForm(unsigned int n_qubits) - : n(n_qubits), Q(0), D1(zer), D2(zer), J(n_qubits, zer) { - if (n > 63) { - throw qubex; - } -} + : n(n_qubits), Q(0), D1(n_qubits, true), D2(n_qubits, true), + J(n_qubits, BitVector(n_qubits, true)) {} QuadraticForm::QuadraticForm(const QuadraticForm &rhs) : J(rhs.n, zer) { n = rhs.n; @@ -246,13 +249,13 @@ std::ostream &operator<<(std::ostream &os, const QuadraticForm &q) { os << "Q: " << q.Q << std::endl; os << "D:"; for (size_t i = 0; i < q.n; i++) { - os << " " << (2 * (2 * ((q.D2 >> i) & one) + ((q.D1 >> i) & one))); + os << " " << (2 * (2 * q.D2[i] + q.D1[i])); } os << std::endl; os << "J:\n"; for (size_t i = 0; i < q.n; i++) { for (size_t j = 0; j < q.n; j++) { - os << ((q.J[i] >> j) & one) << " "; + os << (q.J[i][j]) << " "; } os << "\n"; } @@ -302,28 +305,29 @@ scalar_t QuadraticForm::ExponentialSum() { // Lreal, Limag define the linear part of the Z2-valued forms // M defines the quadratic part of the Z2-valued forms // each column of M is represented by long integer - uint_fast64_t Lreal, c = zer; - std::vector M(n + 1, zer); - Lreal = D2; + BitVector Lreal(m, true), c(m, true); + std::vector M(n + 1, BitVector(m, true)); + for (size_t i = 0; i < D2.length(); i++) + Lreal(i) = D2(i); for (size_t i = 0; i < n; i++) { - if ((D1 >> i) & one) { - c |= (one << i); - M[n] |= (one << i); + if (D1[i]) { + c.set_or(i, true); + M[n].set_or(i, true); } } - uint_fast64_t Limag = Lreal; - Limag ^= (one << n); + BitVector Limag = Lreal; + Limag.set_xor(n, true); for (size_t i = 0; i < n; i++) { for (size_t j = (i + 1); j < n; j++) { - bool x = !!((J[i] >> j) & one); - bool c1 = !!((c >> i) & one); - bool c2 = !!((c >> j) & one); + bool x = !!(J[i][j]); + bool c1 = !!(c[j]); + bool c2 = !!(c[j]); if (x ^ (c1 & c2)) { - M[j] ^= (one << i); + M[j].set_xor(i, true); } } } - uint_fast64_t active = zer; + BitVector active(m, true); active = ~(active); int nActive = m; int firstActive = 0; @@ -331,7 +335,7 @@ scalar_t QuadraticForm::ExponentialSum() { // find the first active variable int i1; for (i1 = firstActive; i1 < m; i1++) { - if ((active >> i1) & one) { + if (active[i1]) { break; } } @@ -340,13 +344,13 @@ scalar_t QuadraticForm::ExponentialSum() { int i2; bool isFound = false; for (i2 = 0; i2 < m; i2++) { - isFound = (((M[i1] >> i2) & one) != ((M[i2] >> i1) & one)); + isFound = (M[i1][i2] != M[i2][i1]); if (isFound) { break; } } - bool L1real = !!(((Lreal >> i1) & one) ^ ((M[i1] >> i1) & one)); - bool L1imag = !!(((Limag >> i1) & one) ^ ((M[i1] >> i1) & one)); + bool L1real = !!(Lreal[i1] ^ M[i1][i1]); + bool L1imag = !!(Limag[i1] ^ M[i1][i1]); // take care of the trivial cases if (!isFound) { // the form is linear in the variable i1 @@ -361,17 +365,15 @@ scalar_t QuadraticForm::ExponentialSum() { pow2_imag++; nActive--; - uint_fast64_t not_shift; - not_shift = ~(one << i1); - active &= not_shift; + active.set(i1, false); // set column i1 to zero - M[i1] = zer; + M[i1].zero(); // set row i1 to zero for (int j = 0; j < m; j++) { - M[j] &= not_shift; + M[j].set(i1, false); } - Lreal &= not_shift; - Limag &= not_shift; + Lreal.set(i1, false); + Limag.set(i1, false); if (isZero_real && isZero_imag) { amp.eps = 0; @@ -379,32 +381,32 @@ scalar_t QuadraticForm::ExponentialSum() { } continue; } // trivial cases - bool L2real = !!(((Lreal >> i2) & one) ^ ((M[i2] >> i2) & one)); - bool L2imag = !!(((Limag >> i2) & one) ^ ((M[i2] >> i2) & one)); - Lreal &= ~(one << i1); - Lreal &= ~(one << i2); - Limag &= ~(one << i1); - Limag &= ~(one << i2); + bool L2real = !!(Lreal[i2] ^ M[i2][i2]); + bool L2imag = !!(Limag[i2] ^ M[i2][i2]); + Lreal.set(i1, false); + Lreal.set(i2, false); + Limag.set(i1, false); + Limag.set(i2, false); // Extract rows i1 and i2 of M - uint_fast64_t m1 = zer; - uint_fast64_t m2 = zer; + BitVector m1(m); + BitVector m2(m); for (int j = 0; j < m; j++) { - m1 ^= ((M[j] >> i1) & one) << j; - m2 ^= ((M[j] >> i2) & one) << j; + m1.set_xor(j, M[j][i1]); + m2.set_xor(j, M[j][i2]); } m1 ^= M[i1]; m2 ^= M[i2]; - m1 &= ~(one << i1); - m1 &= ~(one << i2); - m2 &= ~(one << i1); - m2 &= ~(one << i2); + m1.set(i1, false); + m1.set(i2, false); + m2.set(i1, false); + m2.set(i2, false); // set columns i1, i2 to zero - M[i1] = zer; - M[i2] = zer; + M[i1].zero(); + M[i2].zero(); // set rows i1,i2 to zero for (int j = 0; j < m; j++) { - M[j] &= ~(one << i1); - M[j] &= ~(one << i2); + M[j].set(i1, false); + M[j].set(i2, false); } // update the std::vectors Lreal, Limag @@ -430,12 +432,12 @@ scalar_t QuadraticForm::ExponentialSum() { } // update matrix M for (int j = 0; j < m; j++) { - if ((m2 >> j) & one) { + if (m2[j]) { M[j] ^= m1; } } - active &= ~(one << i1); - active &= ~(one << i2); + active.set(i1, false); + active.set(i2, false); nActive -= 2; } // int J0 = Q; @@ -490,4 +492,7 @@ void Print(std::vector A, unsigned n) { } } // namespace CHSimulator +//------------------------------------------------------------------------------ +} // end namespace AER +//------------------------------------------------------------------------------ #endif diff --git a/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp b/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp index be6a8af609..c4bb422302 100644 --- a/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp +++ b/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp @@ -417,7 +417,7 @@ void State::apply_ops(InputIterator first, InputIterator last, std::vector State::sample_measure(const reg_t &qubits, uint_t shots, RngEngine &rng) { - std::vector output_samples; + std::vector output_samples; if (BaseState::qreg_.get_num_states() == 1) { output_samples = BaseState::qreg_.stabilizer_sampler(shots, rng); } else { @@ -441,10 +441,10 @@ std::vector State::sample_measure(const reg_t &qubits, uint_t shots, } std::vector all_samples; all_samples.reserve(shots); - for (uint_t sample : output_samples) { + for (BitVector &sample : output_samples) { reg_t sample_bits(qubits.size(), 0ULL); for (size_t i = 0; i < qubits.size(); i++) { - if ((sample >> qubits[i]) & 1ULL) { + if (sample[qubits[i]]) { sample_bits[i] = 1ULL; } } @@ -539,7 +539,7 @@ void State::apply_stabilizer_circuit(InputIterator first, InputIterator last, void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, const reg_t &cregister, RngEngine &rng) { - uint_t out_string; + BitVector out_string; // Flag if the Pauli projector is applied already as part of the sampling bool do_projector_correction = true; // Prepare an output register for the qubits we are measurig @@ -565,11 +565,12 @@ void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, } if (do_projector_correction) { // We prepare the Pauli projector corresponding to the measurement result - std::vector paulis(qubits.size(), chpauli_t()); + std::vector paulis(qubits.size(), + chpauli_t(BaseState::qreg_.get_n_qubits())); for (uint_t i = 0; i < qubits.size(); i++) { // Create a Pauli projector onto 1+-Z_{i} on qubit i - paulis[i].Z = (1ULL << qubits[i]); - if ((out_string >> qubits[i]) & 1ULL) { + paulis[i].Z.set(qubits[i], true); + if (out_string[qubits[i]]) { // Additionally, store the output bit for this qubit paulis[i].e = 2; } @@ -579,7 +580,7 @@ void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, } for (uint_t i = 0; i < qubits.size(); i++) { // Create a Pauli projector onto 1+-Z_{i} on qubit i - if ((out_string >> qubits[i]) & 1ULL) { + if (out_string[qubits[i]]) { // Additionally, store the output bit for this qubit outcome[i] = 1ULL; } @@ -589,7 +590,7 @@ void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, } void State::apply_reset(const reg_t &qubits, AER::RngEngine &rng) { - uint_t measure_string; + BitVector measure_string; const int_t NUM_STATES = BaseState::qreg_.get_num_states(); if (BaseState::qreg_.get_num_states() == 1) { measure_string = BaseState::qreg_.stabilizer_sampler(rng); @@ -598,12 +599,12 @@ void State::apply_reset(const reg_t &qubits, AER::RngEngine &rng) { BaseState::qreg_.metropolis_estimation(metropolis_mixing_steps_, rng); } - std::vector paulis(qubits.size(), chpauli_t()); + std::vector paulis(qubits.size(), + chpauli_t(BaseState::qreg_.get_n_qubits())); for (size_t i = 0; i < qubits.size(); i++) { uint_t qubit = qubits[i]; - uint_t shift = 1ULL << qubit; - paulis[i].Z = shift; - if (!!(measure_string & shift)) { + paulis[i].Z.set(qubits[i], true); + if (!!(measure_string[qubits[i]])) { paulis[i].e = 2; } } @@ -613,7 +614,7 @@ void State::apply_reset(const reg_t &qubits, AER::RngEngine &rng) { num_threads(BaseState::threads_) for (int_t i = 0; i < NUM_STATES; i++) { for (auto qubit : qubits) { - if ((measure_string >> qubit) & 1ULL) { + if (measure_string[qubit]) { BaseState::qreg_.apply_x(qubit, i); } } @@ -777,20 +778,20 @@ double State::expval_pauli(const reg_t &qubits, const std::string &pauli, auto state_cpy = BaseState::qreg_; auto phi_norm = state_cpy.norm_estimation(norm_estimation_samples_, norm_estimation_repetitions_, rng); - std::vector paulis(1, chpauli_t()); + std::vector paulis(1, chpauli_t(BaseState::qreg_.get_n_qubits())); for (uint_t pos = 0; pos < qubits.size(); ++pos) { switch (pauli[pauli.size() - 1 - pos]) { case 'I': break; case 'X': - paulis[0].X += (1ULL << qubits[pos]); + paulis[0].X.set(qubits[pos], true); break; case 'Y': - paulis[0].X += (1ULL << qubits[pos]); - paulis[0].Z += (1ULL << qubits[pos]); + paulis[0].X.set(qubits[pos], true); + paulis[0].Z.set(qubits[pos], true); break; case 'Z': - paulis[0].Z += (1ULL << qubits[pos]); + paulis[0].Z.set(qubits[pos], true); break; default: { std::stringstream msg; diff --git a/src/simulators/extended_stabilizer/gates.hpp b/src/simulators/extended_stabilizer/gates.hpp index 3df76b37eb..605d1f0808 100644 --- a/src/simulators/extended_stabilizer/gates.hpp +++ b/src/simulators/extended_stabilizer/gates.hpp @@ -35,6 +35,7 @@ constexpr double pow(double x, int y) { return std::pow(x, y); } } // namespace std #endif +namespace AER { namespace CHSimulator { using uint_t = uint_fast64_t; using complex_t = std::complex; @@ -229,5 +230,8 @@ inline double u1_extent(double lambda) { } } // namespace CHSimulator +//------------------------------------------------------------------------------ +} // end namespace AER +//------------------------------------------------------------------------------ #endif From 414a97abafc6fe4e43aa1d73988715134c4e0b1a Mon Sep 17 00:00:00 2001 From: Jun Doi Date: Tue, 20 Feb 2024 14:11:45 +0900 Subject: [PATCH 2/5] fix qiskit-aer/library/save-instructions/save_data.py --- qiskit_aer/library/save_instructions/save_data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit_aer/library/save_instructions/save_data.py b/qiskit_aer/library/save_instructions/save_data.py index a98e88e12f..7a1b6b77b9 100644 --- a/qiskit_aer/library/save_instructions/save_data.py +++ b/qiskit_aer/library/save_instructions/save_data.py @@ -67,7 +67,7 @@ def assemble(self): instr.label = self._label return instr - def inverse(self): + def inverse(self, annotated=False): """Special case. Return self.""" return copy.copy(self) From e93bfb1e5d51ce4e80b2e49d1ae7e3122c7bb99d Mon Sep 17 00:00:00 2001 From: Jun Doi Date: Tue, 20 Feb 2024 14:26:29 +0900 Subject: [PATCH 3/5] version of qiskit in setup.py to 1.0.0 --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index ba1ef7dcb6..bf08fd7ce1 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ extras_requirements = {"dask": ["dask", "distributed"]} requirements = [ - "qiskit>=0.45.0", + "qiskit>=1.0.0", "numpy>=1.16.3", "scipy>=1.0", "psutil>=5", From 00b76505de11d9599c341acbf6d903b524a114ab Mon Sep 17 00:00:00 2001 From: Jun Doi Date: Tue, 20 Feb 2024 14:48:25 +0900 Subject: [PATCH 4/5] Revert "version of qiskit in setup.py to 1.0.0" This reverts commit e93bfb1e5d51ce4e80b2e49d1ae7e3122c7bb99d. --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index bf08fd7ce1..ba1ef7dcb6 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ extras_requirements = {"dask": ["dask", "distributed"]} requirements = [ - "qiskit>=1.0.0", + "qiskit>=0.45.0", "numpy>=1.16.3", "scipy>=1.0", "psutil>=5", From 23e03e9e1695f85ed0705504966d503023ebdc42 Mon Sep 17 00:00:00 2001 From: Jun Doi Date: Wed, 21 Feb 2024 14:46:35 +0900 Subject: [PATCH 5/5] change versions of wheel builds --- .github/workflows/build.yml | 12 ++++++------ tox.ini | 1 - 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 9b3540e299..e74e3fe1bb 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -123,13 +123,13 @@ jobs: matrix: os: ["macOS-latest", "ubuntu-latest", "windows-2019"] steps: - - uses: actions/checkout@v2 - - name: Set up Python Python 3.8 - uses: actions/setup-python@v2 + - uses: actions/checkout@v3 + - name: Set up Python Python 3.10 + uses: actions/setup-python@v4 with: - python-version: 3.8 + python-version: 3.10 - name: Add msbuild to PATH - uses: microsoft/setup-msbuild@v1.0.2 + uses: microsoft/setup-msbuild@v2 if: runner.os == 'Windows' - name: Install deps run: python -m pip install -U cibuildwheel==2.16.2 @@ -137,7 +137,7 @@ jobs: env: AER_CMAKE_OPENMP_BUILD: 1 run: cibuildwheel --output-dir wheelhouse - - uses: actions/upload-artifact@v2 + - uses: actions/upload-artifact@v3 with: path: ./wheelhouse/*.whl wheel-arm64: diff --git a/tox.ini b/tox.ini index 5a7d818e42..d537174f70 100644 --- a/tox.ini +++ b/tox.ini @@ -51,7 +51,6 @@ commands = black qiskit_aer test tools setup.py deps = -r requirements-dev.txt build - qiskit-ibmq-provider commands = python -I -m build --wheel -C=--build-option=-j4 pip install --find-links={toxinidir}/dist qiskit_aer