diff --git a/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt b/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt index 9540c4edb78..b26fae20dfd 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt +++ b/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt @@ -5,7 +5,7 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) -# Set OpenMp if it is available +# Set OpenMP if it is available find_package(OpenMP REQUIRED) if(OpenMP_CXX_FOUND) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") diff --git a/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h b/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h index 809cd8f16a5..f76437579e7 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h +++ b/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h @@ -1,846 +1,839 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + #pragma once + #include -#include "types.h" -#include "gates.h" #include #include #include -#include "quantum_state.hpp" -#include "basic_quantum_state.hpp" #include #include #include #include - #include - -#define MAX_QUBITS 1024 -#define MIN_QUBITS 64 +#include "quantum_state.hpp" +#include "basic_quantum_state.hpp" +#include "types.h" +#include "gates.h" using namespace std::literals::complex_literals; -namespace Microsoft -{ -namespace Quantum -{ -namespace SPARSESIMULATOR +namespace Microsoft::Quantum::SPARSESIMULATOR { + +constexpr logical_qubit_id MAX_QUBITS = 1024; +constexpr logical_qubit_id MIN_QUBITS = 64; + #ifndef M_PI #define M_PI 3.14159265358979323846 #endif - // Recrusively compiles sizes of QuantumState types between MIN_QUBITS and MAX_QUBITS // qubits large, growing by powers of 2 template std::shared_ptr construct_wfn_helper(logical_qubit_id nqubits) { - return (nqubits > max_num_bits / 2) ? - std::shared_ptr(new QuantumState()) - : (nqubits > MIN_QUBITS ? construct_wfn_helper(nqubits) : - std::shared_ptr(new QuantumState())); + return (nqubits > max_num_bits / 2) ? + std::shared_ptr(new QuantumState()) + : (nqubits > MIN_QUBITS ? construct_wfn_helper(nqubits) : + std::shared_ptr(new QuantumState())); } // Constructs a new quantum state, templated to use enough qubits to hold `nqubits`, // with the same state as `old_sim` template std::shared_ptr expand_wfn_helper(std::shared_ptr old_sim, logical_qubit_id nqubits) { - return (nqubits > max_num_bits / 2) ? std::shared_ptr(new QuantumState(old_sim)): expand_wfn_helper(old_sim, nqubits); + return (nqubits > max_num_bits / 2) ? std::shared_ptr(new QuantumState(old_sim)): expand_wfn_helper(old_sim, nqubits); } - - - - class SparseSimulator - { - public: - - std::set operations_done; +class SparseSimulator +{ +public: + + std::set operations_done; + + SparseSimulator(logical_qubit_id num_qubits) { + // Constructs a quantum state templated to the right number of qubits + // and returns a pointer to it as a basic_quantum_state + _quantum_state = construct_wfn_helper(num_qubits); + // Return the number of qubits this actually produces + num_qubits = _quantum_state->get_num_qubits(); + // Initialize with no qubits occupied + _occupied_qubits = std::vector(num_qubits, 0); + _max_num_qubits_used = 0; + _current_number_qubits_used = 0; + + _Ry_queue = std::vector(num_qubits, 0); + _Rx_queue = std::vector(num_qubits, 0); + _H_queue = std::vector(num_qubits, 0); + _Rx_angles = std::vector(num_qubits, 0.0); + _Ry_angles = std::vector(num_qubits, 0.0); + + } + + + ~SparseSimulator() { + _execute_queued_ops(); + } + + // Outputs the wavefunction to the console, after + // executing any queued operations + void DumpWavefunction(size_t indent = 0){ + _execute_queued_ops(); + _quantum_state->DumpWavefunction(indent); + } + + // Outputs the wavefunction as it is currently, + // without executing any operations + void DumpWavefunctionQuietly(size_t indent = 0){ + _quantum_state->DumpWavefunction(indent); + } + + void set_random_seed(unsigned seed = std::mt19937::default_seed){ + _quantum_state->set_random_seed(seed); + } + + // Returns the number of qubits currently available + // to the simulator, including those already used + logical_qubit_id get_num_qubits() { + return _quantum_state->get_num_qubits(); + } + + // Allocates a qubit at a specific location + // Implies that the caller of this function is tracking + // free qubits + void allocate_specific_qubit(logical_qubit_id qubit) { + size_t num_qubits = _quantum_state->get_num_qubits(); + // Checks that there are enough qubits + if (qubit >= num_qubits){ + // We create a new wavefunction and reallocate + std::shared_ptr old_state = _quantum_state; + _quantum_state = expand_wfn_helper(old_state, qubit+1); - SparseSimulator(logical_qubit_id num_qubits) { - // Constructs a quantum state templated to the right number of qubits - // and returns a pointer to it as a basic_quantum_state - _quantum_state = construct_wfn_helper(num_qubits); - // Return the number of qubits this actually produces num_qubits = _quantum_state->get_num_qubits(); - // Initialize with no qubits occupied - _occupied_qubits = std::vector(num_qubits, 0); - _max_num_qubits_used = 0; - _current_number_qubits_used = 0; - - _Ry_queue = std::vector(num_qubits, 0); - _Rx_queue = std::vector(num_qubits, 0); - _H_queue = std::vector(num_qubits, 0); - _Rx_angles = std::vector(num_qubits, 0.0); - _Ry_angles = std::vector(num_qubits, 0.0); - + _occupied_qubits.resize(num_qubits, 0); + _Ry_queue.resize(num_qubits, 0); + _Rx_queue.resize(num_qubits, 0); + _H_queue.resize(num_qubits, 0); + _Rx_angles.resize(num_qubits, 0.0); + _Ry_angles.resize(num_qubits, 0.0); + } + // The external qubit manager should prevent this, but this checks anyway + if (_occupied_qubits[qubit]) { + throw std::runtime_error("Qubit " + std::to_string(qubit) + " is already occupied"); } + // There is actually nothing to do to "allocate" a qubit, as every qubit + // is already available for use with this data structure + } + - ~SparseSimulator() { - _execute_queued_ops(); - } + // Removes a qubit in the zero state from the list + // of occupied qubits + bool release(logical_qubit_id qubit_id) { + // Quick check if it's zero + if (_occupied_qubits[qubit_id]) { + // If not zero here, we must execute any remaining operations + // Then check if the result is all zero + _execute_queued_ops(qubit_id); + + if (!_quantum_state->is_qubit_zero(qubit_id)){ + throw std::runtime_error("Released qubit not in zero state"); + } - // Outputs the wavefunction to the console, after - // executing any queued operations - void DumpWavefunction(size_t indent = 0){ - _execute_queued_ops(); - _quantum_state->DumpWavefunction(indent); } + _set_qubit_to_zero(qubit_id); + return true; + } - // Outputs the wavefunction as it is currently, - // without executing any operations - void DumpWavefunctionQuietly(size_t indent = 0){ - _quantum_state->DumpWavefunction(indent); - } - void set_random_seed(unsigned seed = std::mt19937::default_seed){ - _quantum_state->set_random_seed(seed); + void X(logical_qubit_id index) { + // XY = - YX + if (_Ry_queue[index]){ + _Ry_angles[index] *= -1.0; } - - // Returns the number of qubits currently available - // to the simulator, including those already used - logical_qubit_id get_num_qubits() { - return _quantum_state->get_num_qubits(); + // Rx trivially commutes + if (_H_queue[index]) { + _queued_operations.push_back(operation(OP::Z, index)); + return; } + _queued_operations.push_back(operation(OP::X, index)); + _set_qubit_to_nonzero(index); + } - // Allocates a qubit at a specific location - // Implies that the caller of this function is tracking - // free qubits - void allocate_specific_qubit(logical_qubit_id qubit) { - size_t num_qubits = _quantum_state->get_num_qubits(); - // Checks that there are enough qubits - if (qubit >= num_qubits){ - // We create a new wavefunction and reallocate - std::shared_ptr old_state = _quantum_state; - _quantum_state = expand_wfn_helper(old_state, qubit+1); - - num_qubits = _quantum_state->get_num_qubits(); - _occupied_qubits.resize(num_qubits, 0); - _Ry_queue.resize(num_qubits, 0); - _Rx_queue.resize(num_qubits, 0); - _H_queue.resize(num_qubits, 0); - _Rx_angles.resize(num_qubits, 0.0); - _Ry_angles.resize(num_qubits, 0.0); - } - // The external qubit manager should prevent this, but this checks anyway - if (_occupied_qubits[qubit]) { - throw std::runtime_error("Qubit " + std::to_string(qubit) + " is already occupied"); - } - // There is actually nothing to do to "allocate" a qubit, as every qubit - // is already available for use with this data structure - } - - - - // Removes a qubit in the zero state from the list - // of occupied qubits - bool release(logical_qubit_id qubit_id) { - // Quick check if it's zero - if (_occupied_qubits[qubit_id]) { - // If not zero here, we must execute any remaining operations - // Then check if the result is all zero - _execute_queued_ops(qubit_id); - - if (!_quantum_state->is_qubit_zero(qubit_id)){ - throw std::runtime_error("Released qubit not in zero state"); - } + // For both CNOT and all types of C*NOT + // If a control index is repeated, it just treats it as one control + // (Q# will throw an error in that condition) + void MCX(std::vector const& controls, logical_qubit_id target) { + // Check for anything on the controls + if (controls.size() > 1){ + _execute_if(controls); + } else { + // An H on the control but not the target forces execution + if (_Ry_queue[controls[0]] || _Rx_queue[controls[0]] || (_H_queue[controls[0]] && !_H_queue[target])){ + _execute_queued_ops(controls, OP::Ry); } - _set_qubit_to_zero(qubit_id); - return true; } - - - void X(logical_qubit_id index) { - // XY = - YX - if (_Ry_queue[index]){ - _Ry_angles[index] *= -1.0; - } - // Rx trivially commutes - if (_H_queue[index]) { - _queued_operations.push_back(operation(OP::Z, index)); - return; - } - _queued_operations.push_back(operation(OP::X, index)); - _set_qubit_to_nonzero(index); + // Ry on the target causes issues + if (_Ry_queue[target]){ + _execute_queued_ops(target, OP::Ry); } + // Rx on the target trivially commutes - - // For both CNOT and all types of C*NOT - // If a control index is repeated, it just treats it as one control - // (Q# will throw an error in that condition) - void MCX(std::vector const& controls, logical_qubit_id target) { - // Check for anything on the controls - if (controls.size() > 1){ - _execute_if(controls); + // An H on the target flips the operation + if (_H_queue[target]){ + // If it is a CNOT and there is also an H on the control, we swap control and target + if (controls.size() == 1 && _H_queue[controls[0]]){ + _queued_operations.push_back(operation(OP::MCX, controls[0], std::vector{target})); + _set_qubit_to_nonzero(controls[0]); } else { - // An H on the control but not the target forces execution - if (_Ry_queue[controls[0]] || _Rx_queue[controls[0]] || (_H_queue[controls[0]] && !_H_queue[target])){ - _execute_queued_ops(controls, OP::Ry); - } - } - // Ry on the target causes issues - if (_Ry_queue[target]){ - _execute_queued_ops(target, OP::Ry); - } - // Rx on the target trivially commutes - - // An H on the target flips the operation - if (_H_queue[target]){ - // If it is a CNOT and there is also an H on the control, we swap control and target - if (controls.size() == 1 && _H_queue[controls[0]]){ - _queued_operations.push_back(operation(OP::MCX, controls[0], std::vector{target})); - _set_qubit_to_nonzero(controls[0]); - } else { - _queued_operations.push_back(operation(OP::MCZ, target, controls)); - } - return; - } - // Queue the operation at this point - _queued_operations.push_back(operation(OP::MCX, target, controls)); - _set_qubit_to_nonzero(target); - } + _queued_operations.push_back(operation(OP::MCZ, target, controls)); + } + return; + } + // Queue the operation at this point + _queued_operations.push_back(operation(OP::MCX, target, controls)); + _set_qubit_to_nonzero(target); + } + + // Same as MCX, but we assert that the target is 0 before execution + void AND(std::vector const& controls, logical_qubit_id target) { + Assert(std::vector{Gates::Basis::PauliZ}, std::vector{target}, 0); + MCX(controls, target); + } + // Same as MCX, but we assert that the target is 0 after execution + void AdjAND(std::vector const& controls, logical_qubit_id target) { + MCX(controls, target); + Assert(std::vector{Gates::Basis::PauliZ}, std::vector{target}, 0); + _set_qubit_to_zero(target); + } + + void Y(logical_qubit_id index) { + // XY = -YX + if (_Rx_queue[index]){ + _Rx_angles[index] *= -1.0; + } + // commutes with H up to phase, so we ignore the H queue + _queued_operations.push_back(operation(OP::Y, index)); + _set_qubit_to_nonzero(index); + } + + void MCY(std::vector const& controls, logical_qubit_id target) { + _execute_if(controls); + // Commutes with Ry on the target, not Rx + if (_Rx_queue[target]){ + _execute_queued_ops(target, OP::Rx); + } + // YH = -YH, so we add a phase to track this + if (_H_queue[target]){ + // The phase added does not depend on the target + // Thus we use one of the controls as a target + _queued_operations.push_back(operation(OP::MCZ, controls[0], controls)); + } + _queued_operations.push_back(operation(OP::MCY, target, controls)); + _set_qubit_to_nonzero(target); + } - // Same as MCX, but we assert that the target is 0 before execution - void AND(std::vector const& controls, logical_qubit_id target) { - Assert(std::vector{Gates::Basis::PauliZ}, std::vector{target}, 0); - MCX(controls, target); + + void Z(logical_qubit_id index) { + // ZY = -YZ + if (_Ry_queue[index]){ + _Ry_angles[index] *= -1; } - // Same as MCX, but we assert that the target is 0 after execution - void AdjAND(std::vector const& controls, logical_qubit_id target) { - MCX(controls, target); - Assert(std::vector{Gates::Basis::PauliZ}, std::vector{target}, 0); - _set_qubit_to_zero(target); + // XZ = -ZX + if (_Rx_queue[index]){ + _Rx_angles[index] *= -1; } - - void Y(logical_qubit_id index) { - // XY = -YX - if (_Rx_queue[index]){ - _Rx_angles[index] *= -1.0; - } - // commutes with H up to phase, so we ignore the H queue - _queued_operations.push_back(operation(OP::Y, index)); + // HZ = XH + if (_H_queue[index]) { + _queued_operations.push_back(operation(OP::X, index)); _set_qubit_to_nonzero(index); + return; } + // No need to modified _occupied_qubits, since if a qubit is 0 + // a Z will not change that + _queued_operations.push_back(operation( OP::Z, index )); + } - void MCY(std::vector const& controls, logical_qubit_id target) { - _execute_if(controls); - // Commutes with Ry on the target, not Rx - if (_Rx_queue[target]){ - _execute_queued_ops(target, OP::Rx); + void MCZ(std::vector const& controls, logical_qubit_id target) { + // If the only thing on the controls is one H, we can switch + // this to an MCX. Any Rx or Ry, or more than 1 H, means we + // must execute. + size_t count = 0; + for (auto control : controls) { + if (_Ry_queue[control] || _Rx_queue[control]){ + count += 2; } - // YH = -YH, so we add a phase to track this - if (_H_queue[target]){ - // The phase added does not depend on the target - // Thus we use one of the controls as a target - _queued_operations.push_back(operation(OP::MCZ, controls[0], controls)); + if (_H_queue[control]){ + count++; } - _queued_operations.push_back(operation(OP::MCY, target, controls)); - _set_qubit_to_nonzero(target); } - - - void Z(logical_qubit_id index) { - // ZY = -YZ - if (_Ry_queue[index]){ - _Ry_angles[index] *= -1; - } - // XZ = -ZX - if (_Rx_queue[index]){ - _Rx_angles[index] *= -1; - } - // HZ = XH - if (_H_queue[index]) { - _queued_operations.push_back(operation(OP::X, index)); - _set_qubit_to_nonzero(index); - return; - } - // No need to modified _occupied_qubits, since if a qubit is 0 - // a Z will not change that - _queued_operations.push_back(operation( OP::Z, index )); - } - - void MCZ(std::vector const& controls, logical_qubit_id target) { - // If the only thing on the controls is one H, we can switch - // this to an MCX. Any Rx or Ry, or more than 1 H, means we - // must execute. - size_t count = 0; - for (auto control : controls) { - if (_Ry_queue[control] || _Rx_queue[control]){ - count += 2; - } + if (_Ry_queue[target] || _Rx_queue[target]){ + count +=2; + } + if (_H_queue[target]) {count++;} + if (count > 1) { + _execute_queued_ops(controls, OP::Ry); + _execute_queued_ops(target, OP::Ry); + } else if (count == 1){ + // Transform to an MCX, but we need to swap one of the controls + // with the target + std::vector new_controls(controls); + for (logical_qubit_id control : controls){ if (_H_queue[control]){ - count++; + std::swap(new_controls[control], target); } } - if (_Ry_queue[target] || _Rx_queue[target]){ - count +=2; - } - if (_H_queue[target]) {count++;} - if (count > 1) { - _execute_queued_ops(controls, OP::Ry); - _execute_queued_ops(target, OP::Ry); - } else if (count == 1){ - // Transform to an MCX, but we need to swap one of the controls - // with the target - std::vector new_controls(controls); - for (logical_qubit_id control : controls){ - if (_H_queue[control]){ - std::swap(new_controls[control], target); - } - } - _queued_operations.push_back(operation(OP::MCX, target, new_controls)); - _set_qubit_to_nonzero(target); - return; - } - _queued_operations.push_back(operation(OP::MCZ, target, controls)); - } - - - // Any phase gate - void Phase(amplitude const& phase, logical_qubit_id index) { - // Rx, Ry, and H do not commute well with arbitrary phase gates - if (_Ry_queue[index] || _Ry_queue[index] || _H_queue[index]){ - _execute_queued_ops(index, OP::Ry); - } - _queued_operations.push_back(operation(OP::Phase, index, phase)); - } - - void MCPhase(std::vector const& controls, amplitude const& phase, logical_qubit_id target){ - _execute_if(controls); - _execute_if(target); - _queued_operations.push_back(operation(OP::MCPhase, target, controls, phase)); + _queued_operations.push_back(operation(OP::MCX, target, new_controls)); + _set_qubit_to_nonzero(target); + return; } + _queued_operations.push_back(operation(OP::MCZ, target, controls)); + } - void T(logical_qubit_id index) { - Phase(amplitude(_normalizer_double, _normalizer_double), index); - } - void AdjT(logical_qubit_id index) { - Phase(amplitude(_normalizer_double, -_normalizer_double), index); + // Any phase gate + void Phase(amplitude const& phase, logical_qubit_id index) { + // Rx, Ry, and H do not commute well with arbitrary phase gates + if (_Ry_queue[index] || _Ry_queue[index] || _H_queue[index]){ + _execute_queued_ops(index, OP::Ry); } - + _queued_operations.push_back(operation(OP::Phase, index, phase)); + } - void R1(double const& angle, logical_qubit_id index) { - Phase(amplitude(std::cos(angle), std::sin(angle)), index); - } + void MCPhase(std::vector const& controls, amplitude const& phase, logical_qubit_id target){ + _execute_if(controls); + _execute_if(target); + _queued_operations.push_back(operation(OP::MCPhase, target, controls, phase)); + } - void MCR1(std::vector const& controls, double const& angle, logical_qubit_id target){ - MCPhase(controls, amplitude(std::cos(angle), std::sin(angle)), target); - } + void T(logical_qubit_id index) { + Phase(amplitude(_normalizer_double, _normalizer_double), index); + } - void R1Frac(std::int64_t numerator, std::int64_t power, logical_qubit_id index) { - R1((double)numerator * pow(0.5, power)*M_PI, index); - } + void AdjT(logical_qubit_id index) { + Phase(amplitude(_normalizer_double, -_normalizer_double), index); + } + - void MCR1Frac(std::vector const& controls, std::int64_t numerator, std::int64_t power, logical_qubit_id target){ - MCR1(controls, (double)numerator * pow(0.5, power) * M_PI, target); - } + void R1(double const& angle, logical_qubit_id index) { + Phase(amplitude(std::cos(angle), std::sin(angle)), index); + } - void S(logical_qubit_id index) { - Phase(1i, index); - } + void MCR1(std::vector const& controls, double const& angle, logical_qubit_id target){ + MCPhase(controls, amplitude(std::cos(angle), std::sin(angle)), target); + } - void AdjS(logical_qubit_id index) { - Phase(-1i, index); - } + void R1Frac(std::int64_t numerator, std::int64_t power, logical_qubit_id index) { + R1((double)numerator * pow(0.5, power)*M_PI, index); + } + void MCR1Frac(std::vector const& controls, std::int64_t numerator, std::int64_t power, logical_qubit_id target){ + MCR1(controls, (double)numerator * pow(0.5, power) * M_PI, target); + } + void S(logical_qubit_id index) { + Phase(1i, index); + } - void R(Gates::Basis b, double phi, logical_qubit_id index) - { - if (b == Gates::Basis::PauliI){ - return; - } + void AdjS(logical_qubit_id index) { + Phase(-1i, index); + } - // Tries to absorb the rotation into the existing queue, - // if it hits a different kind of rotation, the queue executes - if (b == Gates::Basis::PauliY){ - _Ry_queue[index] = true; - _Ry_angles[index] += phi; - _set_qubit_to_nonzero(index); - return; - } else if (_Ry_queue[index]) { - _execute_queued_ops(index, OP::Ry); - } - - if (b == Gates::Basis::PauliX){ - _Rx_queue[index] = true; - _Rx_angles[index] += phi; - _set_qubit_to_nonzero(index); - return; - } else if (_Rx_queue[index]){ - _execute_queued_ops(index, OP::Rz); - } - - // An Rz is just a phase - if (b == Gates::Basis::PauliZ){ - // HRz = RxH, but that's the wrong order for this structure - // Thus we must execute the H queue - if (_H_queue[index]){ - _execute_queued_ops(index, OP::H); - } - // Rz(phi) = RI(phi)*R1(-2*phi) - // Global phase from RI is ignored - R1(phi, index); - } - } - - void MCR (std::vector const& controls, Gates::Basis b, double phi, logical_qubit_id target) { - if (b == Gates::Basis::PauliI){ - // Controlled I rotations are equivalent to controlled phase gates - if (controls.size() > 1){ - MCPhase(controls, amplitude(std::cos(phi),std::sin(phi)), controls[0]); - } else { - Phase(amplitude(std::cos(phi),std::sin(phi)), controls[0]); - } - return; - } - _execute_if(controls); - // The target can commute with rotations of the same type - if (_Ry_queue[target] && b != Gates::Basis::PauliY){ - _execute_queued_ops(target, OP::Ry); - } - if (_Rx_queue[target] && b != Gates::Basis::PauliX){ - _execute_queued_ops(target, OP::Rx); - } - if (_H_queue[target]){ - _execute_queued_ops(target, OP::H); - } - // Execute any phase and permutation gates - // These are not indexed by qubit so it does - // not matter what the qubit argument is - _execute_queued_ops(0, OP::PermuteLarge); - _quantum_state->MCR(controls, b, phi, target); - _set_qubit_to_nonzero(target); - } - - void RFrac(Gates::Basis axis, std::int64_t numerator, std::int64_t power, logical_qubit_id index) { - // Opposite sign convention - R(axis, -(double)numerator * std::pow(0.5, power-1 )*M_PI, index); - } - - void MCRFrac(std::vector const& controls, Gates::Basis axis, std::int64_t numerator, std::int64_t power, logical_qubit_id target) { - // Opposite sign convention - MCR(controls, axis, -(double)numerator * pow(0.5, power - 1) * M_PI, target); - } - - void Exp(std::vector const& axes, double angle, std::vector const& qubits){ - amplitude cosAngle = std::cos(angle); - amplitude sinAngle = 1i*std::sin(angle); - // This does not commute nicely with anything, so we execute everything - _execute_queued_ops(qubits); - _quantum_state->PauliCombination(axes, qubits, cosAngle, sinAngle); - for (auto qubit : qubits){ - _set_qubit_to_nonzero(qubit); - } - } - void MCExp(std::vector const& controls, std::vector const& axes, double angle, std::vector const& qubits){ - amplitude cosAngle = std::cos(angle); - amplitude sinAngle = 1i*std::sin(angle); - // This does not commute nicely with anything, so we execute everything - _execute_queued_ops(qubits); - _execute_queued_ops(controls); - _quantum_state->MCPauliCombination(controls, axes, qubits, cosAngle, sinAngle); - for (auto qubit : qubits){ - _set_qubit_to_nonzero(qubit); - } + void R(Gates::Basis b, double phi, logical_qubit_id index) + { + if (b == Gates::Basis::PauliI){ + return; } - - - void H(logical_qubit_id index) { - // YH = -HY - _Ry_angles[index] *= (_Ry_queue[index] ? -1.0 : 1.0); - // Commuting with Rx creates a phase, but on the wrong side - // So we execute any Rx immediately - if (_Rx_queue[index]){ - _execute_queued_ops(index, OP::Rx); - } - _H_queue[index] = !_H_queue[index]; + // Tries to absorb the rotation into the existing queue, + // if it hits a different kind of rotation, the queue executes + if (b == Gates::Basis::PauliY){ + _Ry_queue[index] = true; + _Ry_angles[index] += phi; _set_qubit_to_nonzero(index); + return; + } else if (_Ry_queue[index]) { + _execute_queued_ops(index, OP::Ry); } - void MCH(std::vector const& controls, logical_qubit_id target) { - // No commutation on controls - _execute_if(controls); - // No Ry or Rx commutation on target - if (_Ry_queue[target] || _Rx_queue[target]){ - _execute_queued_ops(target, OP::Ry); - } - // Commutes through H gates on the target, so it does not check - _execute_phase_and_permute(); - _quantum_state->MCH(controls, target); - _set_qubit_to_nonzero(target); + if (b == Gates::Basis::PauliX){ + _Rx_queue[index] = true; + _Rx_angles[index] += phi; + _set_qubit_to_nonzero(index); + return; + } else if (_Rx_queue[index]){ + _execute_queued_ops(index, OP::Rz); } - - - - void SWAP(logical_qubit_id index_1, logical_qubit_id index_2){ - // This is necessary for the "shift" to make sense - if (index_1 > index_2){ - std::swap(index_2, index_1); - } - // Everything commutes nicely with a swap - _Ry_queue.swap(_Ry_queue[index_1], _Ry_queue[index_2]); - std::swap(_Ry_angles[index_1], _Ry_angles[index_2]); - _Rx_queue.swap(_Rx_queue[index_1], _Rx_queue[index_2]); - std::swap(_Rx_angles[index_1], _Rx_angles[index_2]); - _H_queue.swap(_H_queue[index_1], _H_queue[index_2]); - _occupied_qubits.swap(_occupied_qubits[index_1], _occupied_qubits[index_2]); - logical_qubit_id shift = index_2 - index_1; - _queued_operations.push_back(operation(OP::SWAP, index_1, shift, index_2)); - } - - void CSWAP(std::vector const& controls, logical_qubit_id index_1, logical_qubit_id index_2){ - if (index_1 > index_2){ - std::swap(index_2, index_1); - } - // Nothing commutes nicely with a controlled swap - _execute_if(controls); - _execute_if(index_1); - _execute_if(index_2); - - logical_qubit_id shift = index_2 - index_1; - _queued_operations.push_back(operation(OP::MCSWAP, index_1, shift, controls, index_2)); - // If either qubit is occupied, then set them both to occupied - if(_occupied_qubits[index_1] || _occupied_qubits[index_2]){ - _set_qubit_to_nonzero(index_1); - _set_qubit_to_nonzero(index_2); + // An Rz is just a phase + if (b == Gates::Basis::PauliZ){ + // HRz = RxH, but that's the wrong order for this structure + // Thus we must execute the H queue + if (_H_queue[index]){ + _execute_queued_ops(index, OP::H); } + // Rz(phi) = RI(phi)*R1(-2*phi) + // Global phase from RI is ignored + R1(phi, index); } + } - bool M(logical_qubit_id target) { - // Do nothing if the qubit is known to be 0 - if (!_occupied_qubits[target]){ - return false; - } - // If we get a measurement, we take it as soon as we can - _execute_queued_ops(target, OP::Ry); - // If we measure 0, then this resets the occupied qubit register - if (_quantum_state->M(target)){ - _set_qubit_to_nonzero(target); + void MCR (std::vector const& controls, Gates::Basis b, double phi, logical_qubit_id target) { + if (b == Gates::Basis::PauliI){ + // Controlled I rotations are equivalent to controlled phase gates + if (controls.size() > 1){ + MCPhase(controls, amplitude(std::cos(phi),std::sin(phi)), controls[0]); } else { - _set_qubit_to_zero(target); - } - return _occupied_qubits[target]; - - } - - void Reset(logical_qubit_id target) { - if (!_occupied_qubits[target]){ return; } - // If we get a measurement, we take it as soon as we can - _execute_queued_ops(target, OP::Ry); - _quantum_state->Reset(target); - _set_qubit_to_zero(target); - } - - void Assert(std::vector axes, std::vector const& qubits, bool result) { - // Assertions will not commute well with Rx or Ry - for (auto qubit : qubits) { - if (_Rx_queue[qubit] || _Ry_queue[qubit]){ - _execute_queued_ops(qubits, OP::Ry); - } - } - bool isAllZ = true; - bool isEmpty = true; - // Process each assertion by H commutation - for (int i = 0; i < qubits.size(); i++) { - switch (axes[i]){ - case Gates::Basis::PauliY: - // HY=YH, so we switch the eigenvalue - if (_H_queue[qubits[i]]){ - result ^= _H_queue[qubits[i]]; - } - isAllZ = false; - isEmpty = false; - break; - case Gates::Basis::PauliX: - // HX = ZH - if (_H_queue[qubits[i]]){ - axes[i] = Gates::Basis::PauliZ; - } else { - isAllZ = false; - } - isEmpty = false; - break; - case Gates::Basis::PauliZ: - // HZ = XH - if (_H_queue[qubits[i]]){ - axes[i] = Gates::Basis::PauliX; - isAllZ = false; - } - isEmpty = false; - break; - default: - break; - } - } - if (isEmpty) { - return; - } - // Z assertions are like phase gates - // If it's in release mode, it will queue them - // as a phase/permutation gate - // This means if an assert fails, it will fail - // at some future point, not at the point of failure - #if NDEBUG - if (isAllZ) { - _queued_operations.push_back(operation(OP::Assert, qubits, result)); - return; - } - #endif - // X or Y assertions require execution - _execute_queued_ops(qubits, OP::PermuteLarge); - _quantum_state->Assert(axes, qubits, result); - } - - // Returns the probability of a given measurement in a Pauli basis - // by decomposing each pair of computational basis states into eigenvectors - // and adding the coefficients of the respective components - double MeasurementProbability(std::vector const& axes, std::vector const& qubits) { - _execute_queued_ops(qubits, OP::Ry); - return _quantum_state->MeasurementProbability(axes, qubits); - } - - - - bool Measure(std::vector const& axes, std::vector const& qubits){ - _execute_queued_ops(qubits, OP::Ry); - bool result = _quantum_state->Measure(axes, qubits); - // Switch basis to save space - // Idea being that, e.g., HH = I, but if we know - // that the qubit is in the X-basis, we can apply H - // and execute, and this will send that qubit to all ones - // or all zeros; then we leave the second H in the queue - // Ideally we would also do that with Y, but HS would force execution, - // rendering it pointless - std::vector measurements; - for (int i =0; i < axes.size(); i++){ - if (axes[i]==Gates::Basis::PauliX){ - H(qubits[i]); - measurements.push_back(qubits[i]); - } - } - _execute_queued_ops(measurements, OP::H); - // These operations undo the previous operations, but they will be - // queued - for (int i =0; i < axes.size(); i++){ - if (axes[i]==Gates::Basis::PauliX){ - H(qubits[i]); - } + Phase(amplitude(std::cos(phi),std::sin(phi)), controls[0]); } - return result; + return; } - // Returns the amplitude of a given bitstring - amplitude probe(std::string label) { - _execute_queued_ops(); - return _quantum_state->probe(label); - } - - std::string Sample() { - _execute_queued_ops(); - return _quantum_state->Sample(); - } - - // Dumps the state of a subspace of particular qubits, if they are not entangled - // This requires it to detect if the subspace is entangled, construct a new - // projected wavefunction, then call the `callback` function on each state. - bool dump_qubits(std::vector qubits, void (*callback)(char*, double, double)) { - _execute_queued_ops(qubits, OP::Ry); - return _quantum_state->dump_qubits(qubits, callback); - } - - // Dumps all the states in superposition via a callback function - void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) { - _execute_queued_ops(); - _quantum_state->dump_all(max_qubit_id, callback); - } - - // Updates state to all queued gates - void update_state() { - _execute_queued_ops(); - } - - - private: - - // These indicate whether there are any H, Rx, or Ry gates - // that have yet to be applied to the wavefunction. - // Since HH=I and Rx(theta_1)Rx(theta_2) = Rx(theta_1+theta_2) - // it only needs a boolean to track them. - std::vector _H_queue; - std::vector _Rx_queue; - std::vector _Ry_queue; + _execute_if(controls); + // The target can commute with rotations of the same type + if (_Ry_queue[target] && b != Gates::Basis::PauliY){ + _execute_queued_ops(target, OP::Ry); + } + if (_Rx_queue[target] && b != Gates::Basis::PauliX){ + _execute_queued_ops(target, OP::Rx); + } + if (_H_queue[target]){ + _execute_queued_ops(target, OP::H); + } + // Execute any phase and permutation gates + // These are not indexed by qubit so it does + // not matter what the qubit argument is + _execute_queued_ops(0, OP::PermuteLarge); + _quantum_state->MCR(controls, b, phi, target); + _set_qubit_to_nonzero(target); + } + + void RFrac(Gates::Basis axis, std::int64_t numerator, std::int64_t power, logical_qubit_id index) { + // Opposite sign convention + R(axis, -(double)numerator * std::pow(0.5, power-1 )*M_PI, index); + } + + void MCRFrac(std::vector const& controls, Gates::Basis axis, std::int64_t numerator, std::int64_t power, logical_qubit_id target) { + // Opposite sign convention + MCR(controls, axis, -(double)numerator * pow(0.5, power - 1) * M_PI, target); + } + + void Exp(std::vector const& axes, double angle, std::vector const& qubits){ + amplitude cosAngle = std::cos(angle); + amplitude sinAngle = 1i*std::sin(angle); + // This does not commute nicely with anything, so we execute everything + _execute_queued_ops(qubits); + _quantum_state->PauliCombination(axes, qubits, cosAngle, sinAngle); + for (auto qubit : qubits){ + _set_qubit_to_nonzero(qubit); + } + } + + void MCExp(std::vector const& controls, std::vector const& axes, double angle, std::vector const& qubits){ + amplitude cosAngle = std::cos(angle); + amplitude sinAngle = 1i*std::sin(angle); + // This does not commute nicely with anything, so we execute everything + _execute_queued_ops(qubits); + _execute_queued_ops(controls); + _quantum_state->MCPauliCombination(controls, axes, qubits, cosAngle, sinAngle); + for (auto qubit : qubits){ + _set_qubit_to_nonzero(qubit); + } + } - std::vector _Rx_angles; - std::vector _Ry_angles; + + + void H(logical_qubit_id index) { + // YH = -HY + _Ry_angles[index] *= (_Ry_queue[index] ? -1.0 : 1.0); + // Commuting with Rx creates a phase, but on the wrong side + // So we execute any Rx immediately + if (_Rx_queue[index]){ + _execute_queued_ops(index, OP::Rx); + } + _H_queue[index] = !_H_queue[index]; + _set_qubit_to_nonzero(index); + } + + void MCH(std::vector const& controls, logical_qubit_id target) { + // No commutation on controls + _execute_if(controls); + // No Ry or Rx commutation on target + if (_Ry_queue[target] || _Rx_queue[target]){ + _execute_queued_ops(target, OP::Ry); + } + // Commutes through H gates on the target, so it does not check + _execute_phase_and_permute(); + _quantum_state->MCH(controls, target); + _set_qubit_to_nonzero(target); + } - // Store which qubits are non-zero as a bitstring - std::vector _occupied_qubits; - logical_qubit_id _max_num_qubits_used = 0; - logical_qubit_id _current_number_qubits_used; - // In a situation where we know a qubit is zero, - // this sets the occupied qubit vector and decrements - // the current number of qubits if necessary - void _set_qubit_to_zero(logical_qubit_id index){ - if (_occupied_qubits[index]){ - --_current_number_qubits_used; - } - _occupied_qubits[index] = false; - } + - // In a situation where a qubit may be non-zero, - // we increment which qubits are used, and update the current - // and maximum number of qubits - void _set_qubit_to_nonzero(logical_qubit_id index){ - if (!_occupied_qubits[index]){ - ++_current_number_qubits_used; - _max_num_qubits_used = std::max(_max_num_qubits_used, _current_number_qubits_used); - } - _occupied_qubits[index] = true; + void SWAP(logical_qubit_id index_1, logical_qubit_id index_2){ + // This is necessary for the "shift" to make sense + if (index_1 > index_2){ + std::swap(index_2, index_1); + } + // Everything commutes nicely with a swap + _Ry_queue.swap(_Ry_queue[index_1], _Ry_queue[index_2]); + std::swap(_Ry_angles[index_1], _Ry_angles[index_2]); + _Rx_queue.swap(_Rx_queue[index_1], _Rx_queue[index_2]); + std::swap(_Rx_angles[index_1], _Rx_angles[index_2]); + _H_queue.swap(_H_queue[index_1], _H_queue[index_2]); + _occupied_qubits.swap(_occupied_qubits[index_1], _occupied_qubits[index_2]); + logical_qubit_id shift = index_2 - index_1; + _queued_operations.push_back(operation(OP::SWAP, index_1, shift, index_2)); + } + + void CSWAP(std::vector const& controls, logical_qubit_id index_1, logical_qubit_id index_2){ + if (index_1 > index_2){ + std::swap(index_2, index_1); + } + // Nothing commutes nicely with a controlled swap + _execute_if(controls); + _execute_if(index_1); + _execute_if(index_2); + + logical_qubit_id shift = index_2 - index_1; + _queued_operations.push_back(operation(OP::MCSWAP, index_1, shift, controls, index_2)); + // If either qubit is occupied, then set them both to occupied + if(_occupied_qubits[index_1] || _occupied_qubits[index_2]){ + _set_qubit_to_nonzero(index_1); + _set_qubit_to_nonzero(index_2); + } + } + + bool M(logical_qubit_id target) { + // Do nothing if the qubit is known to be 0 + if (!_occupied_qubits[target]){ + return false; + } + // If we get a measurement, we take it as soon as we can + _execute_queued_ops(target, OP::Ry); + // If we measure 0, then this resets the occupied qubit register + if (_quantum_state->M(target)){ + _set_qubit_to_nonzero(target); + } else { + _set_qubit_to_zero(target); } + return _occupied_qubits[target]; - // Normalizer for T gates: 1/sqrt(2) - const double _normalizer_double = 1.0 / std::sqrt(2.0); - - // Internal quantum state - std::shared_ptr _quantum_state; - - // Queued phase and permutation operations - std::list _queued_operations; + } - // The next three functions execute the H, and/or Rx, and/or Ry - // queues on a single qubit - void _execute_RyRxH_single_qubit(logical_qubit_id const &index){ - if (_H_queue[index]){ - _quantum_state->H(index); - _H_queue[index] = false; - } - if (_Rx_queue[index]){ - _quantum_state->R(Gates::Basis::PauliX, _Rx_angles[index], index); - _Rx_angles[index] = 0.0; - _Rx_queue[index] = false; - } - if (_Ry_queue[index]){ - _quantum_state->R(Gates::Basis::PauliY, _Ry_angles[index], index); - _Ry_angles[index] = 0.0; - _Ry_queue[index] = false; + void Reset(logical_qubit_id target) { + if (!_occupied_qubits[target]){ return; } + // If we get a measurement, we take it as soon as we can + _execute_queued_ops(target, OP::Ry); + _quantum_state->Reset(target); + _set_qubit_to_zero(target); + } + + void Assert(std::vector axes, std::vector const& qubits, bool result) { + // Assertions will not commute well with Rx or Ry + for (auto qubit : qubits) { + if (_Rx_queue[qubit] || _Ry_queue[qubit]){ + _execute_queued_ops(qubits, OP::Ry); + } + } + bool isAllZ = true; + bool isEmpty = true; + // Process each assertion by H commutation + for (int i = 0; i < qubits.size(); i++) { + switch (axes[i]){ + case Gates::Basis::PauliY: + // HY=YH, so we switch the eigenvalue + if (_H_queue[qubits[i]]){ + result ^= _H_queue[qubits[i]]; + } + isAllZ = false; + isEmpty = false; + break; + case Gates::Basis::PauliX: + // HX = ZH + if (_H_queue[qubits[i]]){ + axes[i] = Gates::Basis::PauliZ; + } else { + isAllZ = false; + } + isEmpty = false; + break; + case Gates::Basis::PauliZ: + // HZ = XH + if (_H_queue[qubits[i]]){ + axes[i] = Gates::Basis::PauliX; + isAllZ = false; + } + isEmpty = false; + break; + default: + break; } } - - void _execute_RxH_single_qubit(logical_qubit_id const &index){ - if (_H_queue[index]){ - _quantum_state->H(index); - _H_queue[index] = false; - } - if (_Rx_queue[index]){ - _quantum_state->R(Gates::Basis::PauliX, _Rx_angles[index], index); - _Rx_angles[index] = 0.0; - _Rx_queue[index] = false; - } + if (isEmpty) { + return; } - - void _execute_H_single_qubit(logical_qubit_id const &index){ - if (_H_queue[index]){ - _quantum_state->H(index); - _H_queue[index] = false; + // Z assertions are like phase gates + // If it's in release mode, it will queue them + // as a phase/permutation gate + // This means if an assert fails, it will fail + // at some future point, not at the point of failure + #if NDEBUG + if (isAllZ) { + _queued_operations.push_back(operation(OP::Assert, qubits, result)); + return; } + #endif + // X or Y assertions require execution + _execute_queued_ops(qubits, OP::PermuteLarge); + _quantum_state->Assert(axes, qubits, result); + } + + // Returns the probability of a given measurement in a Pauli basis + // by decomposing each pair of computational basis states into eigenvectors + // and adding the coefficients of the respective components + double MeasurementProbability(std::vector const& axes, std::vector const& qubits) { + _execute_queued_ops(qubits, OP::Ry); + return _quantum_state->MeasurementProbability(axes, qubits); + } + + + + bool Measure(std::vector const& axes, std::vector const& qubits){ + _execute_queued_ops(qubits, OP::Ry); + bool result = _quantum_state->Measure(axes, qubits); + // Switch basis to save space + // Idea being that, e.g., HH = I, but if we know + // that the qubit is in the X-basis, we can apply H + // and execute, and this will send that qubit to all ones + // or all zeros; then we leave the second H in the queue + // Ideally we would also do that with Y, but HS would force execution, + // rendering it pointless + std::vector measurements; + for (int i =0; i < axes.size(); i++){ + if (axes[i]==Gates::Basis::PauliX){ + H(qubits[i]); + measurements.push_back(qubits[i]); + } } - - // Executes all phase and permutation operations, if any exist - void _execute_phase_and_permute(){ - if (_queued_operations.size() != 0){ - _quantum_state->phase_and_permute(_queued_operations); - _queued_operations.clear(); - } + _execute_queued_ops(measurements, OP::H); + // These operations undo the previous operations, but they will be + // queued + for (int i =0; i < axes.size(); i++){ + if (axes[i]==Gates::Basis::PauliX){ + H(qubits[i]); + } } - - // Executes all queued operations (including H and rotations) - // on all qubits - void _execute_queued_ops() { - _execute_phase_and_permute(); - logical_qubit_id num_qubits = _quantum_state->get_num_qubits(); - for (logical_qubit_id index =0; index < num_qubits; index++){ + return result; + } + + // Returns the amplitude of a given bitstring + amplitude probe(std::string label) { + _execute_queued_ops(); + return _quantum_state->probe(label); + } + + std::string Sample() { + _execute_queued_ops(); + return _quantum_state->Sample(); + } + + // Dumps the state of a subspace of particular qubits, if they are not entangled + // This requires it to detect if the subspace is entangled, construct a new + // projected wavefunction, then call the `callback` function on each state. + bool dump_qubits(std::vector qubits, void (*callback)(char*, double, double)) { + _execute_queued_ops(qubits, OP::Ry); + return _quantum_state->dump_qubits(qubits, callback); + } + + // Dumps all the states in superposition via a callback function + void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) { + _execute_queued_ops(); + _quantum_state->dump_all(max_qubit_id, callback); + } + + // Updates state to all queued gates + void update_state() { + _execute_queued_ops(); + } + + +private: + + // These indicate whether there are any H, Rx, or Ry gates + // that have yet to be applied to the wavefunction. + // Since HH=I and Rx(theta_1)Rx(theta_2) = Rx(theta_1+theta_2) + // it only needs a boolean to track them. + std::vector _H_queue; + std::vector _Rx_queue; + std::vector _Ry_queue; + + std::vector _Rx_angles; + std::vector _Ry_angles; + + // Store which qubits are non-zero as a bitstring + std::vector _occupied_qubits; + logical_qubit_id _max_num_qubits_used = 0; + logical_qubit_id _current_number_qubits_used; + + // In a situation where we know a qubit is zero, + // this sets the occupied qubit vector and decrements + // the current number of qubits if necessary + void _set_qubit_to_zero(logical_qubit_id index){ + if (_occupied_qubits[index]){ + --_current_number_qubits_used; + } + _occupied_qubits[index] = false; + } + + // In a situation where a qubit may be non-zero, + // we increment which qubits are used, and update the current + // and maximum number of qubits + void _set_qubit_to_nonzero(logical_qubit_id index){ + if (!_occupied_qubits[index]){ + ++_current_number_qubits_used; + _max_num_qubits_used = std::max(_max_num_qubits_used, _current_number_qubits_used); + } + _occupied_qubits[index] = true; + } + + // Normalizer for T gates: 1/sqrt(2) + const double _normalizer_double = 1.0 / std::sqrt(2.0); + + // Internal quantum state + std::shared_ptr _quantum_state; + + // Queued phase and permutation operations + std::list _queued_operations; + + // The next three functions execute the H, and/or Rx, and/or Ry + // queues on a single qubit + void _execute_RyRxH_single_qubit(logical_qubit_id const &index){ + if (_H_queue[index]){ + _quantum_state->H(index); + _H_queue[index] = false; + } + if (_Rx_queue[index]){ + _quantum_state->R(Gates::Basis::PauliX, _Rx_angles[index], index); + _Rx_angles[index] = 0.0; + _Rx_queue[index] = false; + } + if (_Ry_queue[index]){ + _quantum_state->R(Gates::Basis::PauliY, _Ry_angles[index], index); + _Ry_angles[index] = 0.0; + _Ry_queue[index] = false; + } + } + + void _execute_RxH_single_qubit(logical_qubit_id const &index){ + if (_H_queue[index]){ + _quantum_state->H(index); + _H_queue[index] = false; + } + if (_Rx_queue[index]){ + _quantum_state->R(Gates::Basis::PauliX, _Rx_angles[index], index); + _Rx_angles[index] = 0.0; + _Rx_queue[index] = false; + } + } + + void _execute_H_single_qubit(logical_qubit_id const &index){ + if (_H_queue[index]){ + _quantum_state->H(index); + _H_queue[index] = false; + } + } + + // Executes all phase and permutation operations, if any exist + void _execute_phase_and_permute(){ + if (_queued_operations.size() != 0){ + _quantum_state->phase_and_permute(_queued_operations); + _queued_operations.clear(); + } + } + + // Executes all queued operations (including H and rotations) + // on all qubits + void _execute_queued_ops() { + _execute_phase_and_permute(); + logical_qubit_id num_qubits = _quantum_state->get_num_qubits(); + for (logical_qubit_id index =0; index < num_qubits; index++){ + _execute_RyRxH_single_qubit(index); + } + } + + // Executes all phase and permutation operations, + // then any H, Rx, or Ry gates queued on the qubit index, + // up to the level specified (where H < Rx < Ry) + void _execute_queued_ops(logical_qubit_id index, OP level = OP::Ry){ + _execute_phase_and_permute(); + switch (level){ + case OP::Ry: _execute_RyRxH_single_qubit(index); - } - } - - // Executes all phase and permutation operations, - // then any H, Rx, or Ry gates queued on the qubit index, - // up to the level specified (where H < Rx < Ry) - void _execute_queued_ops(logical_qubit_id index, OP level = OP::Ry){ - _execute_phase_and_permute(); - switch (level){ - case OP::Ry: + break; + case OP::Rx: + _execute_RxH_single_qubit(index); + break; + case OP::H: + _execute_H_single_qubit(index); + break; + default: + break; + } + } + + // Executes all phase and permutation operations, + // then any H, Rx, or Ry gates queued on any of the qubit indices, + // up to the level specified (where H < Rx < Ry) + void _execute_queued_ops(std::vector indices, OP level = OP::Ry){ + _execute_phase_and_permute(); + switch (level){ + case OP::Ry: + for (auto index : indices){ _execute_RyRxH_single_qubit(index); - break; - case OP::Rx: + } + break; + case OP::Rx: + for (auto index : indices){ _execute_RxH_single_qubit(index); - break; - case OP::H: + } + break; + case OP::H: + for (auto index : indices){ _execute_H_single_qubit(index); - break; - default: - break; - } - } - - // Executes all phase and permutation operations, - // then any H, Rx, or Ry gates queued on any of the qubit indices, - // up to the level specified (where H < Rx < Ry) - void _execute_queued_ops(std::vector indices, OP level = OP::Ry){ - _execute_phase_and_permute(); - switch (level){ - case OP::Ry: - for (auto index : indices){ - _execute_RyRxH_single_qubit(index); - } - break; - case OP::Rx: - for (auto index : indices){ - _execute_RxH_single_qubit(index); - } - break; - case OP::H: - for (auto index : indices){ - _execute_H_single_qubit(index); - } - break; - default: - break; - } + } + break; + default: + break; } + } - // Executes if there is anything already queued on the qubit target - // Used when queuing gates that do not commute well - void _execute_if(logical_qubit_id &target){ - if (_Ry_queue[target] || _Rx_queue[target] || _H_queue[target]){ - _execute_queued_ops(target, OP::Ry); - } + // Executes if there is anything already queued on the qubit target + // Used when queuing gates that do not commute well + void _execute_if(logical_qubit_id &target){ + if (_Ry_queue[target] || _Rx_queue[target] || _H_queue[target]){ + _execute_queued_ops(target, OP::Ry); } + } - // Executes if there is anything already queued on the qubits in controls - // Used when queuing gates that do not commute well - void _execute_if(std::vector const &controls) { - for (auto control : controls){ - if (_Ry_queue[control] || _Rx_queue[control] || _H_queue[control]){ - _execute_queued_ops(controls, OP::Ry); - return; - } + // Executes if there is anything already queued on the qubits in controls + // Used when queuing gates that do not commute well + void _execute_if(std::vector const &controls) { + for (auto control : controls){ + if (_Ry_queue[control] || _Rx_queue[control] || _H_queue[control]){ + _execute_queued_ops(controls, OP::Ry); + return; } } + } - }; - +}; -} // namespace SPARSESIMULATOR -} // namespace Quantum -} // namespace Microsoft +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/basic_quantum_state.hpp b/src/Simulation/Simulators/SparseSimulator/Native/basic_quantum_state.hpp index 330261881f5..3f67d6a5211 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/basic_quantum_state.hpp +++ b/src/Simulation/Simulators/SparseSimulator/Native/basic_quantum_state.hpp @@ -1,82 +1,77 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + #pragma once -#include "types.h" -#include "gates.h" + #include +#include "types.h" +#include "gates.h" -namespace Microsoft +namespace Microsoft::Quantum::SPARSESIMULATOR { -namespace Quantum -{ -namespace SPARSESIMULATOR -{ - - - // Virtual class for QuantumState - // This is not templated, so it allows SparseSimulator types to avoid templates - class BasicQuantumState - { - public: - - BasicQuantumState() {} +// Virtual class for QuantumState +// This is not templated, so it allows SparseSimulator types to avoid templates +class BasicQuantumState +{ +public: - virtual logical_qubit_id get_num_qubits() = 0; + BasicQuantumState() {} - virtual void DumpWavefunction(size_t indent = 0) = 0; + virtual logical_qubit_id get_num_qubits() = 0; - virtual void set_random_seed(unsigned seed = std::mt19937::default_seed) = 0; + virtual void DumpWavefunction(size_t indent = 0) = 0; - virtual void set_precision(double new_precision) = 0; + virtual void set_random_seed(unsigned seed = std::mt19937::default_seed) = 0; - virtual float get_load_factor() = 0; + virtual void set_precision(double new_precision) = 0; - virtual void set_load_factor(float new_load_factor) = 0; + virtual float get_load_factor() = 0; - virtual size_t get_wavefunction_size() = 0; + virtual void set_load_factor(float new_load_factor) = 0; - virtual void PauliCombination(std::vector const&, std::vector const&, amplitude, amplitude) = 0; - virtual void MCPauliCombination(std::vector const&, std::vector const&, std::vector const&, amplitude, amplitude) = 0; + virtual size_t get_wavefunction_size() = 0; - virtual bool M(logical_qubit_id) = 0; + virtual void PauliCombination(std::vector const&, std::vector const&, amplitude, amplitude) = 0; + virtual void MCPauliCombination(std::vector const&, std::vector const&, std::vector const&, amplitude, amplitude) = 0; - virtual void Reset(logical_qubit_id) = 0; + virtual bool M(logical_qubit_id) = 0; - + virtual void Reset(logical_qubit_id) = 0; - virtual void Assert(std::vector const&, std::vector const&, bool) = 0; + - virtual double MeasurementProbability(std::vector const&, std::vector const&) = 0; - virtual bool Measure(std::vector const&, std::vector const&) = 0; + virtual void Assert(std::vector const&, std::vector const&, bool) = 0; + virtual double MeasurementProbability(std::vector const&, std::vector const&) = 0; + virtual bool Measure(std::vector const&, std::vector const&) = 0; - virtual amplitude probe(std::string label) = 0; - virtual bool dump_qubits(std::vector qubits, void (*callback)(char*, double, double)) = 0; + virtual amplitude probe(std::string label) = 0; - virtual void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) = 0; + virtual bool dump_qubits(std::vector qubits, void (*callback)(char*, double, double)) = 0; - virtual void phase_and_permute(std::listconst &) = 0; + virtual void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) = 0; - virtual void R(Gates::Basis b, double phi, logical_qubit_id index) = 0; - virtual void MCR (std::vector const&, Gates::Basis, double, logical_qubit_id) = 0; + virtual void phase_and_permute(std::listconst &) = 0; - virtual void H(logical_qubit_id index) = 0; - virtual void MCH(std::vector const& controls, logical_qubit_id index) = 0; + virtual void R(Gates::Basis b, double phi, logical_qubit_id index) = 0; + virtual void MCR (std::vector const&, Gates::Basis, double, logical_qubit_id) = 0; - virtual bool is_qubit_zero(logical_qubit_id) = 0; + virtual void H(logical_qubit_id index) = 0; + virtual void MCH(std::vector const& controls, logical_qubit_id index) = 0; - virtual universal_wavefunction get_universal_wavefunction() = 0; + virtual bool is_qubit_zero(logical_qubit_id) = 0; - virtual std::function get_rng() = 0; + virtual universal_wavefunction get_universal_wavefunction() = 0; - virtual void complete_threads() = 0; - virtual int get_num_threads() = 0; + virtual std::function get_rng() = 0; - virtual std::string Sample() = 0; - }; + virtual void complete_threads() = 0; + virtual int get_num_threads() = 0; + virtual std::string Sample() = 0; +}; -} // namespace SPARSESIMULATOR -} // namespace Quantum -} // namespace Microsoft +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp b/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp index 24bd58358d2..999702d665b 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp +++ b/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp @@ -1,4 +1,4 @@ -// Copyright (c) Microsoft Corporation. All rights reserved. +// Copyright (c) Microsoft Corporation. // Licensed under the MIT License. // Wrapper functions for basic C++ functions @@ -6,18 +6,17 @@ // index into the vector of simulators, // then call a member function +#include +#include +#include +#include + #include "capi.hpp" #include "SparseSimulator.h" #include "factory.hpp" -#include -#include using namespace Microsoft::Quantum::SPARSESIMULATOR; - -#include -#include - std::string sample_string; extern "C" diff --git a/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp b/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp index 85add5e5b84..8f4c3839041 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp +++ b/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp @@ -1,11 +1,10 @@ -// Copyright (c) Microsoft Corporation. All rights reserved. +// Copyright (c) Microsoft Corporation. // Licensed under the MIT License. #pragma once - - #include + #include "types.h" // SAL only defined in windows. diff --git a/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp b/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp index fa42fd3c5fe..e8aa09a7a22 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp +++ b/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp @@ -1,21 +1,18 @@ -// Copyright (c) Microsoft Corporation. All rights reserved. +// Copyright (c) Microsoft Corporation. // Licensed under the MIT License. // Manages simulators in a vector of pointers to simulators -#include "factory.hpp" -#include "SparseSimulator.h" -#include "types.h" #include #include +#include "factory.hpp" +#include "SparseSimulator.h" +#include "types.h" -namespace Microsoft -{ -namespace Quantum -{ -namespace SPARSESIMULATOR +namespace Microsoft::Quantum::SPARSESIMULATOR { + // Ensures exclusive access to _simulators, the vector of simulators std::shared_mutex _mutex; std::vector> _simulators; @@ -23,7 +20,7 @@ std::vector> _simulators; unsigned createSimulator(logical_qubit_id num_qubits) { if (num_qubits > MAX_QUBITS) - throw std::runtime_error("Max number of qubits is 1024!"); + throw std::runtime_error("Max number of qubits is exceeded!"); std::lock_guard lock(_mutex); size_t emptySlot = -1; for (auto const& s : _simulators) @@ -63,6 +60,4 @@ std::shared_ptr& getSimulator(unsigned id) return _simulators[id]; } -} // namespace Simulator -} // namespace Quantum -} // namespace Microsoft +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp b/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp index 52b95f5c64a..ffb54385f22 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp +++ b/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp @@ -1,22 +1,19 @@ -// Copyright (c) Microsoft Corporation. All rights reserved. +// Copyright (c) Microsoft Corporation. // Licensed under the MIT License. // Manages simulators in a vector of pointers to simulators #pragma once + #include "types.h" #include "SparseSimulator.h" -namespace Microsoft -{ -namespace Quantum -{ -namespace SPARSESIMULATOR +namespace Microsoft::Quantum::SPARSESIMULATOR { + unsigned createSimulator(logical_qubit_id); void destroySimulator(unsigned); std::shared_ptr& getSimulator(unsigned); -} // namespace Simulator -} // namespace Quantum -} // namespace Microsoft + +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/gates.h b/src/Simulation/Simulators/SparseSimulator/Native/gates.h index 5b05870fa21..6871f2a9e2d 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/gates.h +++ b/src/Simulation/Simulators/SparseSimulator/Native/gates.h @@ -1,284 +1,279 @@ -#pragma once +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. +#pragma once -namespace Microsoft +namespace Microsoft::Quantum::SPARSESIMULATOR { - namespace Quantum - { - namespace SPARSESIMULATOR - { - // Used to track differnet kind of gates, - // mainly for timing summary data - enum OP { - MCPhase, - Phase, - T, - AdjT, - S, - AdjS, - H, - MCH, - X, - MCX, - Y, - MCY, - Z, - MCZ, - M, - Measure, - Exp, - MCExp, - R1, - MCR1, - Rx, - Rz, - Ry, - MCRx, - MCRz, - MCRy, - SWAP, - MCSWAP, - Assert, - PermuteSmall, - PermuteLarge, - Proj, - MCProj, - Allocate, - Release, - NUM_OPS // counts gate types; do not add gates after this! - }; - // Strings for the names of the various gate types - const static std::string op_name(OP gate_type) { - switch (gate_type) { - case OP::MCPhase: - return "MCPhase"; - case OP::Phase: - return "Phase"; - case OP::T: - return "T"; - case OP::AdjT: - return "AdjT"; - case OP::S: - return "S"; - case OP::AdjS: - return "AdjS"; - case OP::H: - return "H"; - case OP::MCH: - return "MCH"; - case OP::X: - return "X"; - case OP::MCX: - return "MCX"; - case OP::Y: - return "Y"; - case OP::MCY: - return "MCY"; - case OP::Z: - return "Z"; - case OP::MCZ: - return "MCZ"; - case OP::M: - return "M"; - case OP::Measure: - return "Measure"; - case OP::Exp: - return "Exp"; - case OP::MCExp: - return "MCExp"; - case OP::R1: - return "R1"; - case OP::MCR1: - return "MCR1"; - case OP::Rx: - return "Rx"; - case OP::Rz: - return "Rz"; - case OP::Ry: - return "Ry"; - case OP::MCRx: - return "MCRx"; - case OP::MCRz: - return "MCRz"; - case OP::MCRy: - return "MCRy"; - case OP::SWAP: - return "SWAP"; - case OP::MCSWAP: - return "MCSWAP"; - case OP::Assert: - return "Assert"; - case OP::PermuteSmall: - return "Perm_S"; - case OP::PermuteLarge: - return "Perm_L"; - case OP::Proj: - return "Project"; - case OP::MCProj: - return "MCProj"; - case OP::Allocate: - return "Alloc"; - case OP::Release: - return "Release"; - default: - return "Not a gate"; - } - } +// Used to track differnet kind of gates, +// mainly for timing summary data +enum class OP { + MCPhase, + Phase, + T, + AdjT, + S, + AdjS, + H, + MCH, + X, + MCX, + Y, + MCY, + Z, + MCZ, + M, + Measure, + Exp, + MCExp, + R1, + MCR1, + Rx, + Rz, + Ry, + MCRx, + MCRz, + MCRy, + SWAP, + MCSWAP, + Assert, + PermuteSmall, + PermuteLarge, + Proj, + MCProj, + Allocate, + Release, + NUM_OPS // counts gate types; do not add gates after this! +}; + +// Strings for the names of the various gate types +const static std::string op_name(OP gate_type) { + switch (gate_type) { + case OP::MCPhase: + return "MCPhase"; + case OP::Phase: + return "Phase"; + case OP::T: + return "T"; + case OP::AdjT: + return "AdjT"; + case OP::S: + return "S"; + case OP::AdjS: + return "AdjS"; + case OP::H: + return "H"; + case OP::MCH: + return "MCH"; + case OP::X: + return "X"; + case OP::MCX: + return "MCX"; + case OP::Y: + return "Y"; + case OP::MCY: + return "MCY"; + case OP::Z: + return "Z"; + case OP::MCZ: + return "MCZ"; + case OP::M: + return "M"; + case OP::Measure: + return "Measure"; + case OP::Exp: + return "Exp"; + case OP::MCExp: + return "MCExp"; + case OP::R1: + return "R1"; + case OP::MCR1: + return "MCR1"; + case OP::Rx: + return "Rx"; + case OP::Rz: + return "Rz"; + case OP::Ry: + return "Ry"; + case OP::MCRx: + return "MCRx"; + case OP::MCRz: + return "MCRz"; + case OP::MCRy: + return "MCRy"; + case OP::SWAP: + return "SWAP"; + case OP::MCSWAP: + return "MCSWAP"; + case OP::Assert: + return "Assert"; + case OP::PermuteSmall: + return "Perm_S"; + case OP::PermuteLarge: + return "Perm_L"; + case OP::Proj: + return "Project"; + case OP::MCProj: + return "MCProj"; + case OP::Allocate: + return "Alloc"; + case OP::Release: + return "Release"; + default: + return "Not a gate"; + } +} - // Used in operation queues for phases/permutations - // Different constructors correspond to different - // kinds of gates, so some data is not initialized - // for certain types of gates - struct operation { - OP gate_type; - logical_qubit_id target; - operation(OP gate_type_arg, - logical_qubit_id target_arg) : - gate_type(gate_type_arg), - target(target_arg) {} +// Used in operation queues for phases/permutations +// Different constructors correspond to different +// kinds of gates, so some data is not initialized +// for certain types of gates +struct operation { + OP gate_type; + logical_qubit_id target; + operation(OP gate_type_arg, + logical_qubit_id target_arg) : + gate_type(gate_type_arg), + target(target_arg) {} - std::vector controls; - operation(OP gate_type_arg, - logical_qubit_id target_arg, - std::vector controls_arg - ) : gate_type(gate_type_arg), - target(target_arg), - controls(controls_arg){} + std::vector controls; + operation(OP gate_type_arg, + logical_qubit_id target_arg, + std::vector controls_arg + ) : gate_type(gate_type_arg), + target(target_arg), + controls(controls_arg){} - logical_qubit_id shift; - logical_qubit_id target_2; - //swap - operation (OP gate_type_arg, - logical_qubit_id target1_arg, - logical_qubit_id shift_arg, - logical_qubit_id target_2_arg - ) :gate_type(gate_type_arg), - target(target1_arg), - shift(shift_arg), - target_2(target_2_arg){} - //mcswap - operation(OP gate_type_arg, - logical_qubit_id target1_arg, - logical_qubit_id shift_arg, - std::vector controls_arg, - logical_qubit_id target_2_arg - ) : gate_type(gate_type_arg), - target(target1_arg), - shift(shift_arg), - controls(controls_arg), - target_2(target_2_arg){} - amplitude phase; - // Phase - operation(OP gate_type_arg, - logical_qubit_id target_arg, - amplitude phase_arg - ) :gate_type(gate_type_arg), - target(target_arg), - phase(phase_arg){} - // MCPhase - operation(OP gate_type_arg, - logical_qubit_id target_arg, - std::vector controls_arg, - amplitude phase_arg - ) : gate_type(gate_type_arg), - target(target_arg), - controls(controls_arg), - phase(phase_arg) - {} - bool result; - // Assert - operation(OP gate_type_arg, - std::vector Zs_arg, - bool result_arg - ) : target(0), - gate_type(gate_type_arg), - controls(Zs_arg), - result(result_arg) - {} - }; + logical_qubit_id shift; + logical_qubit_id target_2; + //swap + operation (OP gate_type_arg, + logical_qubit_id target1_arg, + logical_qubit_id shift_arg, + logical_qubit_id target_2_arg + ) :gate_type(gate_type_arg), + target(target1_arg), + shift(shift_arg), + target_2(target_2_arg){} + //mcswap + operation(OP gate_type_arg, + logical_qubit_id target1_arg, + logical_qubit_id shift_arg, + std::vector controls_arg, + logical_qubit_id target_2_arg + ) : gate_type(gate_type_arg), + target(target1_arg), + shift(shift_arg), + controls(controls_arg), + target_2(target_2_arg){} + amplitude phase; + // Phase + operation(OP gate_type_arg, + logical_qubit_id target_arg, + amplitude phase_arg + ) :gate_type(gate_type_arg), + target(target_arg), + phase(phase_arg){} + // MCPhase + operation(OP gate_type_arg, + logical_qubit_id target_arg, + std::vector controls_arg, + amplitude phase_arg + ) : gate_type(gate_type_arg), + target(target_arg), + controls(controls_arg), + phase(phase_arg) + {} + bool result; + // Assert + operation(OP gate_type_arg, + std::vector Zs_arg, + bool result_arg + ) : target(0), + gate_type(gate_type_arg), + controls(Zs_arg), + result(result_arg) + {} +}; - // Also represents operations, but uses - // bitsets instead of vectors of qubit ids - // to save time/space - template - struct condensed_operation { - OP gate_type; - logical_qubit_id target; - condensed_operation(OP gate_type_arg, - logical_qubit_id target_arg) : - gate_type(gate_type_arg), - target(target_arg) - {} +// Also represents operations, but uses +// bitsets instead of vectors of qubit ids +// to save time/space +template +struct condensed_operation { + OP gate_type; + logical_qubit_id target; + condensed_operation(OP gate_type_arg, + logical_qubit_id target_arg) : + gate_type(gate_type_arg), + target(target_arg) + {} - std::bitset controls; - condensed_operation(OP gate_type_arg, - logical_qubit_id target_arg, - std::bitset controls_arg - ) : gate_type(gate_type_arg), - target(target_arg), - controls(controls_arg){} + std::bitset controls; + condensed_operation(OP gate_type_arg, + logical_qubit_id target_arg, + std::bitset controls_arg + ) : gate_type(gate_type_arg), + target(target_arg), + controls(controls_arg){} - logical_qubit_id target_2; - //swap - condensed_operation (OP gate_type_arg, - logical_qubit_id target1_arg, - logical_qubit_id target_2_arg - ) :gate_type(gate_type_arg), - target(target1_arg), - target_2(target_2_arg){} - //mcswap - condensed_operation(OP gate_type_arg, - logical_qubit_id target1_arg, - std::bitset controls_arg, - logical_qubit_id target_2_arg - ) : gate_type(gate_type_arg), - target(target1_arg), - controls(controls_arg), - target_2(target_2_arg){} - amplitude phase; - // Phase - condensed_operation(OP gate_type_arg, - logical_qubit_id target_arg, - amplitude phase_arg - ) :gate_type(gate_type_arg), - target(target_arg), - phase(phase_arg){} - // MCPhase - condensed_operation(OP gate_type_arg, - logical_qubit_id target_arg, - std::bitset controls_arg, - amplitude phase_arg - ) : gate_type(gate_type_arg), - target(target_arg), - controls(controls_arg), - phase(phase_arg) - {} - bool result; - // Assert - condensed_operation(OP gate_type_arg, - std::bitset Zs_arg, - bool result_arg - ) : target(0), - gate_type(gate_type_arg), - controls(Zs_arg), - result(result_arg){} - }; - - namespace Gates - { - /// a type for runtime basis specification - enum class Basis - { - PauliI = 0, - PauliX = 1, - PauliY = 3, - PauliZ = 2 - }; - + logical_qubit_id target_2; + //swap + condensed_operation (OP gate_type_arg, + logical_qubit_id target1_arg, + logical_qubit_id target_2_arg + ) :gate_type(gate_type_arg), + target(target1_arg), + target_2(target_2_arg){} + //mcswap + condensed_operation(OP gate_type_arg, + logical_qubit_id target1_arg, + std::bitset controls_arg, + logical_qubit_id target_2_arg + ) : gate_type(gate_type_arg), + target(target1_arg), + controls(controls_arg), + target_2(target_2_arg){} + amplitude phase; + // Phase + condensed_operation(OP gate_type_arg, + logical_qubit_id target_arg, + amplitude phase_arg + ) :gate_type(gate_type_arg), + target(target_arg), + phase(phase_arg){} + // MCPhase + condensed_operation(OP gate_type_arg, + logical_qubit_id target_arg, + std::bitset controls_arg, + amplitude phase_arg + ) : gate_type(gate_type_arg), + target(target_arg), + controls(controls_arg), + phase(phase_arg) + {} + bool result; + // Assert + condensed_operation(OP gate_type_arg, + std::bitset Zs_arg, + bool result_arg + ) : target(0), + gate_type(gate_type_arg), + controls(Zs_arg), + result(result_arg){} +}; +namespace Gates +{ + /// a type for runtime basis specification + enum class Basis + { + PauliI = 0, + PauliX = 1, + PauliY = 3, + PauliZ = 2 + }; +} // namespace Gates - } - } - } -} \ No newline at end of file +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp b/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp index 77f9692b40c..01a1b7d10f7 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp +++ b/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp @@ -1,1511 +1,1506 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + #pragma once + #include #include -#include "types.h" -#include "gates.h" #include #include #include -#include "basic_quantum_state.hpp" -#include "quantum_hash_map.hpp" #include #include #include - #include #include - #ifdef _OPENMP - #include + #include #endif +#include "basic_quantum_state.hpp" +#include "quantum_hash_map.hpp" +#include "types.h" +#include "gates.h" using namespace std::literals::complex_literals; -namespace Microsoft -{ -namespace Quantum +namespace Microsoft::Quantum::SPARSESIMULATOR { -namespace SPARSESIMULATOR -{ - // power of square root of -1 inline amplitude iExp(int power) { - int p = ((power % 4) + 8) % 4; - switch (p) - { - case 0: - return 1; - case 1: - return 1i; - case 2: - return -1; - case 3: - return -1i; - default: - return 0; - } - return 0; + int p = ((power % 4) + 8) % 4; + switch (p) + { + case 0: + return 1; + case 1: + return 1i; + case 2: + return -1; + case 3: + return -1i; + default: + return 0; + } + return 0; } template bool poppar(std::bitset bitstring){ - return bitstring.count() % 2; + return bitstring.count() % 2; } // Compares two bitsets as through they were bitstrings // Used to enforce an ordering on bitsets, though currently not referenced template inline bool operator<(const std::bitset& lhs, const std::bitset& rhs) { - std::bitset mask = lhs ^ rhs; - std::bitset const ull_mask = std::bitset((unsigned long long) -1); - for (int i = static_cast(N - 8*sizeof(unsigned long long)); i > 0; i-= static_cast(8*sizeof(unsigned long long))){ - if (((mask >> i) & ull_mask).to_ullong() > 0){ - return ((lhs >> i) & ull_mask).to_ullong() < ((rhs >> i) & ull_mask).to_ullong(); - } - } - return ((lhs) & ull_mask).to_ullong() < ((rhs) & ull_mask).to_ullong(); + std::bitset mask = lhs ^ rhs; + std::bitset const ull_mask = std::bitset((unsigned long long) -1); + for (int i = static_cast(N - 8*sizeof(unsigned long long)); i > 0; i-= static_cast(8*sizeof(unsigned long long))){ + if (((mask >> i) & ull_mask).to_ullong() > 0){ + return ((lhs >> i) & ull_mask).to_ullong() < ((rhs >> i) & ull_mask).to_ullong(); + } + } + return ((lhs) & ull_mask).to_ullong() < ((rhs) & ull_mask).to_ullong(); } - + // Transforms a vector of indices into a bitset where the indices indicate precisely // which bits are non-zero template std::bitset get_mask(std::vector const& indices){ - std::bitset mask = std::bitset(); - for (logical_qubit_id index : indices) { - mask.set(index); - } - return mask; + std::bitset mask = std::bitset(); + for (logical_qubit_id index : indices) { + mask.set(index); + } + return mask; } template class QuantumState : public BasicQuantumState { public: - // Type for qubit labels, with a specific size built-in - using qubit_label = qubit_label_type; - - // Type of hash maps with the required labels - using wavefunction = abstract_wavefunction; - - QuantumState() { - _qubit_data = wavefunction(); - _qubit_data.max_load_factor(_load_factor); - // Create an initial all-zeros state - _qubit_data.emplace((logical_qubit_id)0, 1); - // Initialize randomness - std::random_device rd; - std::mt19937 gen(rd()); - std::uniform_real_distribution dist(0, 1); - _rng = [gen, dist]() mutable { return dist(gen); }; - - _initialize_threads(); - } - - ~QuantumState() { - complete_threads(); - } - - // Copy data from an existing simulator - // This is used to move between different qubit sizes - // without needing a lot of templated functions - QuantumState(std::shared_ptr old_state) { - // Copy any needed data - _rng = old_state->get_rng(); - // Outputs the previous data with labels as strings - universal_wavefunction old_qubit_data = old_state->get_universal_wavefunction(); - _qubit_data = wavefunction(old_qubit_data.size()); - _load_factor = old_state->get_load_factor(); - _qubit_data.max_load_factor(_load_factor); - // Writes this into the current wavefunction as qubit_label types - for (auto current_state = old_qubit_data.begin(); current_state != old_qubit_data.end(); ++current_state) { - _qubit_data.emplace(qubit_label(current_state->first), current_state->second); - } - // Create local threads for this state - _initialize_threads(); - } - - logical_qubit_id get_num_qubits() { - return (logical_qubit_id)num_qubits; - } - - // Outputs all states and amplitudes to the console - void DumpWavefunction(size_t indent = 0){ - DumpWavefunction(_qubit_data, indent); - } - - // Outputs all states and amplitudes from an input wavefunction to the console - void DumpWavefunction(wavefunction &wfn, size_t indent = 0){ - std::string spacing(indent, ' '); - std::cout << spacing << "Wavefunction:\n"; - auto line_dump = [spacing](qubit_label label, amplitude val){ - std::cout << spacing << " " << label.to_string() << ": "; - std::cout << val.real(); - std::cout << (val.imag() < 0 ? " - " : " + ") << std::abs(val.imag()) << "i\n"; - }; - _DumpWavefunction_base(wfn, line_dump); - std::cout << spacing << "--end wavefunction\n"; - } - - - void set_random_seed(unsigned seed) { - std::mt19937 gen(seed); - std::uniform_real_distribution dist(0, 1); - _rng = [gen, dist]() mutable { return dist(gen); }; - } - - // Used to decide when an amplitude is close enough to 0 to discard - void set_precision(double new_precision) { - _precision = new_precision; - _precision_squared = _precision *_precision; - } - - // Load factor of the underlying hash map - float get_load_factor() { - return _load_factor; - } - - void set_load_factor(float new_load_factor) { - _load_factor = new_load_factor; - } - - // Returns the number of states in superposition - size_t get_wavefunction_size() { - return _qubit_data.size(); - } - - - - // Applies the operator id_coeff*I + pauli_coeff * P - // where P is the Pauli operators defined by axes applied to the qubits in qubits. - void PauliCombination(std::vector const& axes, std::vector const& qubits, amplitude id_coeff, amplitude pauli_coeff) { - // Bit-vectors indexing where gates of each type are applied - qubit_label XYs = 0; - qubit_label YZs = 0; - logical_qubit_id ycount = 0; - for (int i=0; i < axes.size(); i++){ - switch (axes[i]){ - case Gates::Basis::PauliY: - YZs.set(qubits[i]); - XYs.set(qubits[i]); - ycount++; - break; - case Gates::Basis::PauliX: - XYs.set(qubits[i]); - break; - case Gates::Basis::PauliZ: - YZs.set(qubits[i]); - break; - case Gates::Basis::PauliI: - break; - default: - throw std::runtime_error("Bad Pauli basis"); - } - } - - // All identity - if (XYs.none() && YZs.none()) { - return; - } - - // This branch handles purely Z Pauli vectors - // Purely Z has no addition, which would cause - // problems in the comparison in the next section - if (XYs.none()) { - // 0 terms get the sum of the coefficients - // 1 terms get the difference - pauli_coeff += id_coeff; // id_coeff + pauli_coeff - id_coeff *= 2; - id_coeff -= pauli_coeff; // id_coeff - pauli_coeff - - // To avoid saving states of zero amplitude, these if/else - // check for when one of the coefficients is - // close enough to zero to regard as zero - if (std::norm(pauli_coeff) > _rotation_precision_squared ){ - if (std::norm(id_coeff) > _rotation_precision_squared){ - // If both coefficients are non-zero, we can just modify the state in-place - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - current_state->second *= (poppar(current_state->first & YZs) ? id_coeff : pauli_coeff); - } - } else { - // If id_coeff = 0, then we make a new wavefunction and only add in those that will be multiplied - // by the pauli_coeff - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if (!poppar(current_state->first & YZs)){ - new_qubit_data.emplace(current_state->first, current_state->second * pauli_coeff); - } - } - _qubit_data = std::move(new_qubit_data); - } - } else { - // If pauli_coeff=0, don't add states multiplied by the pauli_coeff - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if (poppar(current_state->first & YZs)){ - new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); - } - } - _qubit_data = std::move(new_qubit_data); - } - } else { // There are some X or Y gates - - // Each Y Pauli adds a global phase of i - switch (ycount % 4) { - case 1: - pauli_coeff *= 1i; - break; - case 2: - pauli_coeff *= -1; - break; - case 3: - pauli_coeff *= -1i; - break; - default: - break; - } - // When both the state and flipped state are in superposition, when adding the contribution of - // the flipped state, we add phase depending on the 1s in the flipped state - // This phase would be the parity of (flipped_state->first ^ YZs) - // However, we know that flipped_state->first = current_state->first ^ YXs - // So the parity of the flipped state will be the parity of the current state, plus - // the parity of YZs & YXs, i.e., the parity of the number of Ys - // Since this is constant for all states, we compute it once here and save it - // Then we only compute the parity of the current state - amplitude pauli_coeff_alt = ycount % 2 ? -pauli_coeff : pauli_coeff; - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); - amplitude new_state; - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - auto alt_state = _qubit_data.find(current_state->first ^ XYs); - if (alt_state == _qubit_data.end()) { // no matching value - new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); - new_qubit_data.emplace(current_state->first ^ XYs, current_state->second * (poppar(current_state->first & YZs) ? -pauli_coeff : pauli_coeff)); - } - else if (current_state->first < alt_state->first) { - // Each Y and Z gate adds a phase (since Y=iXZ) - bool parity = poppar(current_state->first & YZs); - new_state = current_state->second * id_coeff + alt_state->second * (parity ? -pauli_coeff_alt : pauli_coeff_alt); - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first, new_state); - } - - new_state = alt_state->second * id_coeff + current_state->second * (parity ? -pauli_coeff : pauli_coeff); - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(alt_state->first, new_state); - } - } - } - _qubit_data = std::move(new_qubit_data); - } - } - - // Applies the operator id_coeff*I + pauli_coeff * P - // where P is the Pauli operators defined by axes applied to the qubits in qubits. - // Controlled version - void MCPauliCombination(std::vector const& controls, std::vector const& axes, std::vector const& qubits, amplitude id_coeff, amplitude pauli_coeff) { - // Bit-vectors indexing where gates of each type are applied - qubit_label cmask = _get_mask(controls); - qubit_label XYs = 0; - qubit_label YZs = 0; - logical_qubit_id ycount = 0; - // Used for comparing pairs - logical_qubit_id any_xy = -1; - for (int i=0; i < axes.size(); i++){ - switch (axes[i]){ - case Gates::Basis::PauliY: - YZs.set(qubits[i]); - XYs.set(qubits[i]); - ycount++; - any_xy = qubits[i]; - break; - case Gates::Basis::PauliX: - XYs.set(qubits[i]); - any_xy = qubits[i]; - break; - case Gates::Basis::PauliZ: - YZs.set(qubits[i]); - break; - case Gates::Basis::PauliI: - break; - default: - throw std::runtime_error("Bad Pauli basis"); - } - } - - // This branch handles purely Z Pauli vectors - // Purely Z has no addition, which would cause - // problems in the comparison in the next section - if (XYs.none()) { - // 0 terms get the sum of the coefficients - // 1 terms get the difference - pauli_coeff += id_coeff; // <- id_coeff + pauli_coeff - id_coeff *= 2; - id_coeff -= pauli_coeff; // <- id_coeff - pauli_coeff - - // To avoid saving states of zero amplitude, these if/else - // check for when one of the coefficients is - // close enough to zero to regard as zero - if (std::norm(pauli_coeff) > _rotation_precision_squared ){ - if (std::norm(id_coeff) > _rotation_precision_squared){ - // If both coefficients are non-zero, we can just modify the state in-place - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if ((current_state->first & cmask)==cmask) { - current_state->second *= (poppar(current_state->first & YZs) ? id_coeff : pauli_coeff); - } - } - } else { - // If id_coeff = 0, then we make a new wavefunction and only add in those that will be multiplied - // by the pauli_coeff - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if (!poppar(current_state->first & YZs) && (current_state->first & cmask)==cmask){ - new_qubit_data.emplace(current_state->first, current_state->second * pauli_coeff); - } - } - _qubit_data = std::move(new_qubit_data); - } - } else { - // If pauli_coeff=0, don't add states multiplied by the pauli_coeff - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if (poppar(current_state->first & YZs) && (current_state->first & cmask)==cmask){ - new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); - } - } - _qubit_data = std::move(new_qubit_data); - } - } else { // There are some X or Y gates - // Each Y Pauli adds a global phase of i - switch (ycount % 4) { - case 1: - pauli_coeff *= 1i; - break; - case 2: - pauli_coeff *= -1; - break; - case 3: - pauli_coeff *= -1i; - break; - default: - break; - } - // When both the state and flipped state are in superposition, when adding the contribution of - // the flipped state, we add phase depending on the 1s in the flipped state - // This phase would be the parity of (flipped_state->first ^ YZs) - // However, we know that flipped_state->first = current_state->first ^ YXs - // So the parity of the flipped state will be the parity of the current state, plus - // the parity of YZs & YXs, i.e., the parity of the number of Ys - // Since this is constant for all states, we compute it once here and save it - // Then we only compute the parity of the current state - amplitude pauli_coeff_alt = ycount % 2 ? -pauli_coeff : pauli_coeff; - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); - amplitude new_state; - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if ((current_state->first & cmask)==cmask) { - auto alt_state = _qubit_data.find(current_state->first ^ XYs); - if (alt_state == _qubit_data.end()) { // no matching value - new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); - new_qubit_data.emplace(current_state->first ^ XYs, current_state->second * (poppar(current_state->first & YZs) ? -pauli_coeff : pauli_coeff)); - } - else if (current_state->first < alt_state->first) { //current_state->first[any_xy]){// - // Each Y and Z gate adds a phase (since Y=iXZ) - bool parity = poppar(current_state->first & YZs); - new_state = current_state->second * id_coeff + alt_state->second * (parity ? -pauli_coeff_alt : pauli_coeff_alt); - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first, new_state); - } - - new_state = alt_state->second * id_coeff + current_state->second * (parity ? -pauli_coeff : pauli_coeff); - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(alt_state->first, new_state); - } - } - } else { - new_qubit_data.emplace(current_state->first, current_state->second); - } - } - _qubit_data = std::move(new_qubit_data); - } - } - - - bool M(logical_qubit_id target) { - qubit_label flip = qubit_label(); - flip.set(target); - - bool result = _qubit_data.begin()->first[target]; - - double zero_probability = 0.0; - double one_probability = 0.0; - - // Writes data into a ones or zeros wavefunction - // as it adds up probability - // Once it's finished, it picks one randomly, normalizes - // then keeps that one as the new wavefunction - wavefunction ones = make_wavefunction(_qubit_data.size()/2); - wavefunction zeros = make_wavefunction(_qubit_data.size()/2); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - double square_amplitude = std::norm(current_state->second); - if (current_state->first[target]) { - one_probability += square_amplitude; - ones.emplace(current_state->first, current_state->second); - } - else { - zero_probability += square_amplitude; - zeros.emplace(current_state->first, current_state->second); - } - } - // Randomly select - result = (_rng() <= one_probability); - - wavefunction &new_qubit_data = result ? ones : zeros; - // Create a new, normalized state - double normalizer = 1.0/std::sqrt((result) ? one_probability : zero_probability); - for (auto current_state = (new_qubit_data).begin(); current_state != (new_qubit_data).end(); ++current_state) { - current_state->second *= normalizer; - } - _qubit_data = std::move(new_qubit_data); - - return result; - } - - void Reset(logical_qubit_id target) { - qubit_label flip = qubit_label(0); - flip.set(target); - - double zero_probability = 0.0; - double one_probability = 0.0; - - // Writes data into a ones or zeros wavefunction - // as it adds up probability - // Once it's finished, it picks one randomly, normalizes - // then keeps that one as the new wavefunction - - // Used to set the qubit to 0 in the measured result - qubit_label new_mask = qubit_label(); - new_mask.set(); // sets all bits to 1 - new_mask.set(target, 0); - wavefunction ones = make_wavefunction(_qubit_data.size()/2); - wavefunction zeros = make_wavefunction(_qubit_data.size()/2); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - double square_amplitude = std::norm(current_state->second); - if (current_state->first[target]) { - one_probability += square_amplitude; - ones.emplace(current_state->first & new_mask, current_state->second); - } - else { - zero_probability += square_amplitude; - zeros.emplace(current_state->first & new_mask, current_state->second); - } - } - // Randomly select - bool result = (_rng() <= one_probability); - - wavefunction &new_qubit_data = result ? ones : zeros; - // Create a new, normalized state - double normalizer = 1.0/std::sqrt((result) ? one_probability : zero_probability); - for (auto current_state = (new_qubit_data).begin(); current_state != (new_qubit_data).end(); ++current_state) { - current_state->second *= normalizer; - } - _qubit_data = std::move(new_qubit_data); - } - - - // Samples a state from the superposition with probably proportion to - // the amplitude, returning a string of the bits of that state. - // Unlike measurement, this does not modify the state - std::string Sample() { - double probability = _rng(); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - double square_amplitude = std::norm(current_state->second); - probability -= square_amplitude; - if (probability <= 0){ - return current_state->first.to_string(); - } - } - return _qubit_data.begin()->first.to_string(); - } - - void Assert(std::vector const& axes, std::vector const& qubits, bool result) { - // Bit-vectors indexing where gates of each type are applied - qubit_label XYs = 0; - qubit_label YZs = 0; - logical_qubit_id ycount = 0; - for (int i=0; i < axes.size(); i++){ - switch (axes[i]){ - case Gates::Basis::PauliY: - YZs.set(qubits[i]); - XYs.set(qubits[i]); - ycount++; - break; - case Gates::Basis::PauliX: - XYs.set(qubits[i]); - break; - case Gates::Basis::PauliZ: - YZs.set(qubits[i]); - break; - case Gates::Basis::PauliI: - break; - default: - throw std::runtime_error("Bad Pauli basis"); - } - } - - amplitude phaseShift = result ? -1 : 1; - // Each Y Pauli adds a global phase of i - switch (ycount % 4) { - case 1: - phaseShift *= 1i; - break; - case 2: - phaseShift *= -1; - break; - case 3: - phaseShift *= -1i; - break; - default: - break; - } - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if (std::norm(_qubit_data.find(current_state->first ^ XYs)->second - current_state->second * (poppar(current_state->first & YZs) ? -phaseShift : phaseShift)) > _precision_squared) { - qubit_label label = current_state->first; - amplitude val = current_state->second; - std::cout << "Problematic state: " << label << "\n"; - std::cout << "Expected " << val * (poppar(label & YZs) ? -phaseShift : phaseShift); - std::cout << ", got " << _qubit_data.find(label ^ XYs)->second << "\n"; - std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; - throw std::runtime_error("Not an eigenstate"); - } - } - } - - // Returns the probability of a given measurement in a Pauli basis - // by decomposing each pair of computational basis states into eigenvectors - // and adding the coefficients of the respective components - double MeasurementProbability(std::vector const& axes, std::vector const& qubits) { - // Bit-vectors indexing where gates of each type are applied - qubit_label XYs = 0; - qubit_label YZs = 0; - logical_qubit_id ycount = 0; - for (int i=0; i < axes.size(); i++){ - switch (axes[i]){ - case Gates::Basis::PauliY: - YZs.set(qubits[i]); - XYs.set(qubits[i]); - ycount++; - break; - case Gates::Basis::PauliX: - XYs.set(qubits[i]); - break; - case Gates::Basis::PauliZ: - YZs.set(qubits[i]); - break; - case Gates::Basis::PauliI: - break; - default: - throw std::runtime_error("Bad Pauli basis"); - } - } - amplitude phaseShift = 1; - - // Each Y Pauli adds a global phase of i - switch (ycount % 4) { - case 1: - phaseShift *= amplitude(0, 1); - break; - case 2: - phaseShift *= -1; - break; - case 3: - phaseShift *= amplitude(0, -1); - break; - default: - break; - } - // Let P be the pauli operation, |psi> the state - // projection = - - // _qubit_data represents |psi> as sum_x a_x |x>, - // where all |x> are orthonormal. Thus, the projection - // will be the product of a_x and a_P(x), where P|x>=|P(x)> - // Thus, for each |x>, we compute P(x) and look for that state - // If there is a match, we add the product of their coefficients - // to the projection, times a phase dependent on how many Ys and Zs match - // the 1 bits of x - amplitude projection = 0.0; - auto flipped_state = _qubit_data.end(); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - flipped_state = _qubit_data.find(current_state->first ^ XYs); // no match returns _qubit_data.end() - projection += current_state->second * (flipped_state == _qubit_data.end() ? 0 : std::conj(flipped_state->second)) * (poppar(current_state->first & YZs) ? -phaseShift : phaseShift); - } - // The projector onto the -1 eigenspace (a result of "One") is 0.5 * (I - P) - // So = 0.5 - 0.5* - // should always be real so this only takes the real part - return 0.5 - 0.5 * projection.real(); - } - - bool Measure(std::vector const& axes, std::vector const& qubits){ - // Find a probability to get a specific result - double probability = MeasurementProbability(axes, qubits); - bool result = _rng() <= probability; - probability = std::sqrt(probability); - // This step executes immediately so that we reduce the number of states in superposition - PauliCombination(axes, qubits, 0.5/probability, (result ? -0.5 : 0.5)/probability); - return result; - } - - - // Probe the amplitude of a single basis state - amplitude probe(qubit_label label) { - auto qubit = _qubit_data.find(label); - // States not in the hash map are assumed to be 0 - if (qubit == _qubit_data.end()) { - return amplitude(0.0, 0.0); - } - else { - return qubit->second; - } - } - - amplitude probe(std::string label) { - qubit_label bit_label = qubit_label(label); - return probe(bit_label); - } - - // Dumps the state of a subspace of particular qubits, if they are not entangled - // This requires it to detect if the subspace is entangled, construct a new - // projected wavefunction, then call the `callback` function on each state. - bool dump_qubits(std::vector qubits, void (*callback)(char*, double, double)) { - // Create two wavefunctions - // check if they are tensor products - wavefunction dump_wfn; - wavefunction leftover_wfn; - if (!_split_wavefunction(_get_mask(qubits), dump_wfn, leftover_wfn)){ - return false; - } else { - _DumpWavefunction_base(dump_wfn, [qubits, callback](qubit_label label, amplitude val){ - std::string label_string(qubits.size(), '0'); - for (size_t i=0; i < qubits.size(); i++){ - label_string[i] = label[qubits[i]] ? '1' : '0'; - } - callback(const_cast(label_string.c_str()), val.real(), val.imag()); - }); - return true; - } - } - - // Dumps all the states in superposition via a callback function - void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) { - _DumpWavefunction_base(_qubit_data, [max_qubit_id, callback](qubit_label label, amplitude val){ - callback(const_cast(label.to_string().substr(num_qubits - 1 - max_qubit_id, max_qubit_id + 1).c_str()), val.real(), val.imag()); - }); - } - - // Execute a queue of phase/permutation gates - void phase_and_permute(std::list const &operation_list){ - if (operation_list.size()==0){return;} - - // Condense the list into a memory-efficient vector with qubit labels - _operation_vector.reserve(operation_list.size()); - for (auto op : operation_list){ - switch (op.gate_type) { - case OP::X: - case OP::Y: - case OP::Z: - _operation_vector.push_back(internal_operation(op.gate_type, op.target)); - break; - case OP::MCX: - case OP::MCY: - _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls))); - break; - case OP::MCZ: - _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls).set(op.target))); - break; - case OP::Phase: - _operation_vector.push_back(internal_operation(op.gate_type, op.target, op.phase)); - break; - case OP::MCPhase: - _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls).set(op.target), op.phase)); - break; - case OP::SWAP: - _operation_vector.push_back(internal_operation(op.gate_type, op.target, op.target_2)); - break; - case OP::MCSWAP: - _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls), op.target_2)); - break; - case OP::Assert: - _operation_vector.push_back(internal_operation(op.gate_type, _get_mask(op.controls), op.result)); - break; - default: - throw std::runtime_error("Unsupported operation"); - break; - } - } - - _new_qubit_data = make_wavefunction(); - - // Threading introduces a lot of overhead so it - // is only worthwhile for large instances - // 64 and 4096 are empirical values - if (_operation_vector.size() < 64 || _qubit_data.size() < 4096){ // small; do not thread - // Iterates through and applies all operations - for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end(); ++current_state){ - qubit_label label = current_state->first; - amplitude val = current_state->second; - // Iterate through vector of operations and apply each gate - for (int i=0; i < _operation_vector.size(); i++) { - auto &op = _operation_vector[i]; - switch (op.gate_type) { - case OP::X: - label.flip(op.target); - break; - case OP::MCX: - if ((op.controls & label) == op.controls){ - label.flip(op.target); - } - break; - case OP::Y: - label.flip(op.target); - val *= (label[op.target]) ? 1i : -1i; - break; - case OP::MCY: - if ((op.controls & label) == op.controls){ - label.flip(op.target); - val *= (label[op.target]) ? 1i : -1i; - } - break; - case OP::Z: - val *= (label[op.target] ? -1 : 1); - break; - case OP::MCZ: - val *= ((op.controls & label) == op.controls) ? -1 : 1; - break; - case OP::Phase: - val *= label[op.target] ? op.phase : 1; - break; - case OP::MCPhase: - val *= ((op.controls & label) == op.controls) ? op.phase : 1; - break; - case OP::SWAP: - if (label[op.target] != label[op.target_2]){ - label.flip(op.target); - label.flip(op.target_2); - } - break; - case OP::MCSWAP: - if (((label & op.controls) == op.controls) && (label[op.target] != label[op.target_2])){ - label.flip(op.target); - label.flip(op.target_2); - } - break; - case OP::Assert: - if (poppar(label & op.controls) != op.result && std::norm(val) > _precision_squared){ - std::cout << "Problematic state: " << label << "\n"; - std::cout << "Amplitude: " << val << "\n"; - std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; - throw std::runtime_error("Assert failed"); - } - break; - default: - throw std::runtime_error("Unsupported operation"); - break; - } - } - // Insert the new state into the new wavefunction - _new_qubit_data.emplace(label, val); - } - _qubit_data = std::move(_new_qubit_data); - _new_qubit_data.clear(); - _operation_vector.clear(); - } else { // Large enough to multi-thread - #ifndef DOMP_GE_V3 // OMP version is too low; uses condition variables - // Lock the mutex so the threads stay asleep during prep - std::unique_lock state_lock(_state_mtx); - _new_qubit_data = wavefunction(_qubit_data.size()); - // jump_size gives a rough guess for how many states each thread should process - // to minimize contention for current_state - _jump_size = std::max((size_t)1 , (size_t)(_qubit_data.size() / (8*_thread_pool.size()))); - // Set the _current state (this allows the threads to wake up) - _current_state = _qubit_data.begin(); - - // Wake up all threads - cv.notify_all(); - // Wait for the number of running threads to be 0 - // This will spuriously wake many times - cv.wait(state_lock, [&](){ - return _current_state == _qubit_data.end() && _running_threads == 0; - }); - // Here all threads are finished - #else - #pragma omp parallel - { - _jump_size = std::max((size_t)1 , (size_t)(_qubit_data.size() / (8*omp_get_num_threads()))); - #pragma omp single - { - auto internal_state = _qubit_data.begin(); - auto internal_end = _qubit_data.end(); - for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end();){ - internal_state = current_state; - current_state.jump_forward(_jump_size); - internal_end = current_state; - #pragma omp task firstprivate(internal_state) firstprivate(internal_end) - { - qubit_label label; - amplitude val; - for (; internal_state != internal_end; ++internal_state){ - label = internal_state->first; - val = internal_state->second; - for (int i=0; i < _operation_vector.size(); i++) { - auto &op = _operation_vector[i]; - switch (op.gate_type) { - case OP::X: - label.flip(op.target); - break; - case OP::MCX: - if ((op.controls & label) == op.controls){ - label.flip(op.target); - } - break; - case OP::Y: - label.flip(op.target); - val *= (label[op.target]) ? 1i : -1i; - break; - case OP::MCY: - if ((op.controls & label) == op.controls){ - label.flip(op.target); - val *= (label[op.target]) ? 1i : -1i; - } - break; - case OP::Z: - val *= (label[op.target] ? -1 : 1); - break; - case OP::MCZ: - val *= ((op.controls & label) == op.controls) ? -1 : 1; - break; - case OP::Phase: - val *= label[op.target] ? op.phase : 1; - break; - case OP::MCPhase: - val *= ((op.controls & label) == op.controls) ? op.phase : 1; - break; - case OP::SWAP: - if (label[op.target] != label[op.target_2]){ - label.flip(op.target); - label.flip(op.target_2); - } - break; - case OP::MCSWAP: - if (((label & op.controls) == op.controls) && (label[op.target] != label[op.target_2])){ - label.flip(op.target); - label.flip(op.target_2); - } - break; - case OP::Assert: - if (poppar(label & op.controls) != op.result && std::norm(val) > _precision_squared){ - std::cout << "Problematic state: " << label << "\n"; - std::cout << "Amplitude: " << val << "\n"; - std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; - throw std::runtime_error("Assert failed"); - } - break; - default: - throw std::runtime_error("Unsupported operation"); - break; - } - } - #pragma omp critical - { - _new_qubit_data.emplace(label, val); - } - } - } - - } - } - #pragma omp barrier - } - #endif - _qubit_data = std::move(_new_qubit_data); - _operation_vector.clear(); - _new_qubit_data.clear(); - } - } - - void R(Gates::Basis b, double phi, logical_qubit_id index){ - // Z rotation can be done in-place - if (b == Gates::Basis::PauliZ) { - amplitude exp_0 = amplitude(std::cos(phi / 2.0), -std::sin(phi / 2.0)); - amplitude exp_1 = amplitude(std::cos(phi / 2.0), std::sin(phi / 2.0)); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - current_state->second *= current_state->first[index] ? exp_1 : exp_0; - } - } - else if (b == Gates::Basis::PauliX || b == Gates::Basis::PauliY) { - amplitude M00 = std::cos(phi / 2.0); - amplitude M01 = -std::sin(phi / 2.0) * (b == Gates::Basis::PauliY ? 1 : 1i); - if (std::norm(M00) < _rotation_precision_squared){ - // This is just a Y or X gate - phase_and_permute(std::list{operation(b==Gates::Basis::PauliY ? OP::Y : OP::X, index)}); - return; - } else if (std::norm(M01) < _rotation_precision_squared){ - // just an identity - return; - } - - amplitude M10 = M01 * amplitude(b == Gates::Basis::PauliY ? -1 : 1); - // Holds the amplitude of the new state to make it easier to check if it's non-zero - amplitude new_state; - qubit_label flip(0); - flip.set(index); - wavefunction new_qubit_data = make_wavefunction(); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - auto flipped_state = _qubit_data.find(current_state->first ^ flip); - if (flipped_state == _qubit_data.end()) { // no matching value - if (current_state->first[index]) {// 1 on that qubit - new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M01); - new_qubit_data.emplace(current_state->first, current_state->second * M00); - } - else { - new_qubit_data.emplace(current_state->first, current_state->second * M00); - new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M10); - } - } - // Add up the two values, only when reaching the zero value - else if (!(current_state->first[index])) { - new_state = current_state->second * M00 + flipped_state->second * M01; // zero state - if (std::norm(new_state) > _rotation_precision_squared ) { - new_qubit_data.emplace(current_state->first, new_state); - } - new_state = current_state->second * M10 + flipped_state->second * M00; // one state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(flipped_state->first, new_state); - } - } - } - _qubit_data = std::move(new_qubit_data); - } - } - - // Multi-controlled rotation - void MCR (std::vector const& controls, Gates::Basis b, double phi, logical_qubit_id target) { - qubit_label checks = _get_mask(controls); - // A Z-rotation can be done without recreating the wavefunction - if (b == Gates::Basis::PauliZ) { - amplitude exp_0 = amplitude(std::cos(phi / 2.0), -std::sin(phi / 2.0)); - amplitude exp_1 = amplitude(std::cos(phi / 2.0), std::sin(phi / 2.0)); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if ((current_state->first & checks)==checks){ - current_state->second *= current_state->first[target] ? exp_1 : exp_0; - } - } - } - // X or Y requires a new wavefunction - else if (b == Gates::Basis::PauliX || b == Gates::Basis::PauliY) { - amplitude M00 = std::cos(phi / 2.0); - amplitude M01 = -std::sin(phi / 2.0) * (b == Gates::Basis::PauliY ? 1 : 1i); - amplitude M10 = (b == Gates::Basis::PauliY ? -1.0 : 1.0) * M01; - - if (std::norm(M00) < _rotation_precision_squared){ - // This is just an MCY or MCX gate, but with a phase - // So we need to preprocess with a multi-controlled phase - if (b==Gates::Basis::PauliY){ - amplitude phase = -1i/M01; - phase_and_permute(std::list{ - operation(OP::MCPhase, controls[0], controls, phase), - operation(OP::MCY, target, controls) - }); - } else { - amplitude phase = 1.0/M01; - phase_and_permute(std::list{ - operation(OP::MCPhase, controls[0], controls, phase), - operation(OP::MCY, target, controls) - }); - } - return; - } else if (std::norm(M01) < _rotation_precision_squared){ - // This is equivalent to a multi-controlled Z if the rotation is -1 - if (std::norm(M01 + 1.0) < _rotation_precision_squared){ - phase_and_permute(std::list{operation(OP::MCZ, controls[0], controls)}); - } - return; - } - - amplitude new_state; - qubit_label flip(0); - flip.set(target); - wavefunction new_qubit_data = make_wavefunction(); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if ((current_state->first & checks)==checks){ - auto flipped_state = _qubit_data.find(current_state->first ^ flip); - if (flipped_state == _qubit_data.end()) { // no matching value - if (current_state->first[target]) {// 1 on that qubit - new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M01); - new_qubit_data.emplace(current_state->first, current_state->second * M00); - } - else { - new_qubit_data.emplace(current_state->first, current_state->second * M00); - new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M10); - } - } - // Add up the two values, only when reaching the zero val - else if (!(current_state->first[target])) { - new_state = current_state->second * M00 + flipped_state->second * M01; // zero state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first, new_state); - } - new_state = current_state->second * M10 + flipped_state->second * M00; // one state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first | flip, new_state); - } - } - } else { - new_qubit_data.emplace(current_state->first, current_state->second); - } - } - _qubit_data = std::move(new_qubit_data); - } - } - - void H(logical_qubit_id index){ - // Initialize a new wavefunction, which will store the modified state - // We initialize with twice as much space as the current one, - // as this is the worst case result of an H gate - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); - // This label makes it easier to find associated labels (where the index is flipped) - qubit_label flip(0); - flip.set(index); - // The amplitude for the new state - amplitude new_state; - // Loops over all states in the wavefunction _qubut_date - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - // An iterator pointing to the state labelled by the flip - auto flipped_state = _qubit_data.find(current_state->first ^ flip); - // Checks for whether it needs to add amplitudes from matching states - // or create two new states - if (flipped_state == _qubit_data.end()) { // no matching value - new_qubit_data.emplace(current_state->first & (~flip), current_state->second * _normalizer); - // Flip the value if the second bit, depending on whether the original had 1 or 0 - new_qubit_data.emplace(current_state->first | flip, current_state->second * (current_state->first[index] ? -_normalizer : _normalizer)); - } - else if (!(current_state->first[index])) { - new_state = current_state->second + flipped_state->second; // zero state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first, new_state * _normalizer); - } - - new_state = current_state->second - flipped_state->second; // one state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first | flip, new_state * _normalizer); - } - } - } - // Moves the new data back into the old one (thus destroying - // the old data) - _qubit_data = std::move(new_qubit_data); - } - - void MCH(std::vector const& controls, logical_qubit_id index){ - wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); - qubit_label flip(0); - flip.set(index); - amplitude new_state; - qubit_label checks = _get_mask(controls); - for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { - if ((checks & current_state->first) == checks){ - auto flipped_state = _qubit_data.find(current_state->first ^ flip); - if (flipped_state == _qubit_data.end()) { // no matching value - new_qubit_data.emplace(current_state->first & (~flip), current_state->second * _normalizer); - // Flip the value if the second bit, depending on whether the original had 1 or 0 - new_qubit_data.emplace(current_state->first | flip, current_state->second * (current_state->first[index] ? -_normalizer : _normalizer)); - } - else if (!(current_state->first[index])) { - new_state = current_state->second + flipped_state->second; // zero state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first, new_state * _normalizer); - } - - new_state = current_state->second - flipped_state->second; // one state - if (std::norm(new_state) > _rotation_precision_squared) { - new_qubit_data.emplace(current_state->first | flip, new_state * _normalizer); - } - } - } else { - new_qubit_data.emplace(current_state->first, current_state->second); - } - } - _qubit_data = std::move(new_qubit_data); - } - - // Checks whether a qubit is 0 in all states in the superposition - bool is_qubit_zero(logical_qubit_id target){ - for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end(); ++current_state){ - if (current_state->first[target] && std::norm(current_state->second) > _precision_squared) { - return false; - } - } - return true; - } - - // Creates a new wavefunction hash map indexed by strings - // Not intended for computations but as a way to transfer between - // simulators templated with different numbers of qubits - universal_wavefunction get_universal_wavefunction() { - universal_wavefunction universal_qubit_data = universal_wavefunction(_qubit_data.bucket_count()); - for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end(); ++current_state) { - universal_qubit_data.emplace(current_state->first.to_string(), current_state->second); - } - return universal_qubit_data; - } + // Type for qubit labels, with a specific size built-in + using qubit_label = qubit_label_type; + + // Type of hash maps with the required labels + using wavefunction = abstract_wavefunction; + + QuantumState() { + _qubit_data = wavefunction(); + _qubit_data.max_load_factor(_load_factor); + // Create an initial all-zeros state + _qubit_data.emplace((logical_qubit_id)0, 1); + // Initialize randomness + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution dist(0, 1); + _rng = [gen, dist]() mutable { return dist(gen); }; + + _initialize_threads(); + } + + ~QuantumState() { + complete_threads(); + } + + // Copy data from an existing simulator + // This is used to move between different qubit sizes + // without needing a lot of templated functions + QuantumState(std::shared_ptr old_state) { + // Copy any needed data + _rng = old_state->get_rng(); + // Outputs the previous data with labels as strings + universal_wavefunction old_qubit_data = old_state->get_universal_wavefunction(); + _qubit_data = wavefunction(old_qubit_data.size()); + _load_factor = old_state->get_load_factor(); + _qubit_data.max_load_factor(_load_factor); + // Writes this into the current wavefunction as qubit_label types + for (auto current_state = old_qubit_data.begin(); current_state != old_qubit_data.end(); ++current_state) { + _qubit_data.emplace(qubit_label(current_state->first), current_state->second); + } + // Create local threads for this state + _initialize_threads(); + } + + logical_qubit_id get_num_qubits() { + return (logical_qubit_id)num_qubits; + } + + // Outputs all states and amplitudes to the console + void DumpWavefunction(size_t indent = 0){ + DumpWavefunction(_qubit_data, indent); + } + + // Outputs all states and amplitudes from an input wavefunction to the console + void DumpWavefunction(wavefunction &wfn, size_t indent = 0){ + std::string spacing(indent, ' '); + std::cout << spacing << "Wavefunction:\n"; + auto line_dump = [spacing](qubit_label label, amplitude val){ + std::cout << spacing << " " << label.to_string() << ": "; + std::cout << val.real(); + std::cout << (val.imag() < 0 ? " - " : " + ") << std::abs(val.imag()) << "i\n"; + }; + _DumpWavefunction_base(wfn, line_dump); + std::cout << spacing << "--end wavefunction\n"; + } + + + void set_random_seed(unsigned seed) { + std::mt19937 gen(seed); + std::uniform_real_distribution dist(0, 1); + _rng = [gen, dist]() mutable { return dist(gen); }; + } + + // Used to decide when an amplitude is close enough to 0 to discard + void set_precision(double new_precision) { + _precision = new_precision; + _precision_squared = _precision *_precision; + } + + // Load factor of the underlying hash map + float get_load_factor() { + return _load_factor; + } + + void set_load_factor(float new_load_factor) { + _load_factor = new_load_factor; + } + + // Returns the number of states in superposition + size_t get_wavefunction_size() { + return _qubit_data.size(); + } + + + + // Applies the operator id_coeff*I + pauli_coeff * P + // where P is the Pauli operators defined by axes applied to the qubits in qubits. + void PauliCombination(std::vector const& axes, std::vector const& qubits, amplitude id_coeff, amplitude pauli_coeff) { + // Bit-vectors indexing where gates of each type are applied + qubit_label XYs = 0; + qubit_label YZs = 0; + logical_qubit_id ycount = 0; + for (int i=0; i < axes.size(); i++){ + switch (axes[i]){ + case Gates::Basis::PauliY: + YZs.set(qubits[i]); + XYs.set(qubits[i]); + ycount++; + break; + case Gates::Basis::PauliX: + XYs.set(qubits[i]); + break; + case Gates::Basis::PauliZ: + YZs.set(qubits[i]); + break; + case Gates::Basis::PauliI: + break; + default: + throw std::runtime_error("Bad Pauli basis"); + } + } - // Returns the rng from this simulator - std::function get_rng() { return _rng; } - - // Finishes the threads which wait for the queue - void complete_threads() { - in_use=false; - // If this spuriously wakes up threads, - // they will still end because in_use is false - _current_state = _qubit_data.begin(); - cv.notify_all(); - - for (auto &thread : _thread_pool){ - thread.join(); - } - } - - int get_num_threads() { - #ifndef DOMP_GE_V3 - return static_cast(_thread_pool.size()); - #else - int threads=0; - #pragma omp parallel - { - #pragma omp master - threads = omp_get_num_threads(); - } - return threads; - #endif - } + // All identity + if (XYs.none() && YZs.none()) { + return; + } + + // This branch handles purely Z Pauli vectors + // Purely Z has no addition, which would cause + // problems in the comparison in the next section + if (XYs.none()) { + // 0 terms get the sum of the coefficients + // 1 terms get the difference + pauli_coeff += id_coeff; // id_coeff + pauli_coeff + id_coeff *= 2; + id_coeff -= pauli_coeff; // id_coeff - pauli_coeff + + // To avoid saving states of zero amplitude, these if/else + // check for when one of the coefficients is + // close enough to zero to regard as zero + if (std::norm(pauli_coeff) > _rotation_precision_squared ){ + if (std::norm(id_coeff) > _rotation_precision_squared){ + // If both coefficients are non-zero, we can just modify the state in-place + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + current_state->second *= (poppar(current_state->first & YZs) ? id_coeff : pauli_coeff); + } + } else { + // If id_coeff = 0, then we make a new wavefunction and only add in those that will be multiplied + // by the pauli_coeff + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if (!poppar(current_state->first & YZs)){ + new_qubit_data.emplace(current_state->first, current_state->second * pauli_coeff); + } + } + _qubit_data = std::move(new_qubit_data); + } + } else { + // If pauli_coeff=0, don't add states multiplied by the pauli_coeff + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if (poppar(current_state->first & YZs)){ + new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); + } + } + _qubit_data = std::move(new_qubit_data); + } + } else { // There are some X or Y gates + + // Each Y Pauli adds a global phase of i + switch (ycount % 4) { + case 1: + pauli_coeff *= 1i; + break; + case 2: + pauli_coeff *= -1; + break; + case 3: + pauli_coeff *= -1i; + break; + default: + break; + } + // When both the state and flipped state are in superposition, when adding the contribution of + // the flipped state, we add phase depending on the 1s in the flipped state + // This phase would be the parity of (flipped_state->first ^ YZs) + // However, we know that flipped_state->first = current_state->first ^ YXs + // So the parity of the flipped state will be the parity of the current state, plus + // the parity of YZs & YXs, i.e., the parity of the number of Ys + // Since this is constant for all states, we compute it once here and save it + // Then we only compute the parity of the current state + amplitude pauli_coeff_alt = ycount % 2 ? -pauli_coeff : pauli_coeff; + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); + amplitude new_state; + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + auto alt_state = _qubit_data.find(current_state->first ^ XYs); + if (alt_state == _qubit_data.end()) { // no matching value + new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); + new_qubit_data.emplace(current_state->first ^ XYs, current_state->second * (poppar(current_state->first & YZs) ? -pauli_coeff : pauli_coeff)); + } + else if (current_state->first < alt_state->first) { + // Each Y and Z gate adds a phase (since Y=iXZ) + bool parity = poppar(current_state->first & YZs); + new_state = current_state->second * id_coeff + alt_state->second * (parity ? -pauli_coeff_alt : pauli_coeff_alt); + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first, new_state); + } + + new_state = alt_state->second * id_coeff + current_state->second * (parity ? -pauli_coeff : pauli_coeff); + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(alt_state->first, new_state); + } + } + } + _qubit_data = std::move(new_qubit_data); + } + } + + // Applies the operator id_coeff*I + pauli_coeff * P + // where P is the Pauli operators defined by axes applied to the qubits in qubits. + // Controlled version + void MCPauliCombination(std::vector const& controls, std::vector const& axes, std::vector const& qubits, amplitude id_coeff, amplitude pauli_coeff) { + // Bit-vectors indexing where gates of each type are applied + qubit_label cmask = _get_mask(controls); + qubit_label XYs = 0; + qubit_label YZs = 0; + logical_qubit_id ycount = 0; + // Used for comparing pairs + logical_qubit_id any_xy = -1; + for (int i=0; i < axes.size(); i++){ + switch (axes[i]){ + case Gates::Basis::PauliY: + YZs.set(qubits[i]); + XYs.set(qubits[i]); + ycount++; + any_xy = qubits[i]; + break; + case Gates::Basis::PauliX: + XYs.set(qubits[i]); + any_xy = qubits[i]; + break; + case Gates::Basis::PauliZ: + YZs.set(qubits[i]); + break; + case Gates::Basis::PauliI: + break; + default: + throw std::runtime_error("Bad Pauli basis"); + } + } + + // This branch handles purely Z Pauli vectors + // Purely Z has no addition, which would cause + // problems in the comparison in the next section + if (XYs.none()) { + // 0 terms get the sum of the coefficients + // 1 terms get the difference + pauli_coeff += id_coeff; // <- id_coeff + pauli_coeff + id_coeff *= 2; + id_coeff -= pauli_coeff; // <- id_coeff - pauli_coeff + + // To avoid saving states of zero amplitude, these if/else + // check for when one of the coefficients is + // close enough to zero to regard as zero + if (std::norm(pauli_coeff) > _rotation_precision_squared ){ + if (std::norm(id_coeff) > _rotation_precision_squared){ + // If both coefficients are non-zero, we can just modify the state in-place + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if ((current_state->first & cmask)==cmask) { + current_state->second *= (poppar(current_state->first & YZs) ? id_coeff : pauli_coeff); + } + } + } else { + // If id_coeff = 0, then we make a new wavefunction and only add in those that will be multiplied + // by the pauli_coeff + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if (!poppar(current_state->first & YZs) && (current_state->first & cmask)==cmask){ + new_qubit_data.emplace(current_state->first, current_state->second * pauli_coeff); + } + } + _qubit_data = std::move(new_qubit_data); + } + } else { + // If pauli_coeff=0, don't add states multiplied by the pauli_coeff + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size()); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if (poppar(current_state->first & YZs) && (current_state->first & cmask)==cmask){ + new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); + } + } + _qubit_data = std::move(new_qubit_data); + } + } else { // There are some X or Y gates + // Each Y Pauli adds a global phase of i + switch (ycount % 4) { + case 1: + pauli_coeff *= 1i; + break; + case 2: + pauli_coeff *= -1; + break; + case 3: + pauli_coeff *= -1i; + break; + default: + break; + } + // When both the state and flipped state are in superposition, when adding the contribution of + // the flipped state, we add phase depending on the 1s in the flipped state + // This phase would be the parity of (flipped_state->first ^ YZs) + // However, we know that flipped_state->first = current_state->first ^ YXs + // So the parity of the flipped state will be the parity of the current state, plus + // the parity of YZs & YXs, i.e., the parity of the number of Ys + // Since this is constant for all states, we compute it once here and save it + // Then we only compute the parity of the current state + amplitude pauli_coeff_alt = ycount % 2 ? -pauli_coeff : pauli_coeff; + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); + amplitude new_state; + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if ((current_state->first & cmask)==cmask) { + auto alt_state = _qubit_data.find(current_state->first ^ XYs); + if (alt_state == _qubit_data.end()) { // no matching value + new_qubit_data.emplace(current_state->first, current_state->second * id_coeff); + new_qubit_data.emplace(current_state->first ^ XYs, current_state->second * (poppar(current_state->first & YZs) ? -pauli_coeff : pauli_coeff)); + } + else if (current_state->first < alt_state->first) { //current_state->first[any_xy]){// + // Each Y and Z gate adds a phase (since Y=iXZ) + bool parity = poppar(current_state->first & YZs); + new_state = current_state->second * id_coeff + alt_state->second * (parity ? -pauli_coeff_alt : pauli_coeff_alt); + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first, new_state); + } + + new_state = alt_state->second * id_coeff + current_state->second * (parity ? -pauli_coeff : pauli_coeff); + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(alt_state->first, new_state); + } + } + } else { + new_qubit_data.emplace(current_state->first, current_state->second); + } + } + _qubit_data = std::move(new_qubit_data); + } + } + + + bool M(logical_qubit_id target) { + qubit_label flip = qubit_label(); + flip.set(target); + + bool result = _qubit_data.begin()->first[target]; + + double zero_probability = 0.0; + double one_probability = 0.0; + + // Writes data into a ones or zeros wavefunction + // as it adds up probability + // Once it's finished, it picks one randomly, normalizes + // then keeps that one as the new wavefunction + wavefunction ones = make_wavefunction(_qubit_data.size()/2); + wavefunction zeros = make_wavefunction(_qubit_data.size()/2); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + double square_amplitude = std::norm(current_state->second); + if (current_state->first[target]) { + one_probability += square_amplitude; + ones.emplace(current_state->first, current_state->second); + } + else { + zero_probability += square_amplitude; + zeros.emplace(current_state->first, current_state->second); + } + } + // Randomly select + result = (_rng() <= one_probability); + + wavefunction &new_qubit_data = result ? ones : zeros; + // Create a new, normalized state + double normalizer = 1.0/std::sqrt((result) ? one_probability : zero_probability); + for (auto current_state = (new_qubit_data).begin(); current_state != (new_qubit_data).end(); ++current_state) { + current_state->second *= normalizer; + } + _qubit_data = std::move(new_qubit_data); + + return result; + } + + void Reset(logical_qubit_id target) { + qubit_label flip = qubit_label(0); + flip.set(target); + + double zero_probability = 0.0; + double one_probability = 0.0; + + // Writes data into a ones or zeros wavefunction + // as it adds up probability + // Once it's finished, it picks one randomly, normalizes + // then keeps that one as the new wavefunction + + // Used to set the qubit to 0 in the measured result + qubit_label new_mask = qubit_label(); + new_mask.set(); // sets all bits to 1 + new_mask.set(target, 0); + wavefunction ones = make_wavefunction(_qubit_data.size()/2); + wavefunction zeros = make_wavefunction(_qubit_data.size()/2); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + double square_amplitude = std::norm(current_state->second); + if (current_state->first[target]) { + one_probability += square_amplitude; + ones.emplace(current_state->first & new_mask, current_state->second); + } + else { + zero_probability += square_amplitude; + zeros.emplace(current_state->first & new_mask, current_state->second); + } + } + // Randomly select + bool result = (_rng() <= one_probability); + + wavefunction &new_qubit_data = result ? ones : zeros; + // Create a new, normalized state + double normalizer = 1.0/std::sqrt((result) ? one_probability : zero_probability); + for (auto current_state = (new_qubit_data).begin(); current_state != (new_qubit_data).end(); ++current_state) { + current_state->second *= normalizer; + } + _qubit_data = std::move(new_qubit_data); + } + + + // Samples a state from the superposition with probably proportion to + // the amplitude, returning a string of the bits of that state. + // Unlike measurement, this does not modify the state + std::string Sample() { + double probability = _rng(); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + double square_amplitude = std::norm(current_state->second); + probability -= square_amplitude; + if (probability <= 0){ + return current_state->first.to_string(); + } + } + return _qubit_data.begin()->first.to_string(); + } + + void Assert(std::vector const& axes, std::vector const& qubits, bool result) { + // Bit-vectors indexing where gates of each type are applied + qubit_label XYs = 0; + qubit_label YZs = 0; + logical_qubit_id ycount = 0; + for (int i=0; i < axes.size(); i++){ + switch (axes[i]){ + case Gates::Basis::PauliY: + YZs.set(qubits[i]); + XYs.set(qubits[i]); + ycount++; + break; + case Gates::Basis::PauliX: + XYs.set(qubits[i]); + break; + case Gates::Basis::PauliZ: + YZs.set(qubits[i]); + break; + case Gates::Basis::PauliI: + break; + default: + throw std::runtime_error("Bad Pauli basis"); + } + } + + amplitude phaseShift = result ? -1 : 1; + // Each Y Pauli adds a global phase of i + switch (ycount % 4) { + case 1: + phaseShift *= 1i; + break; + case 2: + phaseShift *= -1; + break; + case 3: + phaseShift *= -1i; + break; + default: + break; + } + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if (std::norm(_qubit_data.find(current_state->first ^ XYs)->second - current_state->second * (poppar(current_state->first & YZs) ? -phaseShift : phaseShift)) > _precision_squared) { + qubit_label label = current_state->first; + amplitude val = current_state->second; + std::cout << "Problematic state: " << label << "\n"; + std::cout << "Expected " << val * (poppar(label & YZs) ? -phaseShift : phaseShift); + std::cout << ", got " << _qubit_data.find(label ^ XYs)->second << "\n"; + std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; + throw std::runtime_error("Not an eigenstate"); + } + } + } + + // Returns the probability of a given measurement in a Pauli basis + // by decomposing each pair of computational basis states into eigenvectors + // and adding the coefficients of the respective components + double MeasurementProbability(std::vector const& axes, std::vector const& qubits) { + // Bit-vectors indexing where gates of each type are applied + qubit_label XYs = 0; + qubit_label YZs = 0; + logical_qubit_id ycount = 0; + for (int i=0; i < axes.size(); i++){ + switch (axes[i]){ + case Gates::Basis::PauliY: + YZs.set(qubits[i]); + XYs.set(qubits[i]); + ycount++; + break; + case Gates::Basis::PauliX: + XYs.set(qubits[i]); + break; + case Gates::Basis::PauliZ: + YZs.set(qubits[i]); + break; + case Gates::Basis::PauliI: + break; + default: + throw std::runtime_error("Bad Pauli basis"); + } + } + amplitude phaseShift = 1; + + // Each Y Pauli adds a global phase of i + switch (ycount % 4) { + case 1: + phaseShift *= amplitude(0, 1); + break; + case 2: + phaseShift *= -1; + break; + case 3: + phaseShift *= amplitude(0, -1); + break; + default: + break; + } + // Let P be the pauli operation, |psi> the state + // projection = + + // _qubit_data represents |psi> as sum_x a_x |x>, + // where all |x> are orthonormal. Thus, the projection + // will be the product of a_x and a_P(x), where P|x>=|P(x)> + // Thus, for each |x>, we compute P(x) and look for that state + // If there is a match, we add the product of their coefficients + // to the projection, times a phase dependent on how many Ys and Zs match + // the 1 bits of x + amplitude projection = 0.0; + auto flipped_state = _qubit_data.end(); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + flipped_state = _qubit_data.find(current_state->first ^ XYs); // no match returns _qubit_data.end() + projection += current_state->second * (flipped_state == _qubit_data.end() ? 0 : std::conj(flipped_state->second)) * (poppar(current_state->first & YZs) ? -phaseShift : phaseShift); + } + // The projector onto the -1 eigenspace (a result of "One") is 0.5 * (I - P) + // So = 0.5 - 0.5* + // should always be real so this only takes the real part + return 0.5 - 0.5 * projection.real(); + } + + bool Measure(std::vector const& axes, std::vector const& qubits){ + // Find a probability to get a specific result + double probability = MeasurementProbability(axes, qubits); + bool result = _rng() <= probability; + probability = std::sqrt(probability); + // This step executes immediately so that we reduce the number of states in superposition + PauliCombination(axes, qubits, 0.5/probability, (result ? -0.5 : 0.5)/probability); + return result; + } + + + // Probe the amplitude of a single basis state + amplitude probe(qubit_label label) { + auto qubit = _qubit_data.find(label); + // States not in the hash map are assumed to be 0 + if (qubit == _qubit_data.end()) { + return amplitude(0.0, 0.0); + } + else { + return qubit->second; + } + } + + amplitude probe(std::string label) { + qubit_label bit_label = qubit_label(label); + return probe(bit_label); + } + + // Dumps the state of a subspace of particular qubits, if they are not entangled + // This requires it to detect if the subspace is entangled, construct a new + // projected wavefunction, then call the `callback` function on each state. + bool dump_qubits(std::vector qubits, void (*callback)(char*, double, double)) { + // Create two wavefunctions + // check if they are tensor products + wavefunction dump_wfn; + wavefunction leftover_wfn; + if (!_split_wavefunction(_get_mask(qubits), dump_wfn, leftover_wfn)){ + return false; + } else { + _DumpWavefunction_base(dump_wfn, [qubits, callback](qubit_label label, amplitude val){ + std::string label_string(qubits.size(), '0'); + for (size_t i=0; i < qubits.size(); i++){ + label_string[i] = label[qubits[i]] ? '1' : '0'; + } + callback(const_cast(label_string.c_str()), val.real(), val.imag()); + }); + return true; + } + } + + // Dumps all the states in superposition via a callback function + void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) { + _DumpWavefunction_base(_qubit_data, [max_qubit_id, callback](qubit_label label, amplitude val){ + callback(const_cast(label.to_string().substr(num_qubits - 1 - max_qubit_id, max_qubit_id + 1).c_str()), val.real(), val.imag()); + }); + } + + // Execute a queue of phase/permutation gates + void phase_and_permute(std::list const &operation_list){ + if (operation_list.size()==0){return;} + + // Condense the list into a memory-efficient vector with qubit labels + _operation_vector.reserve(operation_list.size()); + for (auto op : operation_list){ + switch (op.gate_type) { + case OP::X: + case OP::Y: + case OP::Z: + _operation_vector.push_back(internal_operation(op.gate_type, op.target)); + break; + case OP::MCX: + case OP::MCY: + _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls))); + break; + case OP::MCZ: + _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls).set(op.target))); + break; + case OP::Phase: + _operation_vector.push_back(internal_operation(op.gate_type, op.target, op.phase)); + break; + case OP::MCPhase: + _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls).set(op.target), op.phase)); + break; + case OP::SWAP: + _operation_vector.push_back(internal_operation(op.gate_type, op.target, op.target_2)); + break; + case OP::MCSWAP: + _operation_vector.push_back(internal_operation(op.gate_type, op.target, _get_mask(op.controls), op.target_2)); + break; + case OP::Assert: + _operation_vector.push_back(internal_operation(op.gate_type, _get_mask(op.controls), op.result)); + break; + default: + throw std::runtime_error("Unsupported operation"); + break; + } + } + + _new_qubit_data = make_wavefunction(); + + // Threading introduces a lot of overhead so it + // is only worthwhile for large instances + // 64 and 4096 are empirical values + if (_operation_vector.size() < 64 || _qubit_data.size() < 4096){ // small; do not thread + // Iterates through and applies all operations + for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end(); ++current_state){ + qubit_label label = current_state->first; + amplitude val = current_state->second; + // Iterate through vector of operations and apply each gate + for (int i=0; i < _operation_vector.size(); i++) { + auto &op = _operation_vector[i]; + switch (op.gate_type) { + case OP::X: + label.flip(op.target); + break; + case OP::MCX: + if ((op.controls & label) == op.controls){ + label.flip(op.target); + } + break; + case OP::Y: + label.flip(op.target); + val *= (label[op.target]) ? 1i : -1i; + break; + case OP::MCY: + if ((op.controls & label) == op.controls){ + label.flip(op.target); + val *= (label[op.target]) ? 1i : -1i; + } + break; + case OP::Z: + val *= (label[op.target] ? -1 : 1); + break; + case OP::MCZ: + val *= ((op.controls & label) == op.controls) ? -1 : 1; + break; + case OP::Phase: + val *= label[op.target] ? op.phase : 1; + break; + case OP::MCPhase: + val *= ((op.controls & label) == op.controls) ? op.phase : 1; + break; + case OP::SWAP: + if (label[op.target] != label[op.target_2]){ + label.flip(op.target); + label.flip(op.target_2); + } + break; + case OP::MCSWAP: + if (((label & op.controls) == op.controls) && (label[op.target] != label[op.target_2])){ + label.flip(op.target); + label.flip(op.target_2); + } + break; + case OP::Assert: + if (poppar(label & op.controls) != op.result && std::norm(val) > _precision_squared){ + std::cout << "Problematic state: " << label << "\n"; + std::cout << "Amplitude: " << val << "\n"; + std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; + throw std::runtime_error("Assert failed"); + } + break; + default: + throw std::runtime_error("Unsupported operation"); + break; + } + } + // Insert the new state into the new wavefunction + _new_qubit_data.emplace(label, val); + } + _qubit_data = std::move(_new_qubit_data); + _new_qubit_data.clear(); + _operation_vector.clear(); + } else { // Large enough to multi-thread + #ifndef DOMP_GE_V3 // OMP version is too low; uses condition variables + // Lock the mutex so the threads stay asleep during prep + std::unique_lock state_lock(_state_mtx); + _new_qubit_data = wavefunction(_qubit_data.size()); + // jump_size gives a rough guess for how many states each thread should process + // to minimize contention for current_state + _jump_size = std::max((size_t)1 , (size_t)(_qubit_data.size() / (8*_thread_pool.size()))); + // Set the _current state (this allows the threads to wake up) + _current_state = _qubit_data.begin(); + + // Wake up all threads + cv.notify_all(); + // Wait for the number of running threads to be 0 + // This will spuriously wake many times + cv.wait(state_lock, [&](){ + return _current_state == _qubit_data.end() && _running_threads == 0; + }); + // Here all threads are finished + #else + #pragma omp parallel + { + _jump_size = std::max((size_t)1 , (size_t)(_qubit_data.size() / (8*omp_get_num_threads()))); + #pragma omp single + { + auto internal_state = _qubit_data.begin(); + auto internal_end = _qubit_data.end(); + for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end();){ + internal_state = current_state; + current_state.jump_forward(_jump_size); + internal_end = current_state; + #pragma omp task firstprivate(internal_state) firstprivate(internal_end) + { + qubit_label label; + amplitude val; + for (; internal_state != internal_end; ++internal_state){ + label = internal_state->first; + val = internal_state->second; + for (int i=0; i < _operation_vector.size(); i++) { + auto &op = _operation_vector[i]; + switch (op.gate_type) { + case OP::X: + label.flip(op.target); + break; + case OP::MCX: + if ((op.controls & label) == op.controls){ + label.flip(op.target); + } + break; + case OP::Y: + label.flip(op.target); + val *= (label[op.target]) ? 1i : -1i; + break; + case OP::MCY: + if ((op.controls & label) == op.controls){ + label.flip(op.target); + val *= (label[op.target]) ? 1i : -1i; + } + break; + case OP::Z: + val *= (label[op.target] ? -1 : 1); + break; + case OP::MCZ: + val *= ((op.controls & label) == op.controls) ? -1 : 1; + break; + case OP::Phase: + val *= label[op.target] ? op.phase : 1; + break; + case OP::MCPhase: + val *= ((op.controls & label) == op.controls) ? op.phase : 1; + break; + case OP::SWAP: + if (label[op.target] != label[op.target_2]){ + label.flip(op.target); + label.flip(op.target_2); + } + break; + case OP::MCSWAP: + if (((label & op.controls) == op.controls) && (label[op.target] != label[op.target_2])){ + label.flip(op.target); + label.flip(op.target_2); + } + break; + case OP::Assert: + if (poppar(label & op.controls) != op.result && std::norm(val) > _precision_squared){ + std::cout << "Problematic state: " << label << "\n"; + std::cout << "Amplitude: " << val << "\n"; + std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; + throw std::runtime_error("Assert failed"); + } + break; + default: + throw std::runtime_error("Unsupported operation"); + break; + } + } + #pragma omp critical + { + _new_qubit_data.emplace(label, val); + } + } + } + + } + } + #pragma omp barrier + } + #endif + _qubit_data = std::move(_new_qubit_data); + _operation_vector.clear(); + _new_qubit_data.clear(); + } + } + + void R(Gates::Basis b, double phi, logical_qubit_id index){ + // Z rotation can be done in-place + if (b == Gates::Basis::PauliZ) { + amplitude exp_0 = amplitude(std::cos(phi / 2.0), -std::sin(phi / 2.0)); + amplitude exp_1 = amplitude(std::cos(phi / 2.0), std::sin(phi / 2.0)); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + current_state->second *= current_state->first[index] ? exp_1 : exp_0; + } + } + else if (b == Gates::Basis::PauliX || b == Gates::Basis::PauliY) { + amplitude M00 = std::cos(phi / 2.0); + amplitude M01 = -std::sin(phi / 2.0) * (b == Gates::Basis::PauliY ? 1 : 1i); + if (std::norm(M00) < _rotation_precision_squared){ + // This is just a Y or X gate + phase_and_permute(std::list{operation(b==Gates::Basis::PauliY ? OP::Y : OP::X, index)}); + return; + } else if (std::norm(M01) < _rotation_precision_squared){ + // just an identity + return; + } + + amplitude M10 = M01 * amplitude(b == Gates::Basis::PauliY ? -1 : 1); + // Holds the amplitude of the new state to make it easier to check if it's non-zero + amplitude new_state; + qubit_label flip(0); + flip.set(index); + wavefunction new_qubit_data = make_wavefunction(); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + auto flipped_state = _qubit_data.find(current_state->first ^ flip); + if (flipped_state == _qubit_data.end()) { // no matching value + if (current_state->first[index]) {// 1 on that qubit + new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M01); + new_qubit_data.emplace(current_state->first, current_state->second * M00); + } + else { + new_qubit_data.emplace(current_state->first, current_state->second * M00); + new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M10); + } + } + // Add up the two values, only when reaching the zero value + else if (!(current_state->first[index])) { + new_state = current_state->second * M00 + flipped_state->second * M01; // zero state + if (std::norm(new_state) > _rotation_precision_squared ) { + new_qubit_data.emplace(current_state->first, new_state); + } + new_state = current_state->second * M10 + flipped_state->second * M00; // one state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(flipped_state->first, new_state); + } + } + } + _qubit_data = std::move(new_qubit_data); + } + } + + // Multi-controlled rotation + void MCR (std::vector const& controls, Gates::Basis b, double phi, logical_qubit_id target) { + qubit_label checks = _get_mask(controls); + // A Z-rotation can be done without recreating the wavefunction + if (b == Gates::Basis::PauliZ) { + amplitude exp_0 = amplitude(std::cos(phi / 2.0), -std::sin(phi / 2.0)); + amplitude exp_1 = amplitude(std::cos(phi / 2.0), std::sin(phi / 2.0)); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if ((current_state->first & checks)==checks){ + current_state->second *= current_state->first[target] ? exp_1 : exp_0; + } + } + } + // X or Y requires a new wavefunction + else if (b == Gates::Basis::PauliX || b == Gates::Basis::PauliY) { + amplitude M00 = std::cos(phi / 2.0); + amplitude M01 = -std::sin(phi / 2.0) * (b == Gates::Basis::PauliY ? 1 : 1i); + amplitude M10 = (b == Gates::Basis::PauliY ? -1.0 : 1.0) * M01; + + if (std::norm(M00) < _rotation_precision_squared){ + // This is just an MCY or MCX gate, but with a phase + // So we need to preprocess with a multi-controlled phase + if (b==Gates::Basis::PauliY){ + amplitude phase = -1i/M01; + phase_and_permute(std::list{ + operation(OP::MCPhase, controls[0], controls, phase), + operation(OP::MCY, target, controls) + }); + } else { + amplitude phase = 1.0/M01; + phase_and_permute(std::list{ + operation(OP::MCPhase, controls[0], controls, phase), + operation(OP::MCY, target, controls) + }); + } + return; + } else if (std::norm(M01) < _rotation_precision_squared){ + // This is equivalent to a multi-controlled Z if the rotation is -1 + if (std::norm(M01 + 1.0) < _rotation_precision_squared){ + phase_and_permute(std::list{operation(OP::MCZ, controls[0], controls)}); + } + return; + } + + amplitude new_state; + qubit_label flip(0); + flip.set(target); + wavefunction new_qubit_data = make_wavefunction(); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if ((current_state->first & checks)==checks){ + auto flipped_state = _qubit_data.find(current_state->first ^ flip); + if (flipped_state == _qubit_data.end()) { // no matching value + if (current_state->first[target]) {// 1 on that qubit + new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M01); + new_qubit_data.emplace(current_state->first, current_state->second * M00); + } + else { + new_qubit_data.emplace(current_state->first, current_state->second * M00); + new_qubit_data.emplace(current_state->first ^ flip, current_state->second * M10); + } + } + // Add up the two values, only when reaching the zero val + else if (!(current_state->first[target])) { + new_state = current_state->second * M00 + flipped_state->second * M01; // zero state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first, new_state); + } + new_state = current_state->second * M10 + flipped_state->second * M00; // one state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first | flip, new_state); + } + } + } else { + new_qubit_data.emplace(current_state->first, current_state->second); + } + } + _qubit_data = std::move(new_qubit_data); + } + } + + void H(logical_qubit_id index){ + // Initialize a new wavefunction, which will store the modified state + // We initialize with twice as much space as the current one, + // as this is the worst case result of an H gate + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); + // This label makes it easier to find associated labels (where the index is flipped) + qubit_label flip(0); + flip.set(index); + // The amplitude for the new state + amplitude new_state; + // Loops over all states in the wavefunction _qubut_date + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + // An iterator pointing to the state labelled by the flip + auto flipped_state = _qubit_data.find(current_state->first ^ flip); + // Checks for whether it needs to add amplitudes from matching states + // or create two new states + if (flipped_state == _qubit_data.end()) { // no matching value + new_qubit_data.emplace(current_state->first & (~flip), current_state->second * _normalizer); + // Flip the value if the second bit, depending on whether the original had 1 or 0 + new_qubit_data.emplace(current_state->first | flip, current_state->second * (current_state->first[index] ? -_normalizer : _normalizer)); + } + else if (!(current_state->first[index])) { + new_state = current_state->second + flipped_state->second; // zero state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first, new_state * _normalizer); + } + + new_state = current_state->second - flipped_state->second; // one state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first | flip, new_state * _normalizer); + } + } + } + // Moves the new data back into the old one (thus destroying + // the old data) + _qubit_data = std::move(new_qubit_data); + } + + void MCH(std::vector const& controls, logical_qubit_id index){ + wavefunction new_qubit_data = make_wavefunction(_qubit_data.size() * 2); + qubit_label flip(0); + flip.set(index); + amplitude new_state; + qubit_label checks = _get_mask(controls); + for (auto current_state = (_qubit_data).begin(); current_state != (_qubit_data).end(); ++current_state) { + if ((checks & current_state->first) == checks){ + auto flipped_state = _qubit_data.find(current_state->first ^ flip); + if (flipped_state == _qubit_data.end()) { // no matching value + new_qubit_data.emplace(current_state->first & (~flip), current_state->second * _normalizer); + // Flip the value if the second bit, depending on whether the original had 1 or 0 + new_qubit_data.emplace(current_state->first | flip, current_state->second * (current_state->first[index] ? -_normalizer : _normalizer)); + } + else if (!(current_state->first[index])) { + new_state = current_state->second + flipped_state->second; // zero state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first, new_state * _normalizer); + } + + new_state = current_state->second - flipped_state->second; // one state + if (std::norm(new_state) > _rotation_precision_squared) { + new_qubit_data.emplace(current_state->first | flip, new_state * _normalizer); + } + } + } else { + new_qubit_data.emplace(current_state->first, current_state->second); + } + } + _qubit_data = std::move(new_qubit_data); + } + + // Checks whether a qubit is 0 in all states in the superposition + bool is_qubit_zero(logical_qubit_id target){ + for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end(); ++current_state){ + if (current_state->first[target] && std::norm(current_state->second) > _precision_squared) { + return false; + } + } + return true; + } + + // Creates a new wavefunction hash map indexed by strings + // Not intended for computations but as a way to transfer between + // simulators templated with different numbers of qubits + universal_wavefunction get_universal_wavefunction() { + universal_wavefunction universal_qubit_data = universal_wavefunction(_qubit_data.bucket_count()); + for (auto current_state = _qubit_data.begin(); current_state != _qubit_data.end(); ++current_state) { + universal_qubit_data.emplace(current_state->first.to_string(), current_state->second); + } + return universal_qubit_data; + } + + // Returns the rng from this simulator + std::function get_rng() { return _rng; } + + // Finishes the threads which wait for the queue + void complete_threads() { + in_use=false; + // If this spuriously wakes up threads, + // they will still end because in_use is false + _current_state = _qubit_data.begin(); + cv.notify_all(); + + for (auto &thread : _thread_pool){ + thread.join(); + } + } + + int get_num_threads() { + #ifndef DOMP_GE_V3 + return static_cast(_thread_pool.size()); + #else + int threads=0; + #pragma omp parallel + { + #pragma omp master + threads = omp_get_num_threads(); + } + return threads; + #endif + } private: - // Internal type used to store operations with bitsets - // instead of vectors of qubit ids - using internal_operation = condensed_operation; - - // Hash table of the wavefunction - wavefunction _qubit_data; - - // Internal random numbers - std::function _rng; - - // Threshold to assert that something is zero when asserting it is 0 - double _precision = 1e-5; - // Threshold at which something is zero when - // deciding whether to add it into the superposition - double _rotation_precision = 1e-6; - // Often we compare to norms, so the precision must be squared - double _precision_squared = _precision * _precision; - double _rotation_precision_squared = _rotation_precision*_rotation_precision; - - // Normalizer for H and T gates (1/sqrt(2) as an amplitude) - const amplitude _normalizer = amplitude(1.0, 0.0) / std::sqrt(2.0); - - // The default for bytell_hash_map - // Used when allocating new wavefunctions - float _load_factor = 0.9375; - - // Makes a wavefunction that is preallocated to the right size - // and has the correct load factor - wavefunction make_wavefunction() { - wavefunction data((size_t)(_qubit_data.size() / _load_factor)); - data.max_load_factor(_load_factor); - return data; - } - wavefunction make_wavefunction(uint64_t n_states) { - wavefunction data((size_t)(n_states / _load_factor)); - data.max_load_factor(_load_factor); - return data; - } - - // Creates a qubit_label as a bit mask from a set of indices - qubit_label _get_mask(std::vector const& indices){ - return get_mask(indices); - } - - //********* Member variables for multithreading permutations ******// - // Shared iterator through the wavefunction - // Should point to _qubit_data.end() unless in use - decltype(_qubit_data.begin()) _current_state; - // Vector of waiting threads - std::vector _thread_pool; - // Condition variable, locked on _state_mtx - std::condition_variable cv; - // Mutex to prevent concurrent access to _current_state - std::mutex _state_mtx; - // Mutex to prevent concurrent access to _new_qubit_data - std::mutex _wfn_mtx; - // Used as a flag to decide when the thread pool is necessary, - // i.e., it is only set to false when the simulator is being destroyed - std::atomic in_use = true; - // Counts running threads, so controlling thread knows when the permutation is finished - std::atomic _running_threads = 0; - // A vector of operations that all threads point to to read which operations to apply - std::vector _operation_vector; - // New wavefunction which will replace the old one - wavefunction _new_qubit_data; - // The number of elements of the hash map that each thread will handle in one iteration - size_t _jump_size; - // Max number of threads - int max_num_threads = std::thread::hardware_concurrency(); - - // Called on creation, starts the thread pool and waits for all threads to - // finish initialization tasks - void _initialize_threads() { - // Unnecessary if OpenMP Is available - #ifndef DOMP_GE_V3 - // Prevents spurious wake-ups - _current_state = _qubit_data.end(); - _thread_pool = std::vector(); - #ifndef _OPENMP - int num_threads = std::thread::hardware_concurrency(); - #else - int num_threads = 1; - #pragma omp parallel - { - #pragma omp single - num_threads = omp_get_num_threads(); - } - #endif - // Lock the state, so that all the new threads get stopped during their execution - std::unique_lock state_lock(_state_mtx); - for (int i = 0; i < num_threads; i++){ - // Each running thread will decrement this just before it waits on the condition variable - ++_running_threads; - _thread_pool.push_back(std::thread([this](){this->_wait_to_permute();})); - } - // Waits until all the running threads have finished initializing - // This ensures the simulation does not actually start - // until the thread pool has finished initializing - cv.wait(state_lock, [&](){return _running_threads==0;}); - #endif - } - - // Function that all threads in _thread_pool run - // Waits on the condition variable, then applies _operation_vector - // to the qubit data - void _wait_to_permute(){ - // Prep local variables - auto local_state = _qubit_data.begin(); - auto local_end = _qubit_data.end(); - // Initialize a wfn lock to avoid recreating this object in every loop - std::unique_lock wfn_lock(_wfn_mtx); - wfn_lock.unlock(); - // First lock the state so we can wait on it - std::unique_lock state_lock(_state_mtx); - // Set the _current_state to ensure it waits - _current_state = _qubit_data.end(); - // Reduce number of running threads - _running_threads--; - // Notify all (i.e., the constructor thread) - // which will wait until the lock is released to check the number of running threads - cv.notify_all(); - // Wait until we set the current state to the start of the wavefunction - cv.wait(state_lock, [&](){ - return _current_state != _qubit_data.end(); - }); - - // Loop continuously - while (in_use){ - ++_running_threads; - // Loop through the wavefunction - while (_current_state != _qubit_data.end()){ - // Update local variables - local_state = _current_state; - _current_state.jump_forward(_jump_size); - local_end = _current_state; - // Unlock state to allow other threads to modify their state - state_lock.unlock(); - // Loop through the chunk that this thread has taken - for (; local_state != local_end; ++local_state){ - qubit_label label = local_state->first; - amplitude val = local_state->second; - - for (auto op : _operation_vector){ - switch (op.gate_type) { - case OP::X: - label.flip(op.target); - break; - case OP::MCX: - if ((op.controls & label) == op.controls){ - label.flip(op.target); - } - break; - case OP::Y: - label.flip(op.target); - val *= (label[op.target]) ? 1i : -1i; - break; - case OP::MCY: - if ((op.controls & label) == op.controls){ - label.flip(op.target); - val *= (label[op.target]) ? 1i : -1i; - } - break; - case OP::Z: - val *= (label[op.target] ? -1 : 1); - break; - case OP::MCZ: - val *= ((op.controls & label) == op.controls) ? -1 : 1; - break; - case OP::Phase: - val *= label[op.target] ? op.phase : 1; - break; - case OP::MCPhase: - val *= ((op.controls & label) == op.controls) ? op.phase : 1; - break; - case OP::SWAP: - if (label[op.target] != label[op.target_2]){ - label.flip(op.target); - label.flip(op.target_2); - } - break; - case OP::MCSWAP: - if (((label & op.controls) == op.controls) && (label[op.target] != label[op.target_2])){ - label.flip(op.target); - label.flip(op.target_2); - } - break; - case OP::Assert: - if (poppar(label & op.controls) != op.result && std::norm(val) > _precision_squared){ - std::cout << "Problematic state: " << label << "\n"; - std::cout << "Amplitude: " << val << "\n"; - std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; - throw std::runtime_error("Assert failed"); - } - break; - default: - throw std::runtime_error("Unsupported operation"); - break; - } - } - - wfn_lock.lock(); - _new_qubit_data.emplace(label, val); - wfn_lock.unlock(); - } - // Lock before checking and modifying the current state - state_lock.lock(); - } - // This thread has finished iterating through the wavefunction - --_running_threads; - // Notify all other threads (i.e., the controlling thread) - // if this is the last thread to finish - // The check on _running_threads avoids spurious wake-ups to - // the controlling thread, but takes a bit of time as _running_threads - // is atomic - if (_running_threads ==0) - cv.notify_all(); - // Wait on the state lock before repeating the loop - // The wait unlocks state_lock - cv.wait(state_lock, [&](){return _current_state != _qubit_data.end();}); - } - } - - - // Split the wavefunction if separable, otherwise return false - // Idea is that if we have a_bb|b1>|b2> as the first state, then for - // any other state a_xx|x1>|x2>, we must also have a_xb|x1>|b2> and a_bx|b1>|x2> - // in superposition. - // Also, the coefficients must separate as a_bb=c_b*d_b and a_xx = c_x*d_x, implying - // that a_xb = c_x*d_b and a_bx = c_b * d_x, and thus we can check this holds if - // a_bb*a_xx = a_bx * a_xb. - // If this holds: we write (a_xx/a_bx)|x1> into the first wavefunction and (a_xx/a_xb)|x2> - // into the second. - // This means the coefficients of wfn1 are all of the form (c_x/c_b); - // Thus, the norm of wfn1 will be 1/|c_b|^2; thus the norm of wfn2 is 1/|d_b|^2 = |c_b|^2/|a_bb|^2 - // So we iterate through the smaller wavefunction, to get the normalizing constant, - // then normalize both + // Internal type used to store operations with bitsets + // instead of vectors of qubit ids + using internal_operation = condensed_operation; + + // Hash table of the wavefunction + wavefunction _qubit_data; + + // Internal random numbers + std::function _rng; + + // Threshold to assert that something is zero when asserting it is 0 + double _precision = 1e-5; + // Threshold at which something is zero when + // deciding whether to add it into the superposition + double _rotation_precision = 1e-6; + // Often we compare to norms, so the precision must be squared + double _precision_squared = _precision * _precision; + double _rotation_precision_squared = _rotation_precision*_rotation_precision; + + // Normalizer for H and T gates (1/sqrt(2) as an amplitude) + const amplitude _normalizer = amplitude(1.0, 0.0) / std::sqrt(2.0); + + // The default for bytell_hash_map + // Used when allocating new wavefunctions + float _load_factor = 0.9375; + + // Makes a wavefunction that is preallocated to the right size + // and has the correct load factor + wavefunction make_wavefunction() { + wavefunction data((size_t)(_qubit_data.size() / _load_factor)); + data.max_load_factor(_load_factor); + return data; + } + wavefunction make_wavefunction(uint64_t n_states) { + wavefunction data((size_t)(n_states / _load_factor)); + data.max_load_factor(_load_factor); + return data; + } + + // Creates a qubit_label as a bit mask from a set of indices + qubit_label _get_mask(std::vector const& indices){ + return get_mask(indices); + } + + //********* Member variables for multithreading permutations ******// + // Shared iterator through the wavefunction + // Should point to _qubit_data.end() unless in use + decltype(_qubit_data.begin()) _current_state; + // Vector of waiting threads + std::vector _thread_pool; + // Condition variable, locked on _state_mtx + std::condition_variable cv; + // Mutex to prevent concurrent access to _current_state + std::mutex _state_mtx; + // Mutex to prevent concurrent access to _new_qubit_data + std::mutex _wfn_mtx; + // Used as a flag to decide when the thread pool is necessary, + // i.e., it is only set to false when the simulator is being destroyed + std::atomic in_use = true; + // Counts running threads, so controlling thread knows when the permutation is finished + std::atomic _running_threads = 0; + // A vector of operations that all threads point to to read which operations to apply + std::vector _operation_vector; + // New wavefunction which will replace the old one + wavefunction _new_qubit_data; + // The number of elements of the hash map that each thread will handle in one iteration + size_t _jump_size; + // Max number of threads + int max_num_threads = std::thread::hardware_concurrency(); + + // Called on creation, starts the thread pool and waits for all threads to + // finish initialization tasks + void _initialize_threads() { + // Unnecessary if OpenMP Is available + #ifndef DOMP_GE_V3 + // Prevents spurious wake-ups + _current_state = _qubit_data.end(); + _thread_pool = std::vector(); + #ifndef _OPENMP + int num_threads = std::thread::hardware_concurrency(); + #else + int num_threads = 1; + #pragma omp parallel + { + #pragma omp single + num_threads = omp_get_num_threads(); + } + #endif + // Lock the state, so that all the new threads get stopped during their execution + std::unique_lock state_lock(_state_mtx); + for (int i = 0; i < num_threads; i++){ + // Each running thread will decrement this just before it waits on the condition variable + ++_running_threads; + _thread_pool.push_back(std::thread([this](){this->_wait_to_permute();})); + } + // Waits until all the running threads have finished initializing + // This ensures the simulation does not actually start + // until the thread pool has finished initializing + cv.wait(state_lock, [&](){return _running_threads==0;}); + #endif + } + + // Function that all threads in _thread_pool run + // Waits on the condition variable, then applies _operation_vector + // to the qubit data + void _wait_to_permute(){ + // Prep local variables + auto local_state = _qubit_data.begin(); + auto local_end = _qubit_data.end(); + // Initialize a wfn lock to avoid recreating this object in every loop + std::unique_lock wfn_lock(_wfn_mtx); + wfn_lock.unlock(); + // First lock the state so we can wait on it + std::unique_lock state_lock(_state_mtx); + // Set the _current_state to ensure it waits + _current_state = _qubit_data.end(); + // Reduce number of running threads + _running_threads--; + // Notify all (i.e., the constructor thread) + // which will wait until the lock is released to check the number of running threads + cv.notify_all(); + // Wait until we set the current state to the start of the wavefunction + cv.wait(state_lock, [&](){ + return _current_state != _qubit_data.end(); + }); + + // Loop continuously + while (in_use){ + ++_running_threads; + // Loop through the wavefunction + while (_current_state != _qubit_data.end()){ + // Update local variables + local_state = _current_state; + _current_state.jump_forward(_jump_size); + local_end = _current_state; + // Unlock state to allow other threads to modify their state + state_lock.unlock(); + // Loop through the chunk that this thread has taken + for (; local_state != local_end; ++local_state){ + qubit_label label = local_state->first; + amplitude val = local_state->second; + + for (auto op : _operation_vector){ + switch (op.gate_type) { + case OP::X: + label.flip(op.target); + break; + case OP::MCX: + if ((op.controls & label) == op.controls){ + label.flip(op.target); + } + break; + case OP::Y: + label.flip(op.target); + val *= (label[op.target]) ? 1i : -1i; + break; + case OP::MCY: + if ((op.controls & label) == op.controls){ + label.flip(op.target); + val *= (label[op.target]) ? 1i : -1i; + } + break; + case OP::Z: + val *= (label[op.target] ? -1 : 1); + break; + case OP::MCZ: + val *= ((op.controls & label) == op.controls) ? -1 : 1; + break; + case OP::Phase: + val *= label[op.target] ? op.phase : 1; + break; + case OP::MCPhase: + val *= ((op.controls & label) == op.controls) ? op.phase : 1; + break; + case OP::SWAP: + if (label[op.target] != label[op.target_2]){ + label.flip(op.target); + label.flip(op.target_2); + } + break; + case OP::MCSWAP: + if (((label & op.controls) == op.controls) && (label[op.target] != label[op.target_2])){ + label.flip(op.target); + label.flip(op.target_2); + } + break; + case OP::Assert: + if (poppar(label & op.controls) != op.result && std::norm(val) > _precision_squared){ + std::cout << "Problematic state: " << label << "\n"; + std::cout << "Amplitude: " << val << "\n"; + std::cout << "Wavefunction size: " << _qubit_data.size() << "\n"; + throw std::runtime_error("Assert failed"); + } + break; + default: + throw std::runtime_error("Unsupported operation"); + break; + } + } + + wfn_lock.lock(); + _new_qubit_data.emplace(label, val); + wfn_lock.unlock(); + } + // Lock before checking and modifying the current state + state_lock.lock(); + } + // This thread has finished iterating through the wavefunction + --_running_threads; + // Notify all other threads (i.e., the controlling thread) + // if this is the last thread to finish + // The check on _running_threads avoids spurious wake-ups to + // the controlling thread, but takes a bit of time as _running_threads + // is atomic + if (_running_threads ==0) + cv.notify_all(); + // Wait on the state lock before repeating the loop + // The wait unlocks state_lock + cv.wait(state_lock, [&](){return _current_state != _qubit_data.end();}); + } + } + + + // Split the wavefunction if separable, otherwise return false + // Idea is that if we have a_bb|b1>|b2> as the first state, then for + // any other state a_xx|x1>|x2>, we must also have a_xb|x1>|b2> and a_bx|b1>|x2> + // in superposition. + // Also, the coefficients must separate as a_bb=c_b*d_b and a_xx = c_x*d_x, implying + // that a_xb = c_x*d_b and a_bx = c_b * d_x, and thus we can check this holds if + // a_bb*a_xx = a_bx * a_xb. + // If this holds: we write (a_xx/a_bx)|x1> into the first wavefunction and (a_xx/a_xb)|x2> + // into the second. + // This means the coefficients of wfn1 are all of the form (c_x/c_b); + // Thus, the norm of wfn1 will be 1/|c_b|^2; thus the norm of wfn2 is 1/|d_b|^2 = |c_b|^2/|a_bb|^2 + // So we iterate through the smaller wavefunction, to get the normalizing constant, + // then normalize both bool _split_wavefunction(qubit_label first_mask, wavefunction &wfn1, wavefunction &wfn2){ - qubit_label second_mask = ~first_mask; - // Guesses size - wfn1 = wavefunction((int)std::sqrt(_qubit_data.size())); - wfn2 = wavefunction((int)std::sqrt(_qubit_data.size())); - // base_label_1 = b1 and base_label_2 = b2 in the notation above - auto base_state = _qubit_data.begin(); - qubit_label base_label_1 = base_state->first & first_mask; - qubit_label base_label_2 = base_state->first & second_mask; - // base_val = a_bb - amplitude base_val = base_state->second; - double normalizer_1 = 0.0; - double normalizer_2 = 0.0; - // From here on, base_state is |x1>|x2> - ++base_state; - for (; base_state != _qubit_data.end(); ++base_state){ - qubit_label label_1 = base_state->first & first_mask; - qubit_label label_2 = base_state->first & second_mask; - // first_state is |x1>|b2>, second_state is |b1>|x2> - auto first_state = _qubit_data.find(label_1 | base_label_2); - auto second_state = _qubit_data.find(base_label_1 | label_2); - // Ensures that both |x1>|b2> and |b1>|x2> are in the superposition - if (first_state == _qubit_data.end() || second_state == _qubit_data.end()){ - // state does not exist - // therefore states are entangled - return false; - } else { // label with base label exists - // Checks that a_bba_xx = a_xb*a_bx - if (std::norm(first_state->second * second_state->second - base_val * base_state->second) > _precision_squared*_precision_squared){ - return false; - } else { - // Not entangled so far, save the two states, with amplitudes a_xx/a_bx and a_xx/a_xb, respectively - wfn1[label_1] = base_state->second / second_state->second; - wfn2[label_2] = base_state->second / first_state->second; - } - } - } - // Normalize - // This cannot be done in the previous loop, as that loop will encounter the same data several times - wavefunction &smaller_wfn = (wfn1.size() < wfn2.size()) ? wfn1 : wfn2; - wavefunction &larger_wfn = (wfn1.size() < wfn2.size()) ? wfn2 : wfn1; - for (auto current_state = smaller_wfn.begin(); current_state != smaller_wfn.end(); ++current_state){ - normalizer_1 += std::norm(current_state->second); - } - normalizer_2 = normalizer_1/std::norm(base_val); - normalizer_1 = 1.0/normalizer_1; - for (auto current_state = smaller_wfn.begin(); current_state != smaller_wfn.end(); ++current_state){ - current_state->second *= normalizer_1; - } - for (auto current_state = larger_wfn.begin(); current_state != larger_wfn.end(); ++current_state){ - current_state->second *= normalizer_2; - } - return true; + qubit_label second_mask = ~first_mask; + // Guesses size + wfn1 = wavefunction((int)std::sqrt(_qubit_data.size())); + wfn2 = wavefunction((int)std::sqrt(_qubit_data.size())); + // base_label_1 = b1 and base_label_2 = b2 in the notation above + auto base_state = _qubit_data.begin(); + qubit_label base_label_1 = base_state->first & first_mask; + qubit_label base_label_2 = base_state->first & second_mask; + // base_val = a_bb + amplitude base_val = base_state->second; + double normalizer_1 = 0.0; + double normalizer_2 = 0.0; + // From here on, base_state is |x1>|x2> + ++base_state; + for (; base_state != _qubit_data.end(); ++base_state){ + qubit_label label_1 = base_state->first & first_mask; + qubit_label label_2 = base_state->first & second_mask; + // first_state is |x1>|b2>, second_state is |b1>|x2> + auto first_state = _qubit_data.find(label_1 | base_label_2); + auto second_state = _qubit_data.find(base_label_1 | label_2); + // Ensures that both |x1>|b2> and |b1>|x2> are in the superposition + if (first_state == _qubit_data.end() || second_state == _qubit_data.end()){ + // state does not exist + // therefore states are entangled + return false; + } else { // label with base label exists + // Checks that a_bba_xx = a_xb*a_bx + if (std::norm(first_state->second * second_state->second - base_val * base_state->second) > _precision_squared*_precision_squared){ + return false; + } else { + // Not entangled so far, save the two states, with amplitudes a_xx/a_bx and a_xx/a_xb, respectively + wfn1[label_1] = base_state->second / second_state->second; + wfn2[label_2] = base_state->second / first_state->second; + } + } + } + // Normalize + // This cannot be done in the previous loop, as that loop will encounter the same data several times + wavefunction &smaller_wfn = (wfn1.size() < wfn2.size()) ? wfn1 : wfn2; + wavefunction &larger_wfn = (wfn1.size() < wfn2.size()) ? wfn2 : wfn1; + for (auto current_state = smaller_wfn.begin(); current_state != smaller_wfn.end(); ++current_state){ + normalizer_1 += std::norm(current_state->second); + } + normalizer_2 = normalizer_1/std::norm(base_val); + normalizer_1 = 1.0/normalizer_1; + for (auto current_state = smaller_wfn.begin(); current_state != smaller_wfn.end(); ++current_state){ + current_state->second *= normalizer_1; + } + for (auto current_state = larger_wfn.begin(); current_state != larger_wfn.end(); ++current_state){ + current_state->second *= normalizer_2; + } + return true; } - // Iterates through a wavefunction and calls the output function on each value - // It first sorts the labels before outputting - void _DumpWavefunction_base(wavefunction &wfn, std::function output){ - if (wfn.size() == 0){ return; } - std::vector sortedLabels; - sortedLabels.reserve(wfn.size()); - for (auto current_state = (wfn).begin(); current_state != (wfn).end(); ++current_state) { - sortedLabels.push_back(current_state->first); - } - std::sort( - sortedLabels.begin(), - sortedLabels.end(), - [](const qubit_label& lhs, const qubit_label& rhs){return lhs < rhs;}); - amplitude val; - for (qubit_label label : sortedLabels){ - output(label, _qubit_data.find(label)->second); - - } - } + // Iterates through a wavefunction and calls the output function on each value + // It first sorts the labels before outputting + void _DumpWavefunction_base(wavefunction &wfn, std::function output){ + if (wfn.size() == 0){ return; } + std::vector sortedLabels; + sortedLabels.reserve(wfn.size()); + for (auto current_state = (wfn).begin(); current_state != (wfn).end(); ++current_state) { + sortedLabels.push_back(current_state->first); + } + std::sort( + sortedLabels.begin(), + sortedLabels.end(), + [](const qubit_label& lhs, const qubit_label& rhs){return lhs < rhs;}); + amplitude val; + for (qubit_label label : sortedLabels){ + output(label, _qubit_data.find(label)->second); + + } + } }; -} // namespace SPARSESIMULATOR -} // namespace Quantum -} // namespace Microsoft +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/types.h b/src/Simulation/Simulators/SparseSimulator/Native/types.h index 17116a51bc5..888ebf2abd3 100644 --- a/src/Simulation/Simulators/SparseSimulator/Native/types.h +++ b/src/Simulation/Simulators/SparseSimulator/Native/types.h @@ -1,51 +1,45 @@ -// Copyright (c) Microsoft Corporation. All rights reserved. +// Copyright (c) Microsoft Corporation. // Licensed under the MIT License. -#pragma once +#pragma once #include #include #include #include -#include "quantum_hash_map.hpp" #include #include -namespace Microsoft -{ - namespace Quantum - { - namespace SPARSESIMULATOR - { - using mutex_type = std::mutex; - using recursive_mutex_type = std::recursive_mutex; - using lock_type = std::lock_guard; - using recursive_lock_type = std::lock_guard; +#include "quantum_hash_map.hpp" - #ifndef USE_SINGLE_PRECISION - using RealType = double; - #else - using RealType = float; - #endif +namespace Microsoft::Quantum::SPARSESIMULATOR +{ +using mutex_type = std::mutex; +using recursive_mutex_type = std::recursive_mutex; +using lock_type = std::lock_guard; +using recursive_lock_type = std::lock_guard; +#ifndef USE_SINGLE_PRECISION + using RealType = double; +#else + using RealType = float; +#endif - // Logical qubit id is visible to the clients and is immutable during the lifetime of the qubit. - using logical_qubit_id = unsigned; +// Logical qubit id is visible to the clients and is immutable during the lifetime of the qubit. +using logical_qubit_id = unsigned; - using amplitude = std::complex; +using amplitude = std::complex; - template - using qubit_label_type = std::bitset; +template +using qubit_label_type = std::bitset; - // Wavefunctions are hash maps of some key (std::bitset or a string) - template - using abstract_wavefunction = ska::bytell_hash_map>; +// Wavefunctions are hash maps of some key (std::bitset or a string) +template +using abstract_wavefunction = ska::bytell_hash_map>; - // Wavefunctions with strings as keys are "universal" in that they do not depend - // on the total number of qubits - using universal_wavefunction = abstract_wavefunction; +// Wavefunctions with strings as keys are "universal" in that they do not depend +// on the total number of qubits +using universal_wavefunction = abstract_wavefunction; - } // namespace SPARSESIMULATOR - } // namespace Quantum -} // namespace Microsoft +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/README.md b/src/Simulation/Simulators/SparseSimulator/README.md index 2705828bee0..4b5a42eb3f0 100644 --- a/src/Simulation/Simulators/SparseSimulator/README.md +++ b/src/Simulation/Simulators/SparseSimulator/README.md @@ -3,11 +3,11 @@ This is a an alternative quantum simulator, compatible with Q\# and Microsoft's # Repository Structure -/SparseQuantumSimulator/ -|--/Native/: C++ code for the simulator -|--/SparseSimulatorCS/: C# library for the simulator -|--/SparseSimQSharpTests/: Several Q# tests of the simulator -|--/SparseSimulatorTests/: C++ tests of the simulator +/SparseQuantumSimulator/ +|--/Native/: C++ code for the simulator +|--/SparseSimulatorCS/: C# library for the simulator +|--/SparseSimQSharpTests/: Several Q# tests of the simulator +|--/SparseSimulatorTests/: C++ tests of the simulator # Dependencies - CMake v1.2.0 @@ -75,7 +75,7 @@ To add a gate to the simulator, you will need to add it to: Optionally, you may need to add a Q\# and C\# file to the SparseSimulatorCS folder, to create a Q\# operation to call the gate if it does not already exist. -## Quantum State (quantum_state.hpp and basec_quantum_state.hpp) +## Quantum State (quantum_state.hpp and basic_quantum_state.hpp) This code will modify the actual wavefunction data structure. The code for the `H` gate will be a good template. Typically, such a function will create a new wavefunction object, then iterate through the existing wavefunction to compute new (label, amplitude) pairs, which are inserted into the new wavefunction. Finally, it moves the new data into the old wavefunction object. Since the hash map does not allow modifications while iterating through it, this is the only way to implement most gates. @@ -104,7 +104,8 @@ Here `Id` is an internal variable for the simulator's ID, which must be passed t If you are implementing a gate that Q\# already expects (i.e., it is initialized in the `QuantumProcessorBase` class) then the previous steps will be enough. However, if you want to create an entirely new gate, you will need to create a Q\# operation to call it. "Probes.qs" and "Probes.cs" provide a template for how this code will look. In Q\#, declare the operation `YourNewGate` with code `body intrinsic;`. Then in a separate C\# file, use the following template: -`public partial class YourNewGate +```C# +public partial class YourNewGate { public class Native : YourNewGate { @@ -113,21 +114,24 @@ If you are implementing a gate that Q\# already expects (i.e., it is initialized { sim = m as SparseSimulator; } - public override Func<\#input types from Q\#\#, \#output types to Q\#\#> __Body__ => sim == null ? base.__Body__ : (args) => { + public override Func<'input types from Q#', 'output types to Q#'> __Body__ => sim == null ? base.__Body__ : (args) => { return sim.Your_new_gate(args.Item1, args.Item2,...); }; } -}` +} +``` You will also need to add -`public partial class SparseSimulator : QuantumProcessorDispatcher, IDisposable +```c# +public partial class SparseSimulator : QuantumProcessorDispatcher, IDisposable { - public \# Q\# return type\# Your_new_gate(\#args\#) + public 'Q# return type' Your_new_gate('args') { - return ((SparseSimulatorProcessor)this.QuantumProcessor).Your_new_gate(\#args\#); + return ((SparseSimulatorProcessor)this.QuantumProcessor).Your_new_gate('args'); } } -` +``` + which tells the `SparseSimulator` class to forward the call to its internal `SparseSimulatorProcessor` class. diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs b/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs index 0893c4440b3..24712e556bb 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs @@ -1,4 +1,7 @@ -namespace Microsoft.Quantum.SparseSimulatorTests { +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.SparseSimulatorTests { open Microsoft.Quantum.Canon; open Microsoft.Quantum.Intrinsic; diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs index d058eb86c5d..a9248bdf84f 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs @@ -1,4 +1,7 @@ -// Functions that access the internal wavefunction for diagnostic purposes +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +// Functions that access the internal wavefunction for diagnostic purposes using Microsoft.Quantum.Simulation; using Microsoft.Quantum.Simulation.Common; diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs index 8c0f5ee3c6e..94304b5407f 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs @@ -1,4 +1,7 @@ -namespace Microsoft.Quantum.SparseSimulation { +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.SparseSimulation { open Microsoft.Quantum.Canon; open Microsoft.Quantum.Intrinsic; diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs index 466d8d6b2a0..8512ef6853b 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs @@ -1,4 +1,7 @@ -// Overrides the basic AND gate for faster execution +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +// Overrides the basic AND gate for faster execution using Microsoft.Quantum.Simulation.Core; using System; diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs index 04ad05ca62f..2478e2d8e26 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs @@ -1,4 +1,7 @@ -using System; +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +using System; using System.Collections.Generic; using Microsoft.Quantum.Simulation.QuantumProcessor; using Microsoft.Quantum.Simulation.Common; diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp index f94750c84e5..409e747a77e 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp @@ -1,10 +1,11 @@ -#pragma once +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. -#include "../SparseQuantumSimulator/SparseSimulator.h" -#include "../SparseQuantumSimulator/capi.hpp" -#include "../SparseQuantumSimulator/capi.cpp" // yes really -#include "../SparseQuantumSimulator/factory.hpp" -#include "../SparseQuantumSimulator/factory.cpp" +#include "../Native/SparseSimulator.h" +#include "../Native/capi.hpp" +#include "../Native/capi.cpp" // yes really +#include "../Native/factory.hpp" +#include "../Native/factory.cpp" #include "TestHelpers.hpp" #include "CppUnitTest.h" diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp index c7b61dd96ad..c43956686e5 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp @@ -1,8 +1,9 @@ -#pragma once +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. #include "pch.h" #include "CppUnitTest.h" -#include "../SparseQuantumSimulator/SparseSimulator.h" +#include "../Native/SparseSimulator.h" #include "TestHelpers.hpp" #include #include diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj index 37fe5ddeac1..18ea6d5e2ad 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj @@ -96,6 +96,7 @@ WIN32;_DEBUG;%(PreprocessorDefinitions) true pch.h + stdcpp17 Windows diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp index 539612fd9fa..e2d8135a0df 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp @@ -1,8 +1,9 @@ -#pragma once +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. #include "pch.h" #include "CppUnitTest.h" -#include "../SparseQuantumSimulator/SparseSimulator.h" +#include "../Native/SparseSimulator.h" #include #include diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp index 433a97413b4..8eb71c0f7b1 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp @@ -1,8 +1,11 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + #pragma once #include "pch.h" #include "CppUnitTest.h" -#include "../SparseQuantumSimulator/SparseSimulator.h" +#include "../Native/SparseSimulator.h" #include #include diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp index 64b7eef6d6b..2d641a7e3b9 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp @@ -1,3 +1,6 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + // pch.cpp: source file corresponding to the pre-compiled header #include "pch.h" diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h index 61302cbcf59..2cd3c85c965 100644 --- a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h @@ -1,3 +1,6 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + // pch.h: This is a precompiled header file. // Files listed below are compiled only once, improving build performance for future builds. // This also affects IntelliSense performance, including code completion and many code browsing features.