diff --git a/src/Simulation/Simulators/Microsoft.Quantum.Simulators.csproj b/src/Simulation/Simulators/Microsoft.Quantum.Simulators.csproj index 73e638180e8..0869cee10f5 100644 --- a/src/Simulation/Simulators/Microsoft.Quantum.Simulators.csproj +++ b/src/Simulation/Simulators/Microsoft.Quantum.Simulators.csproj @@ -46,4 +46,11 @@ + + + + + + + diff --git a/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt b/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt new file mode 100644 index 00000000000..12ce6e42ce3 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/CMakeLists.txt @@ -0,0 +1,37 @@ +cmake_minimum_required(VERSION 3.10) +project(SparseQuantumSimulator) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +set(CMAKE_MACOSX_RPATH 1) +# Main build files +add_library(SparseQuantumSimulator SHARED factory.cpp capi.cpp) + +# Set OpenMP if it is available +find_package(OpenMP REQUIRED) +if(OpenMP_CXX_FOUND) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") + if (OpenMP_CXX_VERSION_MAJOR GREATER_EQUAL 3) + target_compile_definitions(SparseQuantumSimulator PRIVATE DOMP_GE_V3=1) + endif() +endif() + +# Windows adds a special dllexport command which must be defined +if (WIN32) + target_compile_options(SparseQuantumSimulator PUBLIC -fdeclspec) + target_compile_definitions(SparseQuantumSimulator PRIVATE BUILD_DLL=1) +endif() +# Try to optimize with gcc +if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") + target_compile_options(SparseQuantumSimulator PUBLIC -O3 -ftree-vectorize -mavx2 -mfma -fopenmp) +endif() + +# Sets NDEBUG to true for all non-debug settings +if (NOT (CMAKE_BUILD_TYPE:STRING STREQUAL "Debug")) + target_compile_definitions(SparseQuantumSimulator PRIVATE NDEBUG=1) +endif() + +message("Compiler flags: ${CMAKE_CXX_FLAGS_RELEASE}") diff --git a/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h b/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h new file mode 100644 index 00000000000..25e524e9c70 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/SparseSimulator.h @@ -0,0 +1,851 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "quantum_state.hpp" +#include "basic_quantum_state.hpp" +#include "types.h" +#include "gates.h" + +using namespace std::literals::complex_literals; + +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())); +} + +// 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); +} + +// Sparse simulator only stores non-zero coefficients of the quantum state. +// It has good performance only when the number of non-zero coefficients is low. +// If the number of non-zero coefficients is low, the number of qubits may be fairly large. +// Sparse simulator employs std::unordered_map hashtable. +// Keys are basis vectors represented by std::bitset<>. +// Values are non-zero amplitudes represented by std::complex. +// Zero amplitudes are simply not stored. +// Hashtable is reallocated and reconstructed on almost every gate. +// Reallocation is saved for some gates that can be performed in one round. +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; + + _queue_Ry = std::vector(num_qubits, 0); + _queue_Rx = std::vector(num_qubits, 0); + _queue_H = std::vector(num_qubits, 0); + _angles_Rx = std::vector(num_qubits, 0.0); + _angles_Ry = 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); + + num_qubits = _quantum_state->get_num_qubits(); + _occupied_qubits.resize(num_qubits, 0); + _queue_Ry.resize(num_qubits, 0); + _queue_Rx.resize(num_qubits, 0); + _queue_H.resize(num_qubits, 0); + _angles_Rx.resize(num_qubits, 0.0); + _angles_Ry.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"); + } + + } + _set_qubit_to_zero(qubit_id); + return true; + } + + + void X(logical_qubit_id index) { + // XY = - YX + if (_queue_Ry[index]){ + _angles_Ry[index] *= -1.0; + } + // Rx trivially commutes + if (_queue_H[index]) { + _queued_operations.push_back(operation(OP::Z, index)); + return; + } + _queued_operations.push_back(operation(OP::X, index)); + _set_qubit_to_nonzero(index); + } + + + // 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 (_queue_Ry[controls[0]] || _queue_Rx[controls[0]] || (_queue_H[controls[0]] && !_queue_H[target])){ + _execute_queued_ops(controls, OP::Ry); + } + } + // Ry on the target causes issues + if (_queue_Ry[target]){ + _execute_queued_ops(target, OP::Ry); + } + // Rx on the target trivially commutes + + // An H on the target flips the operation + if (_queue_H[target]){ + // If it is a CNOT and there is also an H on the control, we swap control and target + if (controls.size() == 1 && _queue_H[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); + } + + // Same as MCX, but we assert that the target is 0 before execution + void MCApplyAnd(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 MCApplyAndAdj(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 (_queue_Rx[index]){ + _angles_Rx[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 (_queue_Rx[target]){ + _execute_queued_ops(target, OP::Rx); + } + // HY = -YH, so we add a phase to track this + if (_queue_H[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); + } + + + void Z(logical_qubit_id index) { + // ZY = -YZ + if (_queue_Ry[index]){ + _angles_Ry[index] *= -1; + } + // XZ = -ZX + if (_queue_Rx[index]){ + _angles_Rx[index] *= -1; + } + // HZ = XH + if (_queue_H[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 (_queue_Ry[control] || _queue_Rx[control]){ + count += 2; + } + if (_queue_H[control]){ + count++; + } + } + if (_queue_Ry[target] || _queue_Rx[target]){ + count +=2; + } + if (_queue_H[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 (_queue_H[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 (_queue_Ry[index] || _queue_Ry[index] || _queue_H[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)); + } + + 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); + } + + + void R1(double const& angle, logical_qubit_id index) { + Phase(amplitude(std::cos(angle), std::sin(angle)), 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 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 AdjS(logical_qubit_id index) { + Phase(-1i, index); + } + + + + void R(Gates::Basis b, double phi, logical_qubit_id index) + { + if (b == Gates::Basis::PauliI){ + return; + } + + // 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){ + _queue_Ry[index] = true; + _angles_Ry[index] += phi; + _set_qubit_to_nonzero(index); + return; + } else if (_queue_Ry[index]) { + _execute_queued_ops(index, OP::Ry); + } + + if (b == Gates::Basis::PauliX){ + _queue_Rx[index] = true; + _angles_Rx[index] += phi; + _set_qubit_to_nonzero(index); + return; + } else if (_queue_Rx[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 (_queue_H[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 (_queue_Ry[target] && b != Gates::Basis::PauliY){ + _execute_queued_ops(target, OP::Ry); + } + if (_queue_Rx[target] && b != Gates::Basis::PauliX){ + _execute_queued_ops(target, OP::Rx); + } + if (_queue_H[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 * std::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 H(logical_qubit_id index) { + // YH = -HY + _angles_Ry[index] *= (_queue_Ry[index] ? -1.0 : 1.0); + // Commuting with Rx creates a phase, but on the wrong side + // So we execute any Rx immediately + if (_queue_Rx[index]){ + _execute_queued_ops(index, OP::Rx); + } + _queue_H[index] = !_queue_H[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 (_queue_Ry[target] || _queue_Rx[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); + } + + + + + 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 + _queue_Ry.swap(_queue_Ry[index_1], _queue_Ry[index_2]); + std::swap(_angles_Ry[index_1], _angles_Ry[index_2]); + _queue_Rx.swap(_queue_Rx[index_1], _queue_Rx[index_2]); + std::swap(_angles_Rx[index_1], _angles_Rx[index_2]); + _queue_H.swap(_queue_H[index_1], _queue_H[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]; + + } + + 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 (_queue_Rx[qubit] || _queue_Ry[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 (_queue_H[qubits[i]]){ + result ^= _queue_H[qubits[i]]; + } + isAllZ = false; + isEmpty = false; + break; + case Gates::Basis::PauliX: + // HX = ZH + if (_queue_H[qubits[i]]){ + axes[i] = Gates::Basis::PauliZ; + } else { + isAllZ = false; + } + isEmpty = false; + break; + case Gates::Basis::PauliZ: + // HZ = XH + if (_queue_H[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 defined(NDEBUG) + if (isAllZ) { + _queued_operations.push_back(operation(OP::Assert, qubits, result)); + // In release mode we queue Z assertion and return. + return; + } + #endif + // X or Y assertions require immediate execution. + // We also execute Z assertion immediately in debug mode + // to provide feedback at the point of failure. + _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]); + } + } + return result; + } + + // Returns the amplitude of a given bitstring + amplitude probe(std::string const& 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 const& 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 _queue_H; + std::vector _queue_Rx; + std::vector _queue_Ry; + + std::vector _angles_Rx; + std::vector _angles_Ry; + + // 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 (_queue_H[index]){ + _quantum_state->H(index); + _queue_H[index] = false; + } + if (_queue_Rx[index]){ + _quantum_state->R(Gates::Basis::PauliX, _angles_Rx[index], index); + _angles_Rx[index] = 0.0; + _queue_Rx[index] = false; + } + if (_queue_Ry[index]){ + _quantum_state->R(Gates::Basis::PauliY, _angles_Ry[index], index); + _angles_Ry[index] = 0.0; + _queue_Ry[index] = false; + } + } + + void _execute_RxH_single_qubit(logical_qubit_id const &index){ + if (_queue_H[index]){ + _quantum_state->H(index); + _queue_H[index] = false; + } + if (_queue_Rx[index]){ + _quantum_state->R(Gates::Basis::PauliX, _angles_Rx[index], index); + _angles_Rx[index] = 0.0; + _queue_Rx[index] = false; + } + } + + void _execute_H_single_qubit(logical_qubit_id const &index){ + if (_queue_H[index]){ + _quantum_state->H(index); + _queue_H[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); + 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 const& 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; + } + } + + + // 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 (_queue_Ry[target] || _queue_Rx[target] || _queue_H[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 (_queue_Ry[control] || _queue_Rx[control] || _queue_H[control]){ + _execute_queued_ops(controls, OP::Ry); + return; + } + } + } + +}; + +} // 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 new file mode 100644 index 00000000000..f1e68bf3227 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/basic_quantum_state.hpp @@ -0,0 +1,74 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +#include + +#include "types.h" +#include "gates.h" + +namespace Microsoft::Quantum::SPARSESIMULATOR +{ + +// Virtual class for QuantumState +// This is not templated, so it allows SparseSimulator types to avoid templates +class BasicQuantumState +{ +public: + + BasicQuantumState() {} + + virtual logical_qubit_id get_num_qubits() = 0; + + virtual void DumpWavefunction(size_t indent = 0) = 0; + + virtual void set_random_seed(unsigned seed = std::mt19937::default_seed) = 0; + + virtual void set_precision(double new_precision) = 0; + + virtual float get_load_factor() = 0; + + virtual void set_load_factor(float new_load_factor) = 0; + + virtual size_t get_wavefunction_size() = 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 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 amplitude probe(std::string const& label) = 0; + + virtual bool dump_qubits(std::vector const& qubits, void (*callback)(char*, double, double)) = 0; + + virtual void dump_all(logical_qubit_id max_qubit_id, void (*callback)(char*, double, double)) = 0; + + virtual void phase_and_permute(std::listconst &) = 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 H(logical_qubit_id index) = 0; + virtual void MCH(std::vector const& controls, logical_qubit_id index) = 0; + + virtual bool is_qubit_zero(logical_qubit_id) = 0; + + virtual universal_wavefunction get_universal_wavefunction() = 0; + + virtual std::function get_rng() = 0; + + virtual std::string Sample() = 0; +}; + +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp b/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp new file mode 100644 index 00000000000..92c153d4443 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/capi.cpp @@ -0,0 +1,306 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +// Wrapper functions for basic C++ functions +// All have the same logic: use the sim_id argument as an +// index into the vector of simulators, +// then call a member function + +#include +#include +#include +#include + +#include "capi.hpp" +#include "SparseSimulator.h" +#include "factory.hpp" + +using namespace Microsoft::Quantum::SPARSESIMULATOR; + +std::string sample_string; + +extern "C" +{ + + MICROSOFT_QUANTUM_DECL unsigned init_cpp(logical_qubit_id num_qubits) + { + return createSimulator(num_qubits); + } + + + MICROSOFT_QUANTUM_DECL void destroy_cpp(unsigned sim_id) + { + destroySimulator(sim_id); + } + + MICROSOFT_QUANTUM_DECL void seed_cpp(unsigned sim_id, _In_ unsigned s){ + getSimulator(sim_id)->set_random_seed(s); + } + + MICROSOFT_QUANTUM_DECL void allocateQubit_cpp(unsigned sim_id, logical_qubit_id q) + { + getSimulator(sim_id)->allocate_specific_qubit(q); + } + + MICROSOFT_QUANTUM_DECL void releaseQubit_cpp(unsigned sim_id, logical_qubit_id q) + { + getSimulator(sim_id)->release(q); + } + + MICROSOFT_QUANTUM_DECL logical_qubit_id num_qubits_cpp(unsigned sim_id) + { + return getSimulator(sim_id)->get_num_qubits(); + } + +// Generic single-qubit gate +#define FWDGATE1(G) \ + MICROSOFT_QUANTUM_DECL void G##_cpp(unsigned sim_id, _In_ logical_qubit_id q) \ + { \ + getSimulator(sim_id)->G(q); \ + } +// Generic multi-qubit gate +#define FWDCSGATE1(G) \ + MICROSOFT_QUANTUM_DECL void MC##G##_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q) \ + { \ + \ + getSimulator(sim_id)->MC##G(std::vector(c, c + n), q); \ + } +#define FWD(G) FWDGATE1(G) + + // single-qubit gates + FWD(X) + FWD(Y) + FWD(Z) + FWD(H) + + FWD(S) + FWD(T) + FWD(AdjS) + FWD(AdjT) + +#define MFWD(G) FWDCSGATE1(G) + MFWD(H) + MFWD(X) + MFWD(Y) + MFWD(Z) + +#undef FWDGATE1 +#undef FWDGATE2 +#undef FWDGATE3 +#undef FWDCSGATE1 +#undef FWD + + + + MICROSOFT_QUANTUM_DECL void SWAP_cpp(unsigned sim_id, _In_ logical_qubit_id q1, _In_ logical_qubit_id q2){ + getSimulator(sim_id)->SWAP(q1, q2); + } + + MICROSOFT_QUANTUM_DECL void MCSWAP_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q1, _In_ logical_qubit_id q2){ + getSimulator(sim_id)->CSWAP(std::vector(c, c + n), q1, q2); + } + + MICROSOFT_QUANTUM_DECL void MCApplyAnd_cpp(unsigned sim_id,_In_ int length, _In_reads_(length) logical_qubit_id* controls, _In_ logical_qubit_id target){ + getSimulator(sim_id)->MCApplyAnd(std::vector(controls, controls + length), target); + } + + MICROSOFT_QUANTUM_DECL void MCAdjointApplyAnd_cpp(unsigned sim_id,_In_ int length, _In_reads_(length) logical_qubit_id* controls, _In_ logical_qubit_id target){ + getSimulator(sim_id)->MCApplyAndAdj(std::vector(controls, controls + length), target); + } + + // rotations + + MICROSOFT_QUANTUM_DECL void R_cpp(unsigned sim_id, _In_ int b, _In_ double phi, _In_ logical_qubit_id q) + { + getSimulator(sim_id)->R(static_cast(b), phi, q); + } + MICROSOFT_QUANTUM_DECL void Rfrac_cpp(unsigned sim_id, _In_ int b, _In_ std::int64_t numerator, _In_ std::int64_t power, _In_ logical_qubit_id q) + { + getSimulator(sim_id)->RFrac(static_cast(b), numerator, power, q); + } + MICROSOFT_QUANTUM_DECL void R1_cpp(unsigned sim_id,_In_ double phi, _In_ logical_qubit_id q) + { + getSimulator(sim_id)->R1(phi, q); + } + MICROSOFT_QUANTUM_DECL void R1frac_cpp(unsigned sim_id, _In_ std::int64_t numerator, _In_ std::int64_t power, _In_ logical_qubit_id q) + { + getSimulator(sim_id)->R1Frac(numerator, power, q); + } + + // multi-controlled rotations + MICROSOFT_QUANTUM_DECL void MCR_cpp( + unsigned sim_id, + _In_ int b, + _In_ double phi, + _In_ logical_qubit_id nc, + _In_reads_(nc) logical_qubit_id* c, + _In_ logical_qubit_id q) + { + std::vector cv(c, c + nc); + getSimulator(sim_id)->MCR(cv, static_cast(b), phi, q); + } + + MICROSOFT_QUANTUM_DECL void MCRFrac_cpp( + unsigned sim_id, + _In_ int b, + _In_ std::int64_t numerator, + _In_ std::int64_t power, + _In_ logical_qubit_id nc, + _In_reads_(nc) logical_qubit_id* c, + _In_ logical_qubit_id q) + { + std::vector cv(c, c + nc); + getSimulator(sim_id)->MCRFrac(cv, static_cast(b), numerator, power, q); + } + + MICROSOFT_QUANTUM_DECL void MCR1_cpp( + unsigned sim_id, + _In_ double phi, + _In_ int nc, + _In_reads_(nc) logical_qubit_id* c, + _In_ logical_qubit_id q) + { + std::vector cv(c, c + nc); + getSimulator(sim_id)->MCR1(cv, phi, q); + } + + MICROSOFT_QUANTUM_DECL void MCR1Frac_cpp( + unsigned sim_id, + _In_ std::int64_t numerator, + _In_ std::int64_t power, + _In_ int nc, + _In_reads_(nc) logical_qubit_id* c, + _In_ logical_qubit_id q) + { + std::vector cv(c, c + nc); + getSimulator(sim_id)->MCR1Frac(cv, numerator, power, q); + } + + // Exponential of Pauli operators + MICROSOFT_QUANTUM_DECL void Exp_cpp( + unsigned sim_id, + _In_ int n, + _In_reads_(n) int* b, + _In_ double phi, + _In_reads_(n) logical_qubit_id* q) + { + std::vector bv; + bv.reserve(n); + for (int i = 0; i < n; ++i) + bv.push_back(static_cast(*(b + i))); + std::vector qv(q, q + n); + getSimulator(sim_id)->Exp(bv, phi, qv); + } + + MICROSOFT_QUANTUM_DECL void MCExp_cpp( + unsigned sim_id, + _In_ int nc, + _In_ int n, + _In_reads_(nc) logical_qubit_id* c, + _In_reads_(n) int* b, + _In_ double phi, + _In_reads_(n) logical_qubit_id* q) + { + std::vector bv; + bv.reserve(n); + for (int i = 0; i < n; ++i) + bv.push_back(static_cast(*(b + i))); + std::vector qv(q, q + n); + std::vector cv(c, c + nc); + getSimulator(sim_id)->MCExp(cv, bv, phi, qv); + } + + // measurements + MICROSOFT_QUANTUM_DECL bool M_cpp(unsigned sim_id, _In_ logical_qubit_id q) + { + return getSimulator(sim_id)->M(q); + } + + MICROSOFT_QUANTUM_DECL void Reset_cpp(unsigned sim_id, _In_ logical_qubit_id q){ + getSimulator(sim_id)->Reset(q); + } + + MICROSOFT_QUANTUM_DECL bool Measure_cpp( + unsigned sim_id, + _In_ int n, + _In_reads_(n) int* b, + _In_reads_(n) logical_qubit_id* q) + { + std::vector bv; + bv.reserve(n); + for (int i = 0; i < n; ++i) + bv.push_back(static_cast(*(b + i))); + std::vector qv(q, q + n); + return getSimulator(sim_id)->Measure(bv, qv); + } + + // Extracts the probability of measuring a One result on qubits q with basis b + MICROSOFT_QUANTUM_DECL double JointEnsembleProbability_cpp( + unsigned sim_id, + _In_ int n, + _In_reads_(n) int* b, + _In_reads_(n) logical_qubit_id* q) + { + std::vector bv; + bv.reserve(n); + for (int i = 0; i < n; ++i) + bv.push_back(static_cast(*(b + i))); + std::vector qv(q, q + n); + return getSimulator(sim_id)->MeasurementProbability(bv, qv); + } + + + // Iterates through the entire wavefunction and calls `callback` on every state in the superposition + // It will write the label of the state, in binary, from qubit 0 to `max_qubit_id`, into the char* pointer, then call `callback` + // with the real and complex values as the double arguments + MICROSOFT_QUANTUM_DECL void Dump_cpp(unsigned sim_id, _In_ logical_qubit_id max_qubit_id, _In_ void (*callback)(char* , double, double)){ + return getSimulator(sim_id)->dump_all(max_qubit_id, callback); + } + + // Same as Dump_cpp, but only dumps the wavefunction on the qubits in `q`, ensuring they are separable from the rest of the state first + MICROSOFT_QUANTUM_DECL bool DumpQubits_cpp( + unsigned sim_id, + _In_ logical_qubit_id n, + _In_reads_(n) logical_qubit_id* q, + _In_ void (*callback)(char* , double, double)) + { + std::vector qs(q, q + n); + return getSimulator(sim_id)->dump_qubits(qs, callback); + } + + // Writes the real and imaginary parts of the amplitude of the state `label` + // (a binary bistring of length `length`) to the arguments of `real` and `imag` + MICROSOFT_QUANTUM_DECL void GetAmplitude_cpp(_In_ unsigned sim_id, _In_ size_t length, _In_ const char* label, _Out_ double &real, _Out_ double &imag) + { + amplitude val = getSimulator(sim_id)->probe(std::string(label, length)); + real = val.real(); + imag = val.imag(); + } + + + // Asserts that the gates in `b`, measured on the qubits in `q`, return `result` + MICROSOFT_QUANTUM_DECL bool Assert_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) int* b, _In_reads_(n) logical_qubit_id* q, bool result){ + std::vector bv; + bv.reserve(n); + for (int i = 0; i < n; ++i) + bv.push_back(static_cast(*(b + i))); + std::vector qv(q, q + n); + try { + getSimulator(sim_id)->Assert(bv, qv, result); + } + catch(const std::exception&){ + // C# will not call "Dispose" + // after this exception, so this cleans up manually + destroySimulator(sim_id); + return false; + } + return true; + } + + // Returns a character string of the label of a randomly sampled state + MICROSOFT_QUANTUM_DECL const char* Sample_cpp(unsigned sim_id){ + sample_string = getSimulator(sim_id)->Sample(); + return sample_string.c_str(); + } + +} diff --git a/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp b/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp new file mode 100644 index 00000000000..8f4c3839041 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/capi.hpp @@ -0,0 +1,145 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +#include + +#include "types.h" + +// SAL only defined in windows. +#ifndef _In_ +#define _In_ +#define _In_reads_(n) +#endif +#ifndef _Out_ +#define _Out_ +#endif + +#ifdef BUILD_DLL +#define MICROSOFT_QUANTUM_DECL __declspec(dllexport) +#else +#define MICROSOFT_QUANTUM_DECL +#endif +#define MICROSOFT_QUANTUM_DECL_IMPORT __declspec(dllimport) + +using namespace Microsoft::Quantum::SPARSESIMULATOR; + +// All of these are called by the C# SparseSimulator class +extern "C" +{ + MICROSOFT_QUANTUM_DECL unsigned init_cpp(logical_qubit_id num_qubits); + MICROSOFT_QUANTUM_DECL void destroy_cpp(unsigned sim_id); + + MICROSOFT_QUANTUM_DECL void seed_cpp(unsigned sim_id, _In_ unsigned s); + // allocate and release + MICROSOFT_QUANTUM_DECL void allocateQubit_cpp(unsigned sim_id, logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void releaseQubit_cpp(unsigned sim_id, logical_qubit_id q); + MICROSOFT_QUANTUM_DECL logical_qubit_id num_qubits_cpp(unsigned sim_id); + + // single-qubit gates + MICROSOFT_QUANTUM_DECL void X_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void Y_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void Z_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void H_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void S_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void T_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void AdjS_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void AdjT_cpp(unsigned sim_id, _In_ logical_qubit_id q); + + + MICROSOFT_QUANTUM_DECL void MCX_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void MCY_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void MCZ_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void MCH_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q); + + MICROSOFT_QUANTUM_DECL void SWAP_cpp(unsigned sim_id, _In_ logical_qubit_id q1, _In_ logical_qubit_id q2); + MICROSOFT_QUANTUM_DECL void MCSWAP_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) logical_qubit_id* c, _In_ logical_qubit_id q1, _In_ logical_qubit_id q2); + MICROSOFT_QUANTUM_DECL void MCAnd_cpp(unsigned sim_id,_In_ int length, _In_reads_(length) logical_qubit_id* controls, _In_ logical_qubit_id target); + MICROSOFT_QUANTUM_DECL void MCAdjointAnd_cpp(unsigned sim_id,_In_ int length, _In_reads_(length) logical_qubit_id* controls, _In_ logical_qubit_id target); + + // rotations + MICROSOFT_QUANTUM_DECL void R_cpp(unsigned sim_id, _In_ int b, _In_ double phi, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void Rfrac_cpp(unsigned sim_id, _In_ int b, _In_ std::int64_t numerator, std::int64_t power, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void R1_cpp(unsigned sim_id, _In_ double phi, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void R1frac_cpp(unsigned sim_id, _In_ std::int64_t numerator, std::int64_t power, _In_ logical_qubit_id q); + + + + // multi-controlled rotations + MICROSOFT_QUANTUM_DECL void MCR_cpp( + unsigned sim_id, + _In_ int b, + _In_ double phi, + _In_ logical_qubit_id n, + _In_reads_(n) logical_qubit_id* c, + _In_ logical_qubit_id q); + + MICROSOFT_QUANTUM_DECL void MCRFrac_cpp( + unsigned sim_id, + _In_ int b, + _In_ std::int64_t numerator, + _In_ std::int64_t power, + _In_ logical_qubit_id nc, + _In_reads_(nc) logical_qubit_id* c, + _In_ logical_qubit_id q); + + MICROSOFT_QUANTUM_DECL void MCR1_cpp( + unsigned sim_id, + _In_ double phi, + _In_ int n, + _In_reads_(n) logical_qubit_id* c, + _In_ logical_qubit_id q); + + MICROSOFT_QUANTUM_DECL void MCR1Frac_cpp( + unsigned sim_id, + _In_ std::int64_t numerator, + _In_ std::int64_t power, + _In_ int nc, + _In_reads_(nc) logical_qubit_id* c, + _In_ logical_qubit_id q); + + // Exponential of Pauli operators + MICROSOFT_QUANTUM_DECL void Exp_cpp( + unsigned sim_id, + _In_ int n, + _In_reads_(n) int* b, + _In_ double phi, + _In_reads_(n) logical_qubit_id* q); + MICROSOFT_QUANTUM_DECL void MCExp_cpp( + unsigned sim_id, + _In_ int nc, + _In_ int n, + _In_reads_(nc) logical_qubit_id* c, + _In_reads_(n) int* b, + _In_ double phi, + _In_reads_(n) logical_qubit_id* q); + + // measurements + MICROSOFT_QUANTUM_DECL bool M_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL void Reset_cpp(unsigned sim_id, _In_ logical_qubit_id q); + MICROSOFT_QUANTUM_DECL bool Measure_cpp( + unsigned sim_id, + _In_ int n, + _In_reads_(n) int* b, + _In_reads_(n) logical_qubit_id* q); + + + MICROSOFT_QUANTUM_DECL double JointEnsembleProbability_cpp( + unsigned sim_id, + _In_ int n, + _In_reads_(n) int* b, + _In_reads_(n) logical_qubit_id* q); + + MICROSOFT_QUANTUM_DECL void Dump_cpp(unsigned sim_id, _In_ logical_qubit_id max_qubit_id, _In_ void (*callback)(char*, double, double)); + MICROSOFT_QUANTUM_DECL bool DumpQubits_cpp( + unsigned sim_id, + _In_ logical_qubit_id n, + _In_reads_(n) logical_qubit_id* q, + _In_ void (*callback)(char* , double, double)); + + MICROSOFT_QUANTUM_DECL void GetAmplitude_cpp(_In_ unsigned sim_id, _In_ size_t length, _In_ const char* label, _Out_ double &real, _Out_ double &imag); + + MICROSOFT_QUANTUM_DECL bool Assert_cpp(unsigned sim_id, _In_ int n, _In_reads_(n) int* b, _In_reads_(n) logical_qubit_id* q, bool result); + MICROSOFT_QUANTUM_DECL const char* Sample_cpp(unsigned sim_id); +} diff --git a/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp b/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp new file mode 100644 index 00000000000..4ba408b6e7b --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/factory.cpp @@ -0,0 +1,63 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +// Manages simulators in a vector of pointers to simulators + +#include +#include + +#include "factory.hpp" +#include "SparseSimulator.h" +#include "types.h" + +namespace Microsoft::Quantum::SPARSESIMULATOR +{ + +// Ensures exclusive access to _simulators, the vector of simulators +std::shared_mutex _mutex; +std::vector> _simulators; + +unsigned createSimulator(logical_qubit_id num_qubits) +{ + if (num_qubits > MAX_QUBITS) + throw std::runtime_error("Max number of qubits is exceeded!"); + std::lock_guard lock(_mutex); + size_t emptySlot = -1; + for (auto const& s : _simulators) + { + if (s == nullptr) + { + emptySlot = &s - &_simulators[0]; + break; + } + } + if (emptySlot == -1) + { + _simulators.push_back(std::make_shared(num_qubits)); + emptySlot = _simulators.size() - 1; + } + else + { + _simulators[emptySlot] = std::make_shared(num_qubits); + } + + return static_cast(emptySlot); +} + +// Deletes a simulator in the vector +void destroySimulator(unsigned id) +{ + std::lock_guard lock(_mutex); + + _simulators[id].reset(); // Set pointer to nullptr +} + +// Returns a simulator at some id (used for the C++/C# API) +std::shared_ptr& getSimulator(unsigned id) +{ + std::shared_lock shared_lock(_mutex); + + return _simulators[id]; +} + +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp b/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp new file mode 100644 index 00000000000..ffb54385f22 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/factory.hpp @@ -0,0 +1,19 @@ +// 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::Quantum::SPARSESIMULATOR +{ + +unsigned createSimulator(logical_qubit_id); +void destroySimulator(unsigned); + +std::shared_ptr& getSimulator(unsigned); + +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/gates.h b/src/Simulation/Simulators/SparseSimulator/Native/gates.h new file mode 100644 index 00000000000..e0a3d8a7b54 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/gates.h @@ -0,0 +1,279 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +namespace Microsoft::Quantum::SPARSESIMULATOR +{ + +// 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) {} + + 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) + {} +}; + +// 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 const& 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 const& 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 const& 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 const& 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 + +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp b/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp new file mode 100644 index 00000000000..c476f1f9361 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/quantum_state.hpp @@ -0,0 +1,1178 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "basic_quantum_state.hpp" + +#include "types.h" +#include "gates.h" + +using namespace std::literals::complex_literals; + +namespace Microsoft::Quantum::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; +} + +template +bool get_parity(std::bitset const& bitstring){ + 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(); +} + +// 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; + 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); }; + } + + // 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); + } + } + + 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 *= (get_parity(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 (!get_parity(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 (get_parity(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 * (get_parity(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 = get_parity(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 *= (get_parity(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 (!get_parity(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 (get_parity(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 * (get_parity(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 = get_parity(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) { + // The amplitude of current_state should always be non-zero, if the data structure + // is properly maintained. Since the flipped state should match the amplitude (up to phase), + // if the flipped state is not in _qubit_data, it implicitly has an ampltude of 0.0, which + // is *not* a match, so the assertion should fail. + auto flipped_state = _qubit_data.find(current_state->first ^ XYs); + if (flipped_state == _qubit_data.end() || + std::norm(flipped_state->second - current_state->second * (get_parity(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 * (get_parity(label & YZs) ? -phaseShift : phaseShift); + std::cout << ", got " << (flipped_state == _qubit_data.end() ? 0.0 : flipped_state->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)) * (get_parity(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 const& 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 const& 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 const& 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 + // TODO: Is this still needed after multithreading is removed? Can we work off operation_list? + std::vector operation_vector; + 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; + } + } + + wavefunction new_qubit_data = make_wavefunction(); + + // 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 (get_parity(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); + operation_vector.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 _qubit_data + 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; } + +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); + } + + // 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 const& 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; + } + + // 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 Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/Native/types.h b/src/Simulation/Simulators/SparseSimulator/Native/types.h new file mode 100644 index 00000000000..a9f6fb656c4 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/Native/types.h @@ -0,0 +1,36 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +#include +#include +#include +#include + +namespace Microsoft::Quantum::SPARSESIMULATOR +{ + +#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; + +using amplitude = std::complex; + +template +using qubit_label_type = std::bitset; + +// Wavefunctions are hash maps of some key (std::bitset or a string) +template +using abstract_wavefunction = std::unordered_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; + +} // namespace Microsoft::Quantum::SPARSESIMULATOR diff --git a/src/Simulation/Simulators/SparseSimulator/README.md b/src/Simulation/Simulators/SparseSimulator/README.md new file mode 100644 index 00000000000..a5d774569b6 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/README.md @@ -0,0 +1,144 @@ +# Introduction +This is a an alternative quantum simulator, compatible with Q\# and Microsoft's quantum development kit. It is efficient for algorithms with a smaller number of states in superposition, regardless of the number of qubits. + +# 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 + +# Dependencies + - CMake v1.2.0 + - Q\# v0.15.2101126940 + - C++11 + - Dotnet core v 3.1 + - Dotnet cake v1.1.0 + +# Setup +To build the Sparse Simulator, call `dotnet restore` from the main `SparseSimulator` folder, then `dotnet cake`. This builds the C++ backend and the C\# interface. + +To use the SparseSimulator in a Q\# project, ensure that it includes +`` +` ` +`` +in the `.csproj` file, to import the sparse simulator. + +The C\# sparse simulator class and all its Q\# functions are in the `Microsoft.Quantum.SparseSimulation` namespace. + +A sparse simulator object can be initialized in C\# by calling `SparseSimulator()`. In an executable Q\# project, add +`Microsoft.Quantum.SparseSimulation.SparseSimulator` +to the property group in the `.csproj` file. + +To run a test with the sparse simulator, add `@Test("Microsoft.Quantum.SparseSimulation.SparseSimulator")` above the test operation. + +# Other + +## Qubit Capacity +The simulator has a qubit capacity which can be either 64, 128, 256, 512, or 1024. The capacity is only related to the internal data structure of quantum states, and only needs to be at least as large as the number of qubits actually in use. By default it starts with 64 qubits, and allocating more than 64 qubits will simply increase its capacity to accommodate this. It can be initialized with more qubits with an optional argument to C\#, i.e., `var sim = SparseSimulator(128)` will initialize with capacity for 128 qubits. + +## Emulation +Currently the simulator only emulates AND and adjoint AND gates. It asserts that the target qubit is |0> before the AND and after the adjoint AND. + +## Statistics +By default, the simulator tracks the total number of qubits used. It also tracks the number of occupied qubits, which is based on whether they are ever in a non-zero state. That is, if one allocates a qubit and applies Z to it, the simulator will still regard it as unoccupied because it wil still be in the |0> state. + +## Assertion optimizations +Z-assertions can be treated as phase/permutation gates, and thus in a Release build, they are added to the queue and not immediately run. This means if such an assertion fails, it might not give an error until much later in the execution. In a Debug build, they are executed immediately (which will slow down execution, but be more helpful for debugging). + + +# Special Operations +The sparse simulator has some extra C\# Q\# operations for specific purposes. + +## Member functions of `SparseSimulator` +`SetLogging(bool)`: Logging is set to `false` by default; if set to true, the simulator will output the name of each operation as it executes them. + +## Q\# Functions + +`GetAmplitudeFromInt(qubit : Qubit[], label : Int)`: Uses `label` as a bitstring to index the qubits in `qubit`. It returns the amplitude of the state with that label, where the value of all other qubits is 0. + +`GetAmplitude(qubit : Qubit[], label : BigInt) : Complex` : Same as `GetAmplitudeFromInt`, but uses a BigInt as an index. + +`AssertProbBigInt (stateIndex : BigInt, expected : Double, qubits : Microsoft.Quantum.Arithmetic.LittleEndian, tolerance : Double) : Unit`: Asserts that the probability to measure the qubits in `Qubits` and find the value `stateIndex` is within `tolerance` of the expected probability `expected`. + +`Sample(register : qubit [])`: This acts like measuring `register` in the Pauli-Z basis, and returns a boolean array of the result of that measurement, except it is non-destructive. + +# Adding Gates + +To add a gate to the simulator, you will need to add it to: +- SparseSimulatorCS/SparseSimulator.cs +- Native/capi.cpp +- Native/quantum_state.hpp +- Native/basic_quantum_state.hpp +- Native/SparseSimulator.h + +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 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. + +Finally, the `BasicQuantumState` class needs to implement the same function virtually, so the `SparseSimulator` class can call it. + +## SparseSimulator (SparseSimulator.h) +The `SparseSimulator` class will also need a wrapper function for the new gate. This will call `_quantum_state->your_new_gate(args...);`. However, `SparseSimulator` also manages gate queues, so you will need to decide how the new gate commutes with existing ones. The simplest approach will be to start by calling `_execute_queued_ops();`, which will flush all the gate queues and allow you to run your gate without needing any commutation relations. + +## C API (capi.cpp) +Within capi.cpp, you need to add a new wrapper to call the function in `SparseSimulator`. The first argument must be an unsigned integer, representing the ID of the simulator, so you will call it as +`getSimulator(sim_id)->your_new_gate(\# args \#);` +For more complicated arguments from C\#, the other functions provide a template. For consistency, append "\_cpp" to the end of the function name in capi.cpp. + + +## SparseSimulator (SparseSimulator.cs) +Here the C\# code will call your new gate. The wrapper must go in the `SparseSimulatorProcessor` class, and must have the following form: +`private static extern void your_new_gate_cpp(uint sim, \# args\#);` +`public override void Your_new_gate(\#args\#) +{ + your_new_gate_cpp(Id, \# arguments, parsed into a format that C++ can read\#); +}` +Here `Id` is an internal variable for the simulator's ID, which must be passed to C++. + +## (Optional) Q\# Code +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: +```C# +public partial class YourNewGate +{ + public class Native : YourNewGate + { + private SparseSimulator sim = null; + public Native(IOperationFactory m) : base(m) + { + sim = m as SparseSimulator; + } + 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 +```c# +public partial class SparseSimulator : QuantumProcessorDispatcher, IDisposable +{ + public 'Q# return type' 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. + + +# Internal Logic +We desribe the main data structure strategies and optimizations in the paper at https://arxiv.org/abs/2105.01533. + +# Future Optimizations + +## Delayed Release +Currently the simulator executes any queued operations on any qubits it needs to release. This is not strictly necessary: it could add an assertion that the qubits are zero, and continue on. However, it's important to set the `_occupied_qubits` vector to be `0` for that qubit after it is released, but this will not necessarily be true if the gates are not executed. Hence, to delay the release of qubits, it will need a more involved method to track occupied qubits, and qubits which are actually 0. diff --git a/src/Simulation/Simulators/SparseSimulator/SparseQuantumSimulator.sln b/src/Simulation/Simulators/SparseSimulator/SparseQuantumSimulator.sln new file mode 100644 index 00000000000..05a785cd491 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseQuantumSimulator.sln @@ -0,0 +1,65 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio Version 16 +VisualStudioVersion = 16.0.30907.101 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "SparseSimulatorTests", "SparseSimulatorTests\SparseSimulatorTests.vcxproj", "{E1663845-BBA8-41EA-AC31-50A092EEA728}" +EndProject +Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "SparseSimulator", "SparseSimulatorCS\SparseSimulator.csproj", "{1E901047-6E72-4B99-8683-A5B64D755400}" +EndProject +Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{1CECBF17-34F0-41FF-BE04-88659415074A}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "SparseSimQSharpTests", "SparseSimQSharpTests\SparseSimQSharpTests.csproj", "{2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Any CPU = Debug|Any CPU + Debug|x64 = Debug|x64 + Debug|x86 = Debug|x86 + Release|Any CPU = Release|Any CPU + Release|x64 = Release|x64 + Release|x86 = Release|x86 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Debug|Any CPU.ActiveCfg = Debug|Win32 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Debug|x64.ActiveCfg = Debug|x64 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Debug|x64.Build.0 = Debug|x64 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Debug|x86.ActiveCfg = Debug|Win32 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Debug|x86.Build.0 = Debug|Win32 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Release|Any CPU.ActiveCfg = Release|Win32 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Release|x64.ActiveCfg = Release|x64 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Release|x64.Build.0 = Release|x64 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Release|x86.ActiveCfg = Release|Win32 + {E1663845-BBA8-41EA-AC31-50A092EEA728}.Release|x86.Build.0 = Release|Win32 + {1E901047-6E72-4B99-8683-A5B64D755400}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Debug|Any CPU.Build.0 = Debug|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Debug|x64.ActiveCfg = Debug|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Debug|x64.Build.0 = Debug|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Debug|x86.ActiveCfg = Debug|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Debug|x86.Build.0 = Debug|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Release|Any CPU.ActiveCfg = Release|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Release|Any CPU.Build.0 = Release|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Release|x64.ActiveCfg = Release|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Release|x64.Build.0 = Release|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Release|x86.ActiveCfg = Release|Any CPU + {1E901047-6E72-4B99-8683-A5B64D755400}.Release|x86.Build.0 = Release|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Debug|Any CPU.Build.0 = Debug|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Debug|x64.ActiveCfg = Debug|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Debug|x64.Build.0 = Debug|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Debug|x86.ActiveCfg = Debug|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Debug|x86.Build.0 = Debug|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Release|Any CPU.ActiveCfg = Release|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Release|Any CPU.Build.0 = Release|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Release|x64.ActiveCfg = Release|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Release|x64.Build.0 = Release|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Release|x86.ActiveCfg = Release|Any CPU + {2D6F4915-A10D-4109-8B83-DA6E5C36AEA4}.Release|x86.Build.0 = Release|Any CPU + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection + GlobalSection(ExtensibilityGlobals) = postSolution + SolutionGuid = {756B5CCE-BACE-4478-A9A1-25027C14669B} + EndGlobalSection +EndGlobal diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs b/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs new file mode 100644 index 00000000000..c0886fe7ca1 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/Program.qs @@ -0,0 +1,631 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.SparseSimulatorTests { + + open Microsoft.Quantum.Canon; + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Diagnostics; + open Microsoft.Quantum.Convert; + open Microsoft.Quantum.Arithmetic; + open Microsoft.Quantum.Arrays; + open Microsoft.Quantum.Preparation; + open Microsoft.Quantum.SparseSimulation; + open Microsoft.Quantum.Math; + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation ControlledHTest() : Unit { + use qubits = Qubit[4]{ + H(qubits[0]); + H(qubits[1]); + CNOT(qubits[0], qubits[2]); + CNOT(qubits[1], qubits[3]); + (Controlled H)([qubits[2]], (qubits[3])); + let expected_vals = [ + Complex(0.5, 0.0), + Complex(0.353553,0.0), + Complex(0.5,0.0), + Complex(0.353553,0.0), + Complex(0.353553,0.0), + Complex(-0.353553,0.0)]; + let expected_labels = [ + 0L, + 5L, + 10L, + 7L, + 13L, + 15L]; + AssertAmplitudes(expected_labels, expected_vals, LittleEndian(qubits), 0.001); + ResetAll(qubits); + } + } + + + operation DumpIdRotation() : Unit { + use qubits = Qubit[2] { + H(qubits[0]); + H(qubits[1]); + (Controlled Exp)([qubits[0]], ([PauliI], 0.5, [qubits[1]])); + DumpMachine(); + ResetAll(qubits); + } + } + + operation _R(pauli : Pauli, theta : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + R(pauli, theta, qubits[0]); + } + operation _MCR(pauli : Pauli, theta : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + (Controlled R)(qubits[1..Length(qubits)-1], (pauli, theta, qubits[0])); + } + operation _MCExp(pauli : Pauli, theta : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + (Controlled Exp)(qubits[1..Length(qubits)-1], ([pauli], theta, [qubits[0]])); + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation ExpRCompare() : Unit { + for angle in 0..2*314 { + AssertOperationsEqualReferenced(1, _R(PauliX, -IntAsDouble(angle)/100.0, _), Exp([PauliX], IntAsDouble(angle)/200.0, _)); + AssertOperationsEqualReferenced(1, _R(PauliZ, -IntAsDouble(angle)/100.0, _), Exp([PauliZ], IntAsDouble(angle)/200.0, _)); + AssertOperationsEqualReferenced(1, _R(PauliY, -IntAsDouble(angle)/100.0, _), Exp([PauliY], IntAsDouble(angle)/200.0, _)); + AssertOperationsEqualReferenced(1, _R(PauliI, -IntAsDouble(angle)/100.0, _), Exp([PauliI], IntAsDouble(angle)/200.0, _)); + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation MCExpMCRCompare() : Unit { + for angle in 0..2*314 { + AssertOperationsEqualReferenced(1, _MCR(PauliX, -IntAsDouble(angle)/100.0, _), _MCExp(PauliX, IntAsDouble(angle)/200.0, _)); + AssertOperationsEqualReferenced(1, _MCR(PauliZ, -IntAsDouble(angle)/100.0, _), _MCExp(PauliZ, IntAsDouble(angle)/200.0, _)); + AssertOperationsEqualReferenced(1, _MCR(PauliY, -IntAsDouble(angle)/100.0, _), _MCExp(PauliY, IntAsDouble(angle)/200.0, _)); + AssertOperationsEqualReferenced(1, _MCR(PauliI, -IntAsDouble(angle)/100.0, _), _MCExp(PauliI, IntAsDouble(angle)/200.0, _)); + } + } + + operation DumpMultiplexZ() : Unit { + use qubits = Qubit[12] { + ApplyToEach(H, qubits[0..5]); + ApplyToEach(CNOT, Zipped(qubits[0..5], qubits[6..11])); + ApproximatelyMultiplexZ(0.001, [0.12, 0.34, -0.26, 0.5, 1.8], LittleEndian(qubits[6..10]), qubits[11]); + DumpMachine(); + ResetAll(qubits); + } + } + + + + internal operation ControlledRz(angle : Double, (control : Qubit, target : Qubit)) : Unit is Adj { + Controlled Rz([control], (angle, target)); + DumpMachine(); + } + + internal operation ControlledRzAsR1(angle : Double, (control : Qubit, target : Qubit)) : Unit is Adj { + Controlled R1([control], (angle, target)); + R1(-angle / 2.0, control); + DumpMachine(); + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation TestEqualityOfControlledRz() : Unit { + for _ in 1..10 { + let angle = Microsoft.Quantum.Random.DrawRandomDouble(0.0, 2.0 * PI()); + AssertOperationsEqualReferenced(2, ApplyToFirstTwoQubits(ControlledRzAsR1(angle, _), _), ApplyToFirstTwoQubitsA(ControlledRz(angle, _), _)); + } + } + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation LargeStateTests() : Unit { + let nqubits = 12; + LargeStateTestWrapper(Rotation1CompareTest, nqubits); + LargeStateTestWrapper(RotationCompareTest, nqubits); + LargeStateTestWrapper(SWAPTest, nqubits); + LargeStateTestWrapper(CSWAPTest, nqubits); + LargeStateTestWrapper(CNOTTest, nqubits); + LargeStateTestWrapper(ResetTest, nqubits); + LargeStateTestWrapper(AssertTest, nqubits); + LargeStateTestWrapper(AndChainTest, nqubits); + LargeStateTestWrapper(CZTest, nqubits); + LargeStateTestWrapper(AllocationTest, nqubits); + LargeStateTestWrapper(Rotation1CompareTest, nqubits); + LargeStateTestWrapper(RotationFracCompareTest, nqubits); + + } + + operation _ToffoliCSwap(targets : Qubit[]) : Unit is Ctl + Adj { + CCNOT(targets[0], targets[1], targets[2]); + CCNOT(targets[0], targets[2], targets[1]); + CCNOT(targets[0], targets[1], targets[2]); + } + operation _ArrayCSwap(targets : Qubit[]) : Unit is Ctl + Adj { + (Controlled SWAP)([targets[0]], (targets[1], targets[2])); + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation CCNOTvToffoliTest() : Unit { + AssertOperationsEqualReferenced(3, _ArrayCSwap, _ToffoliCSwap); + } + + // Creates a big state so an existing test forces parallel execution + operation LargeStateTestWrapper(test : (Unit => Unit), nqubits : Int) : Unit { + let prob = PowD(0.5, IntAsDouble(nqubits)); + use qubits = Qubit[nqubits]; + for idx in 0..nqubits - 1 { + H(qubits[idx]); + } + for idy in 1..128 { + CCNOT(qubits[idy % nqubits], qubits[(idy+1) % nqubits], qubits[(idy+2) % nqubits]); + } + test(); + ResetAll(qubits); + } + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation PartialDumpTest() : Unit { + use qubits = Qubit[4] { + ApplyToEach(H, qubits); + CNOT(qubits[2], qubits[3]); + DumpRegister("Test_file_1", qubits[0..1]); + DumpRegister("Test_file_2", qubits[2..3]); + CNOT(qubits[2], qubits[3]); + H(qubits[2]); + CNOT(qubits[0], qubits[2]); + DumpRegister("Test_file_3", qubits[0..1]); + DumpRegister("Test_file_4", qubits); + ResetAll(qubits); + } + } + + + internal operation FakeR1Frac(numerator : Int, denominator : Int, qubit : Qubit[]) : Unit is Adj + Ctl { + RFrac(PauliZ, -numerator, denominator + 1, qubit[0]); + RFrac(PauliI, numerator, denominator + 1, qubit[0]); + } + + internal operation R1FracWithArray(numerator : Int, denominator : Int, qubit : Qubit[]) : Unit is Adj + Ctl { + R1Frac(numerator, denominator, qubit[0]); + } + + internal operation FakeR1(angle : Double, qubit : Qubit[]) : Unit is Adj + Ctl { + R(PauliZ, angle, qubit[0]); + R(PauliI, -angle, qubit[0]); + } + internal operation R1WithArray(angle : Double, qubit : Qubit[]) : Unit is Adj + Ctl { + R1(angle, qubit[0]); + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation Rotation1CompareTest() : Unit { + for denom in 0..5{ + for num in 1..2..(2^denom - 1){ + AssertOperationsEqualReferenced(1, R1FracWithArray(num, denom, _), FakeR1Frac(num, denom, _)); + AssertOperationsEqualReferenced(1, FakeR1Frac(num, denom, _), R1FracWithArray(num, denom, _)); + } + } + for angle in 0..314 { + AssertOperationsEqualReferenced(1, R1WithArray(IntAsDouble(angle)/100.0, _), FakeR1(IntAsDouble(angle)/100.0, _)); + AssertOperationsEqualReferenced(1, FakeR1(IntAsDouble(angle)/100.0, _), R1WithArray(IntAsDouble(angle)/100.0, _)); + } + } + + internal operation RFracWithArray(axis : Pauli, num : Int, denom : Int, qubit : Qubit[]) : Unit is Adj + Ctl { + RFrac(axis, num, denom, qubit[0]); + } + internal operation RWithArray(axis : Pauli, angle : Double, qubit : Qubit[]) : Unit is Adj + Ctl { + R(axis, angle, qubit[0]); + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation RotationFracCompareTest() : Unit { + for denom in 0..5{ + for num in 1..2..(2^denom - 1){ + let angle = -3.14159265*IntAsDouble(num)/(2.0 ^IntAsDouble(denom-1)); + AssertOperationsEqualReferenced(1, RFracWithArray(PauliX, num, denom, _), RWithArray(PauliX, angle, _)); + AssertOperationsEqualReferenced(1, RWithArray(PauliX, angle, _), RFracWithArray(PauliX, num, denom, _)); + AssertOperationsEqualReferenced(1, RFracWithArray(PauliY, num, denom, _), RWithArray(PauliY, angle, _)); + AssertOperationsEqualReferenced(1, RWithArray(PauliY, angle, _), RFracWithArray(PauliY, num, denom, _)); + AssertOperationsEqualReferenced(1, RFracWithArray(PauliZ, num, denom, _), RWithArray(PauliZ, angle, _)); + AssertOperationsEqualReferenced(1, RWithArray(PauliZ, angle, _), RFracWithArray(PauliZ, num, denom, _)); + } + } + } + + operation _HadamardByRotations(qubit : Qubit[]) : Unit is Adj + Ctl { + RFrac(PauliY, -1, 2, qubit[0]); + RFrac(PauliX, -1, 1, qubit[0]); + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation RotationCompareTest() : Unit { + AssertOperationsEqualReferenced(1, _HadamardByRotations(_), ApplyToEachCA(H, _)); + AssertOperationsEqualReferenced(1, ApplyToEachCA(R1Frac(1,0,_), _), ApplyToEachCA(Z, _)); + AssertOperationsEqualReferenced(1, ApplyToEachCA(R1Frac(1,1,_), _), ApplyToEachCA(S, _)); + AssertOperationsEqualReferenced(1, ApplyToEachCA(R1Frac(1,2,_), _), ApplyToEachCA(T, _)); + AssertOperationsEqualReferenced(1, ApplyToEachCA(RFrac(PauliX, 1, 1, _), _), ApplyToEachCA(X, _)); + AssertOperationsEqualReferenced(1, ApplyToEachCA(RFrac(PauliY, 1, 1, _), _), ApplyToEachCA(Y, _)); + AssertOperationsEqualReferenced(1, ApplyToEachCA(RFrac(PauliZ, 1, 1, _), _), ApplyToEachCA(Z, _)); + } + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation AllocationTest() : Unit { + use qubits = Qubit[512]{ + for idx in 0..511 { + X(qubits[idx]); + } + ResetAll(qubits); + } + } + + // Taken from the PurifiedMixedState documentation + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation QROMPrepareTest() : Unit { + let coefficients = [1.0, 2.0, 3.0, 4.0, 5.0]; + let targetError = 1e-3; + let purifiedState = PurifiedMixedState(targetError, coefficients); + use indexRegister = Qubit[purifiedState::Requirements::NIndexQubits] { + use garbageRegister = Qubit[purifiedState::Requirements::NGarbageQubits] { + purifiedState::Prepare(LittleEndian(indexRegister), new Qubit[0], garbageRegister); + ResetAll(garbageRegister); + } + ResetAll(indexRegister); + } + } + + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation CZTest() : Unit { + let num_qubits = 5; + use qubits = Qubit[num_qubits]{ + H(qubits[0]); + for idx in 0..(2^(num_qubits-1) - 1) { + let result = (idx == 2^(num_qubits - 1) - 1); + within { + ApplyXorInPlace(idx, LittleEndian(qubits[1..num_qubits - 1])); + (Controlled Z)(qubits[1..num_qubits -1], (qubits[0])); + } apply { + if (result){ + AssertMeasurement([PauliX], qubits[0..0], One, "CZ failed to add phase"); + } else { + AssertMeasurement([PauliX], qubits[0..0], Zero, "CZ added unexpected phase"); + } + } + } + H(qubits[0]); + } + } + + operation ApplyAndChain(andOp : ((Qubit, Qubit, Qubit)=>Unit is Adj + Ctl), auxRegister : Qubit[], ctrlRegister : Qubit[], target : Qubit) + : Unit is Adj { + if (Length(ctrlRegister) == 0) { + X(target); + } elif (Length(ctrlRegister) == 1) { + CNOT(Head(ctrlRegister), target); + } else { + EqualityFactI(Length(auxRegister), Length(ctrlRegister) - 2, "Unexpected number of auxiliary qubits"); + let controls1 = ctrlRegister[0..0] + auxRegister; + let controls2 = Rest(ctrlRegister); + let targets = auxRegister + [target]; + ApplyToEachA(andOp, Zipped3(controls1, controls2, targets)); + } + } + + operation AndChainDump() : Unit { + let num_qubits = 5; + use qubits = Qubit[num_qubits]{ + use aux = Qubit[num_qubits - 3]{ + within{ + for idx in 1..num_qubits - 1 { + H(qubits[idx]); + if (idx % 3 == 0){ Z(qubits[idx]);} + } + ApplyAndChain(ApplyAnd, aux, qubits[1..num_qubits -1], qubits[0]); + } apply { + DumpMachine(); + } + } + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation AndChainTest() : Unit { + let num_qubits = 5; + use qubits = Qubit[num_qubits]{ + use aux = Qubit[num_qubits - 3]{ + for idx in 0..(2^(num_qubits-1) - 1) { + let result = (idx == 2^(num_qubits - 1) - 1); + within { + ApplyXorInPlace(idx, LittleEndian(qubits[1..num_qubits - 1])); + ApplyAndChain(ApplyAnd, aux, qubits[1..num_qubits -1], qubits[0]); + } apply { + let after = M(qubits[0]); + if (result){ + Fact(after == One, "Did not apply AND"); + } else { + Fact(after == Zero, "Applied AND unexpectedly"); + } + } + } + } + } + } + + + operation DumpMCXFrac() : Unit { + use qubits = Qubit[2] { + H(qubits[0]); + H(qubits[1]); + (Controlled RFrac)([qubits[0]], (PauliX, -1, 3, qubits[1])); + DumpMachine(); + (Adjoint Controlled RFrac)([qubits[0]], (PauliX, -1, 3, qubits[1])); + H(qubits[0]); + H(qubits[1]); + H(qubits[0]); + H(qubits[1]); + (Controlled R)([qubits[0]], (PauliX, 0.25*Microsoft.Quantum.Math.PI(), qubits[1])); + DumpMachine(); + (Adjoint Controlled R)([qubits[0]], (PauliX, 0.25*Microsoft.Quantum.Math.PI(), qubits[1])); + H(qubits[0]); + H(qubits[1]); + } + } + + operation DumpMCZFrac() : Unit { + use qubits = Qubit[2] { + H(qubits[0]); + H(qubits[1]); + (Controlled RFrac)([qubits[0]], (PauliZ, -1, 3, qubits[1])); + DumpMachine(); + (Adjoint Controlled RFrac)([qubits[0]], (PauliZ, -1, 3, qubits[1])); + H(qubits[0]); + H(qubits[1]); + H(qubits[0]); + H(qubits[1]); + (Controlled R)([qubits[0]], (PauliZ, 0.25*Microsoft.Quantum.Math.PI(), qubits[1])); + DumpMachine(); + (Adjoint Controlled R)([qubits[0]], (PauliZ, 0.25*Microsoft.Quantum.Math.PI(), qubits[1])); + H(qubits[0]); + H(qubits[1]); + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation AssertProbTest() : Unit { + let tolerance = 0.000001; + use qubit = Qubit() { + AssertMeasurementProbability([PauliZ], [qubit], Zero, 1.0, "Failed assert Z on |0>", tolerance); + AssertMeasurementProbability([PauliX], [qubit], Zero, 0.5, "Failed assert X on |0>", tolerance); + AssertMeasurementProbability([PauliY], [qubit], Zero, 0.5, "Failed assert Y on |0>", tolerance); + H(qubit); + AssertMeasurementProbability([PauliZ], [qubit], Zero, 0.5, "Failed assert Z on |+>", tolerance); + AssertMeasurementProbability([PauliX], [qubit], Zero, 1.0, "Failed assert X on |+>", tolerance); + AssertMeasurementProbability([PauliY], [qubit], Zero, 0.5, "Failed assert X on |+>", tolerance); + S(qubit); + AssertMeasurementProbability([PauliZ], [qubit], Zero, 0.5, "Failed assert Z on |+>", tolerance); + AssertMeasurementProbability([PauliX], [qubit], Zero, 0.5, "Failed assert X on |+>", tolerance); + AssertMeasurementProbability([PauliY], [qubit], Zero, 1.0, "Failed assert Y on |i>", tolerance); + S(qubit); + use buddy = Qubit() { + CNOT(qubit, buddy); + AssertMeasurementProbability([PauliZ, PauliZ], [qubit, buddy], Zero, 1.0, "Failed assert ZZ on |++>", tolerance); + AssertMeasurementProbability([PauliZ, PauliI], [qubit, buddy], Zero, 0.5, "Failed assert ZI on |++>", tolerance); + AssertMeasurementProbability([PauliX, PauliX], [qubit, buddy], Zero, 0.0, "Failed assert XX on |++>", tolerance); + AssertMeasurementProbability([PauliX, PauliI], [qubit, buddy], Zero, 0.5, "Failed assert XI on |++>", tolerance); + AssertMeasurementProbability([PauliY, PauliY], [qubit, buddy], Zero, 1.0, "Failed assert YY on |++>", tolerance); + AssertMeasurementProbability([PauliY, PauliI], [qubit, buddy], Zero, 0.5, "Failed assert YI on |++>", tolerance); + AssertMeasurementProbability([PauliX, PauliZ], [qubit, buddy], Zero, 0.5, "Failed assert XZ on |++>", tolerance); + AssertMeasurementProbability([PauliX, PauliY], [qubit, buddy], Zero, 0.5, "Failed assert XY on |++>", tolerance); + AssertMeasurementProbability([PauliY, PauliZ], [qubit, buddy], Zero, 0.5, "Failed assert YZ on |++>", tolerance); + Z(qubit); + AssertMeasurementProbability([PauliZ, PauliZ], [qubit, buddy], Zero, 1.0, "Failed assert ZZ on |-->", tolerance); + AssertMeasurementProbability([PauliZ, PauliI], [qubit, buddy], Zero, 0.5, "Failed assert ZI on |-->", tolerance); + AssertMeasurementProbability([PauliX, PauliX], [qubit, buddy], Zero, 1.0, "Failed assert XX on |-->", tolerance); + AssertMeasurementProbability([PauliX, PauliI], [qubit, buddy], Zero, 0.5, "Failed assert XI on |-->", tolerance); + AssertMeasurementProbability([PauliY, PauliY], [qubit, buddy], Zero, 0.0, "Failed assert YY on |-->", tolerance); + AssertMeasurementProbability([PauliY, PauliI], [qubit, buddy], Zero, 0.5, "Failed assert YI on |-->", tolerance); + AssertMeasurementProbability([PauliX, PauliZ], [qubit, buddy], Zero, 0.5, "Failed assert XZ on |-->", tolerance); + AssertMeasurementProbability([PauliX, PauliY], [qubit, buddy], Zero, 0.5, "Failed assert XY on |-->", tolerance); + AssertMeasurementProbability([PauliY, PauliZ], [qubit, buddy], Zero, 0.5, "Failed assert YZ on |-->", tolerance); + Reset(buddy); + } + Reset(qubit); + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation AssertTest() : Unit { + use qubits = Qubit[3] { + AssertMeasurement([PauliZ, PauliZ, PauliZ], qubits, Zero, "Assert fails Pauli Z"); + X(qubits[0]); + X(qubits[1]); + X(qubits[2]); + AssertMeasurement([PauliZ, PauliZ, PauliZ], qubits, One, "Assert fails Pauli Z"); + AssertMeasurement([PauliZ, PauliZ, PauliI], qubits, Zero, "Assert fails Pauli Z"); + X(qubits[2]); + X(qubits[1]); + X(qubits[0]); + H(qubits[0]); + CNOT(qubits[0], qubits[1]); + CNOT(qubits[0], qubits[2]); + AssertMeasurement([PauliX, PauliX, PauliX], qubits, Zero, "Assert fails Pauli X"); + Z(qubits[0]); + AssertMeasurement([PauliX, PauliX, PauliX], qubits, One, "Assert fails Pauli X"); + S(qubits[0]); + AssertMeasurement([PauliY, PauliY, PauliY], qubits, Zero, "Assert fails Pauli Y"); + Z(qubits[0]); + AssertMeasurement([PauliY, PauliY, PauliY], qubits, One, "Assert fails Pauli Y"); + ResetAll(qubits); + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation MTest() : Unit { + use qubit = Qubit() { + X(qubit); + let test = M(qubit); + Fact(M(qubit) == test, "M does not preserves state"); + X(qubit); + let test2 = M(qubit); + Fact(M(qubit) == test2, "M does not preserve state"); + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation ResetTest() : Unit { + use qubit = Qubit() { + Reset(qubit); + Fact(M(qubit) == Zero, "Failed reset"); + X(qubit); + Reset(qubit); + Fact(M(qubit) == Zero, "Failed reset"); + use buddy = Qubit() { + X(qubit); + CNOT(qubit, buddy); + Reset(qubit); + Fact(M(qubit) == Zero, "Failed entangled reset"); + Fact(M(buddy) == One, "Failed entangled reset"); + Reset(buddy); + Fact(M(buddy) == Zero, "Failed entangled reset"); + } + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation SWAPTest () : Unit { + use qubits = Qubit[2] { + SWAP(qubits[0], qubits[1]); + Fact(M(qubits[0]) == Zero, "Bad swap on 00"); + Fact(M(qubits[1]) == Zero, "Bad swap on 00"); + X(qubits[0]); + SWAP(qubits[0], qubits[1]); + Fact(M(qubits[0]) == Zero, "Bad swap on 01"); + Fact(M(qubits[1]) == One, "Bad swap on 01"); + X(qubits[0]); + SWAP(qubits[0], qubits[1]); + Fact(M(qubits[0]) == One, "Bad swap on 11"); + Fact(M(qubits[1]) == One, "Bad swap on 11"); + X(qubits[0]); + SWAP(qubits[0], qubits[1]); + Fact(M(qubits[0]) == One, "Bad swap on 10"); + Fact(M(qubits[1]) == Zero, "Bad swap on 10"); + X(qubits[0]); + } + } + + operation AreAllQubitsOne(qubits : Qubit[]) : Bool { + mutable qubits_all_true = true; + for idx in 0..Length(qubits) - 1 { + if (M(qubits[idx]) == Zero){ + set qubits_all_true = false; + } + } + return qubits_all_true; + } + + operation CNOTTestInternal(target : Qubit) : Unit { + body (...){ + + } + controlled (controls, ...){ + if (AreAllQubitsOne(controls)){ + (Controlled X)(controls, target); + Fact(M(target) == One, "Bad CNOT"); + (Controlled X)(controls, target); + Fact(M(target) == Zero, "Bad CNOT"); + } else { + (Controlled X)(controls, target); + Fact(M(target) == Zero, "Bad CNOT"); + X(target); + (Controlled X)(controls, target); + Fact(M(target) == One, "Bad CNOT"); + X(target); + } + } + } + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation CNOTTest() : Unit { + use qubits = Qubit[3] { + (Controlled CNOTTestInternal)(qubits[0..1], qubits[2]); + X(qubits[0]); + (Controlled CNOTTestInternal)(qubits[0..1], qubits[2]); + X(qubits[1]); + (Controlled CNOTTestInternal)(qubits[0..1], qubits[2]); + X(qubits[0]); + (Controlled CNOTTestInternal)(qubits[0..1], qubits[2]); + X(qubits[1]); + } + } + + operation CSwapTestInternal (qubits : Qubit[]) : Unit { + body (...){ + // nothing + } + controlled (controls, ...) { + if (not AreAllQubitsOne(controls)) { + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == Zero, "Bad swap on 00"); + Fact(M(qubits[1]) == Zero, "Bad swap on 00"); + X(qubits[0]); + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == One, "Bad swap on 01"); + Fact(M(qubits[1]) == Zero, "Bad swap on 01"); + X(qubits[1]); + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == One, "Bad swap on 11"); + Fact(M(qubits[1]) == One, "Bad swap on 11"); + X(qubits[0]); + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == Zero, "Bad swap on 10"); + Fact(M(qubits[1]) == One, "Bad swap on 10"); + X(qubits[1]); + } else { + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == Zero, "Bad swap on 00"); + Fact(M(qubits[1]) == Zero, "Bad swap on 00"); + X(qubits[0]); + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == Zero, "Bad swap on 01"); + Fact(M(qubits[1]) == One, "Bad swap on 01"); + X(qubits[0]); + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == One, "Bad swap on 11"); + Fact(M(qubits[1]) == One, "Bad swap on 11"); + X(qubits[0]); + (Controlled SWAP)(controls, (qubits[0], qubits[1])); + Fact(M(qubits[0]) == One, "Bad swap on 10"); + Fact(M(qubits[1]) == Zero, "Bad swap on 10"); + X(qubits[0]); + } + } + } + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation CSWAPTest () : Unit { + use qubits = Qubit[4] { + (Controlled CSwapTestInternal)(qubits[0..1], qubits[2..3]); + X(qubits[0]); + (Controlled CSwapTestInternal)(qubits[0..1], qubits[2..3]); + X(qubits[1]); + (Controlled CSwapTestInternal)(qubits[0..1], qubits[2..3]); + X(qubits[0]); + (Controlled CSwapTestInternal)(qubits[0..1], qubits[2..3]); + X(qubits[1]); + } + use qubits = Qubit[7] { + (Controlled CSwapTestInternal)([qubits[0]], [qubits[6], qubits[3]]); + X(qubits[0]); + (Controlled CSwapTestInternal)([qubits[0]], [qubits[6], qubits[3]]); + X(qubits[0]); + } + } + + @Test("Microsoft.Quantum.SparseSimulation.SparseSimulator") + operation BasicTest () : Unit { + use new_qubit = Qubit() { + H(new_qubit); + X(new_qubit); + Z(new_qubit); + Y(new_qubit); + T(new_qubit); + (Adjoint T)(new_qubit); + S(new_qubit); + (Adjoint S)(new_qubit); + Y(new_qubit); + Z(new_qubit); + X(new_qubit); + H(new_qubit); + } + use new_qubits = Qubit[4] { + for idx in 0..3 { + H(new_qubits[idx]); + } + for idx in 0..3 { + H(new_qubits[idx]); + } + } + } +} diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/SparseSimQSharpTests.csproj b/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/SparseSimQSharpTests.csproj new file mode 100644 index 00000000000..4b8e18a115c --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimQSharpTests/SparseSimQSharpTests.csproj @@ -0,0 +1,15 @@ + + + + Exe + netcoreapp3.1 + Microsoft.Quantum.SparseSimulation.SparseSimulator + + + + + + + + + diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs new file mode 100644 index 00000000000..e8d54d0bffc --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.cs @@ -0,0 +1,360 @@ +// 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; +using Microsoft.Quantum.Simulation.Core; +using Microsoft.Quantum.Simulation.QuantumProcessor; +using System; +using System.IO; +using System.Linq; +using System.Runtime.InteropServices; +using System.Text; + +namespace Microsoft.Quantum.SparseSimulation +{ + + partial class GetAmplitudeFromInt + { + public class Native : GetAmplitudeFromInt + { + private SparseSimulator sim = null; + public Native(IOperationFactory m) : base(m) + { + sim = m as SparseSimulator; + } + + public override Func<(IQArray, long), Microsoft.Quantum.Math.Complex> __Body__ => (args) => { + return sim.GetAmplitude(args.Item1, args.Item2); + }; + + } + } + + partial class GetAmplitude + { + public class Native : GetAmplitude + { + private SparseSimulator sim = null; + public Native(IOperationFactory m) : base(m) + { + sim = m as SparseSimulator; + } + + public override Func<(IQArray, System.Numerics.BigInteger), Microsoft.Quantum.Math.Complex> __Body__ => args => { + return sim.GetAmplitude(args.Item1, args.Item2); + }; + + } + } + public partial class SparseSimulator : QuantumProcessorDispatcher, IDisposable + { + // // Converts a big integer into a string label + // where the label is read over the qubits passed as an argument, in order + // The label for all other qubits is assumed to be 0 + // That is: Suppose we have 11 = 1011 on qubits 0, 6,9,3. + // That will output the string 0001001001 + public static string InterleaveLabelByArray(IQArray qubits, System.Numerics.BigInteger label) + { + // Converts a big integer into a string label + // where the label is read over the qubits passed as an argument, in order + // The label for all other qubits is assumed to be 0 + char[] newlabel = "".PadRight((int)qubits.GetIds().Max() + 1, '0').ToCharArray(); + string test = new string(newlabel); + for (int idx = 0; idx < (int)qubits.Length; idx++) + { + if (((label >> idx) & 1) == 1) + { + int id = qubits[idx].Id; + newlabel[id] = '1'; + } + } + // std::bitset parses strings backwards + Array.Reverse(newlabel); + return new string(newlabel); + } + + // Returns the amplitude of a single state + public Microsoft.Quantum.Math.Complex GetAmplitude(string label) + { + return ((SparseSimulatorProcessor)this.QuantumProcessor).GetAmplitude(label); + } + public Microsoft.Quantum.Math.Complex GetAmplitude(uint label) + { + return GetAmplitude(label.ToString()); + } + public Microsoft.Quantum.Math.Complex GetAmplitude(System.Numerics.BigInteger label) + { + return GetAmplitude(label.ToString()); + } + + public Microsoft.Quantum.Math.Complex GetAmplitude(IQArray qubits, System.Numerics.BigInteger label) + { + return GetAmplitude(InterleaveLabelByArray(qubits, label)); + } + + + // Dumps the state in qubits to `target` (a filename), if the qubits are + // either null or not entangled to the rest of the machine; + // prints to console if `target` is null + protected virtual QVoid Dump(T target, IQArray qubits = null) + { + + var filename = (target is QVoid) ? "" : target.ToString(); + var logMessage = this.Get, Microsoft.Quantum.Intrinsic.Message>(); + // Once it turns `target` into a string action, it calls this process to push output + // into the string channel + QVoid process(Action channel) + { + // Just to output the ids + var ids = qubits?.Select(q => q.Id).ToArray() ?? this.QubitManager.AllocatedIds().Select(q => (int)q).ToArray(); + + // This will get passed, more or less, to the C++ + SparseSimulatorProcessor.StringCallback callback = (string label, double real, double imag) => + { + channel(Format(label, real, imag)); + }; + + channel($"# wave function for qubits with ids (least to most significant): {string.Join(";", ids)}"); + + // Calls the internal SparseSimulatorProcessor, which calls C++. + // Returns false if the state was entangled. + if (!((SparseSimulatorProcessor)this.QuantumProcessor).Dump(callback, ids.Max(), qubits)) + { + channel("## Qubits were entangled with an external qubit. Cannot dump corresponding wave function. ##"); + } + + return QVoid.Instance; + } + + // This generates the channel, then calls it + // If no file provided, use `Message` to generate the message into the console; + if (string.IsNullOrWhiteSpace(filename)) + { + var op = this.Get, Microsoft.Quantum.Intrinsic.Message>(); + return process((msg) => op.Apply(msg)); + } + else + { + try + { + using (var file = new StreamWriter(filename)) + { + return process(file.WriteLine); + } + } + catch (Exception e) + { + logMessage.Apply($"[warning] Unable to write state to '{filename}' ({e.Message})"); + return QVoid.Instance; + } + } + } + + + // Called by the diagnostics functions + public class SparseSimDumpMachine : Quantum.Diagnostics.DumpMachine + { + private SparseSimulator Simulator { get; } + + public SparseSimDumpMachine(SparseSimulator m) : base(m) + { + this.Simulator = m; + } + + public override Func __Body__ => (location) => + { + if (location == null) { throw new ArgumentNullException(nameof(location)); } + + return Simulator.Dump(location); + }; + } + + // Called by the diagnostics functions + public class QSimDumpRegister : Quantum.Diagnostics.DumpRegister + { + private SparseSimulator Simulator { get; } + + + public QSimDumpRegister(SparseSimulator m) : base(m) + { + this.Simulator = m; + } + + public override Func<(T, IQArray), QVoid> __Body__ => (__in) => + { + var (location, qubits) = __in; + + if (location == null) { throw new ArgumentNullException(nameof(location)); } + + return Simulator.Dump(location, qubits); + }; + } + + private string FormatMagnitude(double magnitude, double phase) => + (new String('*', (int)System.Math.Ceiling(20.0 * magnitude))).PadRight(20) + $" [ {magnitude:F6} ]"; + + /// + /// Returns a string that represents the phase of the amplitude. + /// + private string FormatAngle(double magnitude, double angle) + { + var PI = System.Math.PI; + var offset = PI / 16.0; + if (magnitude == 0.0) + { + return " "; + } + + var chart = " ---"; + if (angle > 0) + { + if (angle >= (0 * PI / 8) + offset && angle < ((1 * PI / 8) + offset)) { chart = " /-"; } + if (angle >= (1 * PI / 8) + offset && angle < ((2 * PI / 8) + offset)) { chart = " / "; } + if (angle >= (2 * PI / 8) + offset && angle < ((3 * PI / 8) + offset)) { chart = " +/ "; } + if (angle >= (3 * PI / 8) + offset && angle < ((4 * PI / 8) + offset)) { chart = " ↑ "; } + if (angle >= (4 * PI / 8) + offset && angle < ((5 * PI / 8) + offset)) { chart = " \\- "; } + if (angle >= (5 * PI / 8) + offset && angle < ((6 * PI / 8) + offset)) { chart = " \\ "; } + if (angle >= (6 * PI / 8) + offset && angle < ((7 * PI / 8) + offset)) { chart = "+\\ "; } + if (angle >= (7 * PI / 8) + offset) { chart = "--- "; } + } + else if (angle < 0) + { + var abs_angle = System.Math.Abs(angle); + if (abs_angle >= (0 * PI / 8) + offset && abs_angle < ((1 * PI / 8) + offset)) { chart = " \\+"; } + if (abs_angle >= (1 * PI / 8) + offset && abs_angle < ((2 * PI / 8) + offset)) { chart = " \\ "; } + if (abs_angle >= (2 * PI / 8) + offset && abs_angle < ((3 * PI / 8) + offset)) { chart = " -\\ "; } + if (abs_angle >= (3 * PI / 8) + offset && abs_angle < ((4 * PI / 8) + offset)) { chart = " ↓ "; } + if (abs_angle >= (4 * PI / 8) + offset && abs_angle < ((5 * PI / 8) + offset)) { chart = " /+ "; } + if (abs_angle >= (5 * PI / 8) + offset && abs_angle < ((6 * PI / 8) + offset)) { chart = " / "; } + if (abs_angle >= (6 * PI / 8) + offset && abs_angle < ((7 * PI / 8) + offset)) { chart = "-/ "; } + } + + return $" {chart} [ {angle,8:F5} rad ]"; + } + + /// + /// Returns a string for the amplitude's polar representation (magnitude/angle). + /// + private string FormatPolar(double magnitude, double angle) => + $"{FormatMagnitude(magnitude, angle)}{FormatAngle(magnitude, angle)}"; + + /// + /// Returns a string for the amplitude's cartesian representation (real + imagnary). + /// + private string FormatCartesian(double real, double img) => + $"{real,9:F6} + {img,9:F6} i"; + + /// + /// The method to use to format the amplitude into a string. + /// + private string Format(string label, double real, double img) + { + var amplitude = (real * real) + (img * img); + var angle = System.Math.Atan2(img, real); + + return $"{label}:\t" + + $"{FormatCartesian(real, img)}\t == \t" + + $"{FormatPolar(amplitude, angle)}"; + + } + + } + + + /// + /// Samples a random label from the wavefunction, proportional to its amplitude squared, + /// and returns it as a boolean array + /// + partial class Sample + { + private SparseSimulator simulator; + public class Native : Sample + { + public Native(IOperationFactory m) : base(m) + { + simulator = m as SparseSimulator; + } + public override Func, IQArray> __Body__ => (__in__) => + { + if (simulator != null) + { + return new QArray(((SparseSimulatorProcessor)simulator.QuantumProcessor).Sample(__in__)); + } + return new QArray( new bool[0]); + }; + } + } + + public partial class SparseSimulatorProcessor : QuantumProcessorBase { + + + // These two delegates exist so that C# can use StringCallback (which has a nice + // string object) and C++ can use DumpCallback (which can handle char*) + public delegate void StringCallback(string label, double real, double img); + + private delegate void DumpCallback(StringBuilder label, double real, double img); + [DllImport(simulator_dll)] + private static extern void Dump_cpp(uint sim, uint max_qubits, DumpCallback callback); + + // This gets called by the base class + public bool Dump(StringCallback callback, int max_id, IQArray qubits = null) + { + DumpCallback _callback = (StringBuilder label_builder, double real, double imag) => + { + callback(label_builder.ToString(), real, imag); + }; + + // It's expensive for the C++ simulator to separate two parts of a wavefunction, + // so if we just want the full state, it calls a different function + if (qubits == null) + { + Dump_cpp(Id, (uint)max_id, _callback); + return true; + } else + { + return DumpQubits_cpp(Id, qubits.Count(), qubits.Select(x => x.Id).ToArray(), _callback) ; + } + } + + [DllImport(simulator_dll)] + [return: MarshalAs(UnmanagedType.I1)] // necessary because C++ and C# represent bools differently + private static extern bool DumpQubits_cpp(uint sim, int length, int[] qubit_ids, DumpCallback callback); + + [DllImport(simulator_dll)] + private static extern uint num_qubits_cpp(uint sim); + + // Returns the amplitude of a specific state, given by a string representing its label + [DllImport(simulator_dll)] + private static extern void GetAmplitude_cpp(uint sim, uint label_length, char[] label, ref double real, ref double imag); + + public Microsoft.Quantum.Math.Complex GetAmplitude(string label) + { + double real = 0; + double imag = 0; + GetAmplitude_cpp(Id, (uint)label.Length, label.ToCharArray(), ref real, ref imag); + return new Microsoft.Quantum.Math.Complex((real, imag)); + } + + [DllImport(simulator_dll)] + private static extern IntPtr Sample_cpp(uint sim); + + public bool[] Sample(IQArray register) + { + var result = Marshal.PtrToStringAnsi(Sample_cpp(Id)).ToCharArray(); + Array.Reverse(result); + var bool_result = new bool[register.Length]; + for (int i=0; i < register.Length; i++) + { + bool_result[i] = (result[register[i].Id] == '1' ? true : false); + } + return bool_result; + } + } + + + + +} diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs new file mode 100644 index 00000000000..7d1f93e1139 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Probes.qs @@ -0,0 +1,57 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +namespace Microsoft.Quantum.SparseSimulation { + + open Microsoft.Quantum.Canon; + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Math; + open Microsoft.Quantum.Diagnostics; + + // These internal operations are physically unrealistic diagnostic + // tools that can read amplitudes and probabilities + // from a quantum state directly without modifyng it. + + // Returns an output with the same distribution + // as measuring all the qubits, but does not + // modify the quantum state. + internal function Sample(register : Qubit[]) : Bool[] { + body intrinsic; + } + + // Returns the amplitude of the quantum state which has a specific qubit + // label (e.g., "|101>" has label 5) + internal function GetAmplitudeFromInt(qubit : Qubit[], label : Int) : Complex { + body intrinsic; + } + internal function GetAmplitude(qubit : Qubit[], label : BigInt) : Complex { + body intrinsic; + } + + // Asserts (within some tolerance) that the probability of measuring the result passed as + // stateIndex is a specific value. + operation AssertProbBigInt (stateIndex : BigInt, expected : Double, qubits : Microsoft.Quantum.Arithmetic.LittleEndian, tolerance : Double) : Unit { + Fact(AbsD(AbsSquaredComplex (GetAmplitude(qubits!, stateIndex)) - expected) < tolerance, "Probability failed"); + } + + // Asserts that the amplitudes of an array of states in superposition on a collection of qubits are equal + // (within tolerance) to a specified array, up to a global phase. + operation AssertAmplitudes(stateIndices : BigInt[], expected : Complex[], qubits : Microsoft.Quantum.Arithmetic.LittleEndian, tolerance : Double) : Unit { + let basicC = GetAmplitude(qubits!, stateIndices[0]); + let invC = Complex(basicC::Real/AbsSquaredComplex(basicC), basicC::Imag/AbsSquaredComplex(basicC)); + let expInv = Complex(expected[0]::Real/AbsSquaredComplex(expected[0]), expected[0]::Imag/AbsSquaredComplex(expected[0])); + for idx in 0..Length(stateIndices)-1 { + let c = TimesC(GetAmplitude(qubits!, stateIndices[idx]), invC); + let expC = TimesC(expected[idx], expInv); + Fact(AbsD(c::Real - expC::Real) + AbsD(c::Imag - expC::Imag) < 0.001, "Amplitude incorrect"); + } + } + + // Asserts that if the qubits in `register` were measured, they would pass some check. + // The measurement would product an array of boolean values, so the argument `assertion` + // should check whether an input array of boolean values, representing a measurement + // output, satisfies some criteria. + operation AssertSample(assertion : (Bool[] -> Bool), register : Qubit[]) : Unit { + Fact(assertion(Sample(register)), "Sampling failed"); + } +} diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs new file mode 100644 index 00000000000..93e2345ff88 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/Simulator.ApplyAnd.cs @@ -0,0 +1,60 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +// Overrides the basic AND gate for faster execution + +using Microsoft.Quantum.Simulation.Core; +using System; +using System.Linq; +using Microsoft.Quantum.Canon; + +namespace Microsoft.Quantum.SparseSimulation +{ + public class ApplyAndWrapper : ApplyAnd + { + private SparseSimulator simulator; + public ApplyAndWrapper(IOperationFactory m) : base(m) + { + simulator = m as SparseSimulator; + } + + public override Func<(Qubit, Qubit, Qubit), QVoid> __Body__ => simulator == null ? base.__Body__ : (args) => + { + simulator.ApplyAnd(new QArray(args.Item1, args.Item2), args.Item3); + return QVoid.Instance; + }; + + public override Func<(Qubit, Qubit, Qubit), QVoid> __AdjointBody__ => simulator == null ? base.__AdjointBody__ : (args) => + { + simulator.AdjointApplyAnd(new QArray(args.Item1, args.Item2), args.Item3); + return QVoid.Instance; + }; + + public override Func<(IQArray, (Qubit, Qubit, Qubit)), QVoid> __ControlledBody__ => simulator == null ? base.__ControlledBody__ : (args) => + { + simulator.ApplyAnd(new QArray(args.Item1.Concat(new QArray(args.Item2.Item1, args.Item2.Item2))), args.Item2.Item3); + return QVoid.Instance; + }; + + public override Func<(IQArray, (Qubit, Qubit, Qubit)), QVoid> __ControlledAdjointBody__ => simulator == null ? base.__ControlledAdjointBody__ : (args) => + { + simulator.AdjointApplyAnd(new QArray(args.Item1.Concat(new QArray(args.Item2.Item1, args.Item2.Item2))), args.Item2.Item3); + return QVoid.Instance; + }; + } + + public partial class SparseSimulator + { + // Wrappers for the relevant functions of SparseSimulatorProcessor + // Since And/AdjointAnd are inherently multi-controlled, there is only a multi-controlled + // emulator + public void ApplyAnd(IQArray controls, Qubit target) + { + ((SparseSimulatorProcessor)this.QuantumProcessor).MCApplyAnd(controls, target); ; + } + public void AdjointApplyAnd(IQArray controls, Qubit target) + { + ((SparseSimulatorProcessor)this.QuantumProcessor).MCAdjointApplyAnd(controls, target); + } + } +} \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs new file mode 100644 index 00000000000..7b256c8cbcc --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.cs @@ -0,0 +1,374 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +using System; +using Microsoft.Quantum.Simulation.QuantumProcessor; +using Microsoft.Quantum.Simulation.Common; +using Microsoft.Quantum.Simulation.Core; +using System.Runtime.InteropServices; +using Microsoft.Quantum.Simulation.Simulators; +using System.Linq; +using Microsoft.Quantum.Canon; + +using Microsoft.Quantum.SparseSimulation; + +namespace Microsoft.Quantum.SparseSimulation +{ + + // Class that can be used as a simulator, but is mostly a wrapper + // for the underlying SparseSimulatorProcessor + public partial class SparseSimulator : QuantumProcessorDispatcher, IDisposable + { + public uint Id { get; } + + // If set to true, the simulator will output the name of all operations + // it executes to the console, as it executes them + public void SetLogging(bool isLogging) + { + ((SparseSimulatorProcessor)this.QuantumProcessor).isLogging = isLogging; + } + + public SparseSimulator() : base(new SparseSimulatorProcessor() ) + { + // Emulates AND + Register(typeof(ApplyAnd), typeof(ApplyAndWrapper), typeof(IUnitary)); + Id = ((SparseSimulatorProcessor)this.QuantumProcessor).Id; + } + + public void Dispose() + { + Dispose(true); + } + + // Clears memory, specifically by telling the C++ to clear the memory + public void Dispose(bool Disposing) + { + ((SparseSimulatorProcessor)this.QuantumProcessor).Dispose(); + } + + // Sets the random seed in the C++ code, used for the randomness + // of measurements. + // Default value is the default for std::mt19937 + public void SetSeed(uint newSeed = 5489) + { + ((SparseSimulatorProcessor)this.QuantumProcessor).SetSeed(newSeed); + } + + } + + + // Class that performs all gate operations + // Mostly a wrapper for C++ code + partial class SparseSimulatorProcessor : QuantumProcessorBase + { + // References the C++ dll that is copied in during the build + private const string simulator_dll = "SparseQuantumSimulator.dll"; + // Controls whether it outputs all operations for debugging + public bool isLogging = false; + // For debugging + private string stackDepth; + + // C++ code uses the Id to refer to this simulator in a vector of simulators + public uint Id { get; } + + public string Name => "Sparse Simulator"; + + [DllImport(simulator_dll)] + private static extern uint init_cpp(uint num_qubits); + public SparseSimulatorProcessor(uint num_qubits = 64) + { + Id = init_cpp(num_qubits); + stackDepth = ""; + } + + + [DllImport(simulator_dll)] + private static extern void seed_cpp(uint sim, uint new_seed); + public void SetSeed(uint newSeed = 5489) + { + seed_cpp(Id, newSeed); + } + + + + // Basic gates + [DllImport(simulator_dll)] + private static extern void MCApplyAnd_cpp(uint sim, int length, int[] controls, int target); + public void MCApplyAnd(IQArray controls, Qubit qubit) + { + MCApplyAnd_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + + [DllImport(simulator_dll)] + private static extern void MCAdjointApplyAnd_cpp(uint sim, int length, int[] controls, int target); + public void MCAdjointApplyAnd(IQArray controls, Qubit qubit) + { + MCAdjointApplyAnd_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + + + [DllImport(simulator_dll)] + private static extern void X_cpp(uint sim, int qubit_id); + public override void X(Qubit qubit) + { X_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void MCX_cpp(uint sim, int length, int[] controls, int target); + public override void ControlledX(IQArray controls, Qubit qubit) + { + MCX_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void Z_cpp(uint sim, int qubit_id); + public override void Z(Qubit qubit) + { Z_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void MCZ_cpp(uint sim, int length, int[] controls, int target); + public override void ControlledZ(IQArray controls, Qubit qubit) + { + MCZ_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void Y_cpp(uint sim, int qubit_id); + public override void Y(Qubit qubit) + { Y_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void MCY_cpp(uint sim, int length, int[] controls, int target); + public override void ControlledY(IQArray controls, Qubit qubit) + { + MCY_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern int H_cpp(uint sim, int qubit_id); + public override void H(Qubit qubit) + { H_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void MCH_cpp(uint sim, int length, int[] controls, int target); + public override void ControlledH(IQArray controls, Qubit qubit) + { + MCH_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void S_cpp(uint sim, int qubit_id); + public override void S(Qubit qubit) + { S_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void AdjS_cpp(uint sim, int qubit_id); + public override void SAdjoint(Qubit qubit) + { AdjS_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void T_cpp(uint sim, int qubit_id); + public override void T(Qubit qubit) + { T_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + private static extern void AdjT_cpp(uint sim, int qubit_id); + public override void TAdjoint(Qubit qubit) + { AdjT_cpp(Id, qubit.Id); } + [DllImport(simulator_dll)] + [return: MarshalAs(UnmanagedType.I1)] // necessary because C++ and C# represent bools differently + private static extern bool M_cpp(uint sim, int qubit_id); + public override Result M(Qubit qubit) + { + if (M_cpp(Id, qubit.Id)) + { + return Result.One; + } + else + { + return Result.Zero; + } + } + [DllImport(simulator_dll)] + [return: MarshalAs(UnmanagedType.I1)] // necessary because C++ and C# represent bools differently + private static extern bool Measure_cpp(uint sim, int length, int[] basis, int[] qubits); + public override Result Measure(IQArray bases, IQArray qubits) + { + if (Measure_cpp(Id, bases.Count(), bases.Select(x => (int)x).ToArray(), qubits.Select(x => x.Id).ToArray())) { + return Result.One; + } else { return Result.Zero; } + } + + [DllImport(simulator_dll)] + private static extern void Reset_cpp(uint sim, int qubit_id); + public override void Reset(Qubit qubit) + { + Reset_cpp(Id, qubit.Id); + } + + + [DllImport(simulator_dll)] + private static extern void R_cpp(uint sim, int axis, double theta, int qubit_id); + public override void R(Pauli axis, double theta, Qubit qubit) + { + R_cpp(Id, (int)axis, theta, qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void MCR_cpp(uint sim, int basis, double angle, int length, int[] controls, int target); + public override void ControlledR(IQArray controls, Pauli axis, double theta, Qubit qubit) + { + MCR_cpp(Id, (int)axis, theta, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void Rfrac_cpp(uint sim, int axis, long numerator, long power, int qubit_id); + public override void RFrac(Pauli axis, long numerator, long power, Qubit qubit) + { + Rfrac_cpp(Id, (int)axis, numerator, power, qubit.Id); + } + + [DllImport(simulator_dll)] + private static extern void MCRFrac_cpp(uint sim, int basis, long numerator, long power, int length, int[] controls, int target); + public override void ControlledRFrac(IQArray controls, Pauli axis, long numerator, long power, Qubit qubit) + { + MCRFrac_cpp(Id, (int)axis, numerator, power, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void R1_cpp(uint sim, double theta, int qubit_id); + public override void R1(double theta, Qubit qubit) + { + R1_cpp(Id, theta, qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void MCR1_cpp(uint sim, double angle, int length, int[] controls, int target); + public override void ControlledR1(IQArray controls, double theta, Qubit qubit) + { + MCR1_cpp(Id, theta, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void R1frac_cpp(uint sim, long numerator, long power, int qubit_id); + public override void R1Frac(long numerator, long power, Qubit qubit) + { + R1frac_cpp(Id, numerator, power, qubit.Id); + } + [DllImport(simulator_dll)] + private static extern void MCR1Frac_cpp(uint sim, long numerator, long power, int length, int[] controls, int target); + public override void ControlledR1Frac(IQArray controls, long numerator, long power, Qubit qubit) + { + MCR1Frac_cpp(Id, numerator, power, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit.Id); + } + + [DllImport(simulator_dll)] + private static extern void SWAP_cpp(uint sim, int qubit_id_1, int qubit_id_2); + public override void SWAP(Qubit qubit1, Qubit qubit2) + { + SWAP_cpp(Id, qubit1.Id, qubit2.Id); + } + [DllImport(simulator_dll)] + private static extern void MCSWAP_cpp(uint sim, int length, int[] controls, int qubit_id_1, int qubit_id_2); + public override void ControlledSWAP(IQArray controls, Qubit qubit1, Qubit qubit2) + { + MCSWAP_cpp(Id, controls.Count(), controls.Select(x => x.Id).ToArray(), qubit1.Id, qubit2.Id); + } + + [DllImport(simulator_dll)] + private static extern void Exp_cpp(uint sim, int length, int[] b, double phi, int[] q); + public override void Exp(IQArray paulis, double theta, IQArray qubits) + { + Exp_cpp(Id, paulis.Count(), paulis.Select(x => (int)x).ToArray(), theta, qubits.Select(x => x.Id).ToArray()); + } + [DllImport(simulator_dll)] + private static extern void MCExp_cpp(uint sim, int controls_length, int length, int[] c, int[] b, double phi, int[] q); + public override void ControlledExp(IQArray controls, IQArray paulis, double theta, IQArray qubits) + { + MCExp_cpp(Id, controls.Count(), paulis.Count(), controls.Select(x => x.Id).ToArray(), paulis.Select(x => (int)x).ToArray(), theta, qubits.Select(x => x.Id).ToArray()); + } + + [DllImport(simulator_dll)] + [return: MarshalAs(UnmanagedType.I1)] // necessary because C++ and C# represent bools differently + private static extern bool Assert_cpp(uint sim, int length, int[] b, int[] q, bool result); + public override void Assert(IQArray bases, IQArray qubits, Result result, string msg) + { + // Relies on C++ catching any assertion errors and returning false + if (!Assert_cpp( + Id, + bases.Count(), + bases.Select(x => (int)x).ToArray(), + qubits.Select(x => x.Id).ToArray(), + result.Equals(Result.One))) + { + throw new Exception(msg); + } + } + + [DllImport(simulator_dll)] + private static extern double JointEnsembleProbability_cpp(uint sim, int length, int[] basis, int[] qubits); + public override void AssertProb(IQArray bases, IQArray qubits, double probabilityOfZero, string msg, double tol) + { + double result = JointEnsembleProbability_cpp(Id, bases.Count(), bases.Select(x => (int)x).ToArray(), qubits.Select(x => x.Id).ToArray()); + if (System.Math.Abs(1.0 - result - probabilityOfZero) > tol) + { + var extendedMsg = $"{msg}\n\tExpected:\t{probabilityOfZero}\n\tActual:\t{result}"; + IgnorableAssert.Assert(false, extendedMsg); + throw new ExecutionFailException(extendedMsg); + } + } + + + [DllImport(simulator_dll)] + private static extern void allocateQubit_cpp(uint sim, int qubit_id); + public override void OnAllocateQubits(IQArray qubits) + { + try + { + foreach (Qubit qubit in qubits) + { + allocateQubit_cpp(Id, qubit.Id); + } + } catch (ExternalException ex) + { + // Possibly C++ errors: qubit already occupied, or not enough qubits + throw ex; + } + } + + [DllImport(simulator_dll)] + private static extern void releaseQubit_cpp(uint sim, int qubit_id); + public override void OnReleaseQubits(IQArray qubits) + { + try + { + foreach (Qubit qubit in qubits) + { + releaseQubit_cpp(Id, qubit.Id); + } + } catch (ExternalException) + { + // Possible C++ errors: Qubit not in zero state + throw new Microsoft.Quantum.Simulation.Simulators.Exceptions.ReleasedQubitsAreNotInZeroState(); + } + } + + // Tells C++ to delete any objects, to manage memory + public void Dispose() + { + Dispose(true); + } + + [DllImport(simulator_dll)] + private static extern void destroy_cpp(uint sim); + public void Dispose(bool Disposing) + { + destroy_cpp(Id); + } + + public override void OnOperationStart(ICallable operation, IApplyData arguments) + { + // Writes the operation to console if it is logging + if (isLogging) + { + Console.WriteLine(stackDepth + operation.FullName + "{"); + stackDepth += " "; + } + } + + public override void OnOperationEnd(ICallable operation, IApplyData arguments) + { + + if (isLogging) + { + stackDepth = stackDepth[0..(stackDepth.Length - 1)]; + Console.WriteLine(stackDepth + "}"); + } + } + + + } +} diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.csproj b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.csproj new file mode 100644 index 00000000000..ce7ad730387 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorCS/SparseSimulator.csproj @@ -0,0 +1,28 @@ + + + + Library + netcoreapp3.1 + true + + + + + ..\Native\build\libSparseQuantumSimulator.dylib + ..\Native\build\libSparseQuantumSimulator.so + ..\Native\build\$(Configuration)\SparseQuantumSimulator.dll + $(QsimDllMac) + $(QsimDllLinux) + $(QsimDllWindows) + + + + + SparseQuantumSimulator.dll + PreserveNewest + true + + + + + diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CMakeLists.txt b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CMakeLists.txt new file mode 100644 index 00000000000..a208bb9df56 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 3.10) +project(SparseQuantumSimulator) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +include(CTest) +enable_testing() + +set(CMAKE_MACOSX_RPATH 1) +add_executable(SparseSimulatorTests SparseSimulatorTests.cpp TestHelpers.cpp) +add_test(SparseSimulatorTests SparseSimulatorTests) \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp new file mode 100644 index 00000000000..409e747a77e --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/CSharpIntegrationTests.cpp @@ -0,0 +1,307 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#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" +#include +#include + +using namespace Microsoft::VisualStudio::CppUnitTestFramework; +using namespace Microsoft::Quantum::SPARSESIMULATOR; +using namespace SparseSimulatorTestHelpers; + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace CSharpIntegrationTests +{ + TEST_CLASS(CSharpIntegrationTests) + { + template + void MultiExpReferenceTest( + std::function qubit_prep, + std::function qubit_clear + ) { + const qubit_label_type zero(0); + logical_qubit_id* qubits = new logical_qubit_id[3]; + qubits[0] = 0; + qubits[1] = 1; + qubits[2] = 2; + int* Paulis = new int[3]; + for (int intPaulis = 0; intPaulis < 4 * 4 * 4; intPaulis++) { + Paulis[0] = intPaulis % 4; + Paulis[1] = (intPaulis / 4 ) % 4; + Paulis[2] = intPaulis / 16; + + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + unsigned sim = init_cpp(32); + qubit_prep(sim); + + std::vector vector_rep(8, 0.0); + for (unsigned i = 0; i < 8; i++) { + vector_rep[i] = getSimulator(sim)->probe(std::bitset<3>(i).to_string()); + } + // New simulator Exp + Exp_cpp(sim, 3, Paulis, angle, qubits); + // Old simulator Exp + std::vector actualPaulis = { (Gates::Basis)Paulis[0],(Gates::Basis)Paulis[1], (Gates::Basis)Paulis[2] }; + apply_exp(vector_rep, actualPaulis, angle, std::vector{ 0, 1, 2 }); + for (unsigned i = 0; i < 8; i++) { + amplitude result = getSimulator(sim)->probe(std::bitset<3>(i).to_string()); + assert_amplitude_equality(vector_rep[i], result); + } + Exp_cpp(sim, 3, Paulis, -angle, qubits); + qubit_clear(sim); + } + } + } + public: + TEST_METHOD(initializationTest) { + unsigned sim = init_cpp(32); + } + + TEST_METHOD(AllocationTest) { + unsigned sim = init_cpp(32); + allocateQubit_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + } + TEST_METHOD(AllocateRebuildTest) { + unsigned sim = init_cpp(64); + for (int i = 0; i < 1024; i++) { + allocateQubit_cpp(sim, i); + getSimulator(sim)->X(i); + getSimulator(sim)->update_state(); + } + for (int i = 0; i < 1024; i++) { + getSimulator(sim)->X(i); + releaseQubit_cpp(sim, i); + } + } + + TEST_METHOD(XTest) { + unsigned sim = init_cpp(32); + allocateQubit_cpp(sim, 0); + X_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0, 0.0); + X_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + } + TEST_METHOD(ZTest) { + unsigned sim = init_cpp(32); + allocateQubit_cpp(sim, 0); + Z_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 1.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 0.0, 0.0); + Z_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 1.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 0.0, 0.0); + X_cpp(sim, 0); + Z_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), -1.0, 0.0); + Z_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0, 0.0); + X_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + } + TEST_METHOD(HTest) { + unsigned sim = init_cpp(32); + allocateQubit_cpp(sim, 0); + H_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 1.0 / sqrt(2.0), 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0 / sqrt(2.0), 0.0); + H_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 1.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 0.0, 0.0); + X_cpp(sim, 0); + H_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 1.0 / sqrt(2.0), 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), -1.0 / sqrt(2.0), 0.0); + H_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0, 0.0); + X_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + } + + TEST_METHOD(TGateTest) { + unsigned sim = init_cpp(32); + allocateQubit_cpp(sim, 0); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 1.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 0.0, 0.0); + X_cpp(sim, 0); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0 / sqrt(2.0), 1.0 / sqrt(2.0)); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 0.0, 1.0); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), -1.0 / sqrt(2.0), 1.0 / sqrt(2.0)); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), -1.0, 0.0); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), -1.0 / sqrt(2.0), -1.0 / sqrt(2.0)); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 0.0, -1.0); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0 / sqrt(2.0), -1.0 / sqrt(2.0)); + T_cpp(sim, 0); + assert_amplitude_equality(getSimulator(sim)->probe("0"), 0.0, 0.0); + assert_amplitude_equality(getSimulator(sim)->probe("1"), 1.0, 0.0); + X_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + } + + TEST_METHOD(HCancellationTest) + { + int n_qubits = 16; + unsigned sim = init_cpp(n_qubits); + size_t buckets = 0; + for (int i = 0; i < n_qubits; i++) { + allocateQubit_cpp(sim, i); + H_cpp(sim, i); + } + for (int i = n_qubits - 1; i >= 0; i--) { + H_cpp(sim, i); + // If the H do not cancel out, release will fail in an opaque way + releaseQubit_cpp(sim, i); + } + } + + TEST_METHOD(HXZCommutationTest) + { + const int n_qubits = 16; + unsigned sim = init_cpp(n_qubits); + for (int i = 0; i < n_qubits; i++) { + allocateQubit_cpp(sim, i); + H_cpp(sim, i); + } + std::bitset one_state = 0; + for (int i = 0; i < n_qubits - 1; i += 2) { + Z_cpp(sim, i); + X_cpp(sim, i+1); + one_state.set(i); + } + for (int i = n_qubits - 1; i >= 0; i--) { + H_cpp(sim, i); + } + for (__int64 i = 0; i< pow(2, n_qubits); i++) { + amplitude state = getSimulator(sim)->probe(std::bitset(i).to_string()); + if (i == one_state.to_ulong()) { + assert_amplitude_equality(state, 1.0, 0.0); + } + else { + assert_amplitude_equality(state, 0.0, 0.0); + } + } + } + + TEST_METHOD(ResetTest) + { + const int n_qubits = 16; + unsigned sim = init_cpp(n_qubits); + allocateQubit_cpp(sim, 0); + Reset_cpp(sim, 0); + amplitude state = getSimulator(sim)->probe("0"); + assert_amplitude_equality(state, 1.0, 0.0); + X_cpp(sim, 0); + Reset_cpp(sim, 0); + state = getSimulator(sim)->probe("0"); + // No qubit exists; should have amplitude 0 + assert_amplitude_equality(state, 1.0, 0.0); + allocateQubit_cpp(sim, 1); + X_cpp(sim, 0); + logical_qubit_id* controls = new logical_qubit_id{ 0 }; + MCX_cpp(sim, 1, controls, 1); + Reset_cpp(sim, 0); + state = getSimulator(sim)->probe("00"); + assert_amplitude_equality(state, 0.0, 0.0); + state = getSimulator(sim)->probe("10"); + assert_amplitude_equality(state, 1.0, 0.0); + Reset_cpp(sim, 1); + state = getSimulator(sim)->probe("00"); + assert_amplitude_equality(state, 1.0, 0.0); + state = getSimulator(sim)->probe("10"); + assert_amplitude_equality(state, 0.0, 0.0); + releaseQubit_cpp(sim, 1); + releaseQubit_cpp(sim, 0); + } + + TEST_METHOD(MultiExpWithHTest) { + const int num_qubits = 32; + auto qubit_prep = [](unsigned sim ) { + H_cpp(sim, 0); + H_cpp(sim, 1); + H_cpp(sim, 2); + }; + auto qubit_clear = [](unsigned sim) { + H_cpp(sim, 2); + releaseQubit_cpp(sim, 2); + H_cpp(sim, 1); + releaseQubit_cpp(sim, 1); + H_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + }; + MultiExpReferenceTest(qubit_prep, qubit_clear); + } + + TEST_METHOD(MultiExpBasisTest) { + const int num_qubits = 32; + auto qubit_prep = [](unsigned sim, int index) { + if ((index & 1) == 0) { X_cpp(sim, 0); } + if ((index & 2) == 0) { X_cpp(sim, 1); } + if ((index & 4) == 0) { X_cpp(sim, 2); } + }; + auto qubit_clear = [](unsigned sim, int index) { + if ((index & 1) == 0) { X_cpp(sim, 0); } + releaseQubit_cpp(sim, 0); + if ((index & 2) == 0) { X_cpp(sim, 1); } + releaseQubit_cpp(sim, 1); + if ((index & 4) == 0) { X_cpp(sim, 2); } + releaseQubit_cpp(sim, 2); + }; + for (int i = 0; i < 8; i++) { + MultiExpReferenceTest([=](unsigned sim) {qubit_prep(sim, i); }, [=](unsigned sim) {qubit_clear(sim, i); }); + } + } + + TEST_METHOD(R1Test) { + const int num_qubits = 32; + amplitude result0; + amplitude result1; + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + unsigned sim = init_cpp(num_qubits); + H_cpp(sim, 0); + R1_cpp(sim, angle, 0); + result0 = getSimulator(sim)->probe("0"); + result1 = getSimulator(sim)->probe("1"); + assert_amplitude_equality(result0, 1.0 / sqrt(2.0)); + assert_amplitude_equality(result1, amplitude(cos(angle), sin(angle))/sqrt(2.0)); + R1_cpp(sim, -angle, 0); + result0 = getSimulator(sim)->probe("0"); + result1 = getSimulator(sim)->probe("1"); + assert_amplitude_equality(result0, 1.0 / sqrt(2.0)); + assert_amplitude_equality(result1, 1.0 / sqrt(2.0)); + H_cpp(sim, 0); + releaseQubit_cpp(sim, 0); + destroy_cpp(sim); + } + } + }; +} \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp new file mode 100644 index 00000000000..d3674644420 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.cpp @@ -0,0 +1,772 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#include "pch.h" +#include "CppUnitTest.h" +#include "../Native/SparseSimulator.h" +#include "TestHelpers.hpp" +#include +#include + +using namespace Microsoft::VisualStudio::CppUnitTestFramework; +using namespace Microsoft::Quantum::SPARSESIMULATOR; +using namespace SparseSimulatorTestHelpers; + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + + + + +namespace SparseSimulatorTests +{ + TEST_CLASS(SparseSimulatorTests) + { + + template + void MultiExpTest( + std::function qubit_prep, + std::function qubit_clear + ) { + for (int intPaulis = 0; intPaulis < 4 * 4 * 4; intPaulis++) { + std::vector Paulis{ + (Gates::Basis)(intPaulis % 4), + (Gates::Basis)((intPaulis / 4) % 4), + (Gates::Basis)(intPaulis / 16) + }; + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + SparseSimulator sim = SparseSimulator(num_qubits); + std::vector qubits{ 0,1,2 }; + qubit_prep(sim); + std::vector vector_rep(8, 0.0); + for (unsigned i = 0; i < 8; i++) { + vector_rep[i] = sim.probe(std::bitset<3>(i).to_string()); + } + // New simulator Exp + sim.Exp(Paulis, angle, qubits); + // Old simulator Exp + apply_exp(vector_rep, Paulis, angle, std::vector{ 0, 1, 2 }); + for (unsigned i = 0; i < 8; i++) { + amplitude result = sim.probe(std::bitset<3>(i).to_string()); + assert_amplitude_equality(vector_rep[i], result); + } + sim.Exp(Paulis, -angle, qubits); + qubit_clear(sim); + } + } + } + + public: + // Tests comparisons of bitstrings + TEST_METHOD(LabelComparisonTest) { + const logical_qubit_id num_qubits = 1024; + SparseSimulator sim = SparseSimulator(num_qubits); + uint64_t i = 0; + uint64_t j; + uint64_t k = 0; + qubit_label_type label1(0); + qubit_label_type label2(1); + + for (i = 0; i < 500; i++){ + k += i * i * i * i; + uint64_t m = 0; + label1 = qubit_label_type(k); + for (j = 0; j < 500; j++){ + m += j * j * j * j; + label2 = qubit_label_type(m); + wchar_t message[100]; + swprintf(message, sizeof(message)/sizeof(*message), L"Comparing %llu to %llu\n", k, m); + Assert::AreEqual(k < m, label1 < label2, message); + } + } + } + // Tests that the X gate flips the computational basis states + TEST_METHOD(XGateTest) { + const logical_qubit_id num_qubits = 32; + SparseSimulator sim = SparseSimulator(num_qubits); + sim.X(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + sim.X(0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + } + + // Tests Z on computational basis states + TEST_METHOD(ZGateTest) { + const logical_qubit_id num_qubits = 32; + SparseSimulator sim = SparseSimulator(num_qubits); + sim.Z(0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.Z(0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.X(0); + sim.Z(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), -1.0, 0.0); + sim.Z(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + } + + // Tests H on computational basis states + TEST_METHOD(HGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + sim.H(0); + assert_amplitude_equality(sim.probe("0"), 1.0 / sqrt(2.0), 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0 / sqrt(2.0), 0.0); + sim.H(0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.X(0); + sim.H(0); + assert_amplitude_equality(sim.probe("0"), 1.0 / sqrt(2.0), 0.0); + assert_amplitude_equality(sim.probe("1"), -1.0 / sqrt(2.0), 0.0); + sim.H(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + } + + // Tests powers of T on computational basis states + TEST_METHOD(TGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.X(0); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0 / sqrt(2.0), 1.0 / sqrt(2.0)); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 1.0); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), -1.0 / sqrt(2.0), 1.0 / sqrt(2.0)); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), -1.0, 0.0); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), -1.0 / sqrt(2.0), -1.0 / sqrt(2.0)); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, -1.0); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0 / sqrt(2.0), -1.0 / sqrt(2.0)); + sim.T(0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + } + + // Tests Rx on computational basis states, for angles between 0 and pi/2 + TEST_METHOD(RxGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + SparseSimulator sim = SparseSimulator(num_qubits); + sim.R(Gates::Basis::PauliX, angle, 0); + assert_amplitude_equality(sim.probe("0"), cos(angle / 2.0), 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, -sin(angle / 2.0)); + sim.R(Gates::Basis::PauliX, -angle, 0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.X(0); + sim.R(Gates::Basis::PauliX, angle, 0); + assert_amplitude_equality(sim.probe("0"), 0.0, -sin(angle / 2.0)); + assert_amplitude_equality(sim.probe("1"), cos(angle / 2.0), 0.0); + sim.R(Gates::Basis::PauliX, -angle, 0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + } + } + + // Tests Ry on computational basis states, for angles between 0 and pi/2 + TEST_METHOD(RyGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + SparseSimulator sim = SparseSimulator(num_qubits); + sim.R(Gates::Basis::PauliY, angle, 0); + assert_amplitude_equality(sim.probe("0"), cos(angle / 2.0), 0.0); + assert_amplitude_equality(sim.probe("1"), sin(angle / 2.0), 0.0); + sim.R(Gates::Basis::PauliY, -angle, 0); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.X(0); + sim.R(Gates::Basis::PauliY, angle, 0); + assert_amplitude_equality(sim.probe("0"), -sin(angle / 2.0), 0.0); + assert_amplitude_equality(sim.probe("1"), cos(angle / 2.0), 0.0); + sim.R(Gates::Basis::PauliY, -angle, 0); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + } + } + + // Tests Rz on computational basis states, for angles between 0 and pi/2 + TEST_METHOD(RzGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + SparseSimulator sim = SparseSimulator(num_qubits); + logical_qubit_id qubit = 0; + sim.R(Gates::Basis::PauliZ, angle, qubit); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.R(Gates::Basis::PauliZ, -angle, qubit); + assert_amplitude_equality(sim.probe("0"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 0.0, 0.0); + sim.X(qubit); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + sim.R(Gates::Basis::PauliZ, angle, qubit); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), cos(angle), sin(angle)); + sim.R(Gates::Basis::PauliZ, -angle, qubit); + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + } + } + + // Tests CNOT on all 2-qubit computational basis stats + TEST_METHOD(CNOTGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + logical_qubit_id qubits[2]{ 0, 1 }; + sim.MCX({ qubits[0] }, qubits[1]); + assert_amplitude_equality(sim.probe("00"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 0.0, 0.0); + sim.X(qubits[0]); + sim.MCX({ qubits[0] }, qubits[1]); + assert_amplitude_equality(sim.probe("00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 1.0, 0.0); + sim.MCX({ qubits[0] }, qubits[1]); + assert_amplitude_equality(sim.probe("00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 0.0, 0.0); + sim.X(qubits[0]); + sim.X(qubits[1]); + sim.MCX({ qubits[0] }, qubits[1]); + assert_amplitude_equality(sim.probe("00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 0.0, 0.0); + } + + + + // Tests all possible computational basis states + // for some number of controls and one target + TEST_METHOD(MCXGateTest) { + const size_t n_qubits = 7; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(n_qubits); + std::vector qubits(n_qubits); + std::generate(qubits.begin(), qubits.end(), [] { static int i{ 0 }; return i++; }); + std::vector controls(qubits.begin() + 1, qubits.end()); + logical_qubit_id target = qubits[0]; + for (logical_qubit_id i = 0; i < pow(2, n_qubits - 1); i++) { // the bitstring of the controls + sim.MCX(controls, target); + for (logical_qubit_id j = 0; j < pow(2, n_qubits - 1); j++) { // bitstring to test + std::bitset j_bits = j; // a bitset for the string to test, with the target as 0 + j_bits = j_bits << 1; + std::bitset j_odd_bits = j_bits; // same as j, but with the target as 1 + j_odd_bits.set(0); + if (j != i) { // The controls are not in this state, so the amplitude should be 0 + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + else if (i == pow(2, n_qubits - 1) - 1) { // All controls are 1, so this should flip the output + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 1.0, 0.0); + } + else { // This is the state of the controls, but they are not all 1, so nothing should have happened + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + } + // Since MCX^2 = I, this should undo anything previously + sim.MCX(controls, target); + for (logical_qubit_id j = 0; j < pow(2, n_qubits - 1); j++) { + std::bitset j_bits = j; + j_bits = j_bits << 1; + std::bitset j_odd_bits = j_bits; + j_odd_bits.set(0); + if (j != i) { // The controls are not in this state, so the amplitude should be 0 + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + else { // This is the state of the controls, but the final qubit should be 0 + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + } + // Update the controls + std::bitset diff = i ^ (i + 1); + for (logical_qubit_id j = 0; j < n_qubits - 1; j++) { + if (diff[j]) sim.X(controls[j]); + } + } + } + + + // Tests a controlled Y + // Same logic as the MCXGateTest + TEST_METHOD(MCYGateTest) { + const size_t n_qubits = 7; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(n_qubits); + std::vector qubits(n_qubits); + std::generate(qubits.begin(), qubits.end(), [] { static int i{ 0 }; return i++; }); + std::vector controls(qubits.begin() + 1, qubits.end()); + logical_qubit_id target = qubits[0]; + for (logical_qubit_id i = 0; i < pow(2, n_qubits - 1); i++) { + sim.MCY(controls, target); + for (logical_qubit_id j = 0; j < pow(2, n_qubits - 1); j++) { + std::bitset j_bits = j; + j_bits = j_bits << 1; + std::bitset j_odd_bits = j_bits; + j_odd_bits.set(0); + if (j != i) { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + else if (i == pow(2, n_qubits - 1) - 1) { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 1.0); + } + else { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + } + sim.MCY(controls, target); + for (logical_qubit_id j = 0; j < pow(2, n_qubits - 1); j++) { + std::bitset j_bits = j; + j_bits = j_bits << 1; + std::bitset j_odd_bits = j_bits; + j_odd_bits.set(0); + if (j != i) { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + else { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + } + std::bitset diff = i ^ (i + 1); + for (logical_qubit_id j = 0; j < n_qubits - 1; j++) { + if (diff[j]) sim.X(controls[j]); + } + } + } + // Tests a controlled Z + // Same logic as the MCXGateTest + TEST_METHOD(MCZGateTest) { + const size_t n_qubits = 7; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(n_qubits); + std::vector qubits(n_qubits); + std::generate(qubits.begin(), qubits.end(), [] { static int i{ 0 }; return i++; }); + std::vector controls(qubits.begin() + 1, qubits.end()); + logical_qubit_id target = qubits[0]; + sim.H(target); + for (logical_qubit_id i = 0; i < pow(2, n_qubits - 1); i++) { + sim.MCZ(controls, target); + for (logical_qubit_id j = 0; j < pow(2, n_qubits - 1); j++) { + std::bitset j_bits = j; + j_bits = j_bits << 1; + std::bitset j_odd_bits = j_bits; + j_odd_bits.set(0); + if (j != i) { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + else if (i == pow(2, n_qubits - 1) - 1) { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0/sqrt(2.0), 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), -1.0 / sqrt(2.0), 0.0); + } + else { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0 / sqrt(2.0), 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 1.0 / sqrt(2.0), 0.0); + } + } + sim.MCZ(controls, target); + for (logical_qubit_id j = 0; j < pow(2, n_qubits - 1); j++) { + std::bitset j_bits = j; + j_bits = j_bits << 1; + std::bitset j_odd_bits = j_bits; + j_odd_bits.set(0); + if (j != i) { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 0.0, 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 0.0, 0.0); + } + else { + assert_amplitude_equality(sim.probe(j_bits.to_string()), 1.0 / sqrt(2.0), 0.0); + assert_amplitude_equality(sim.probe(j_odd_bits.to_string()), 1.0 / sqrt(2.0), 0.0); + } + } + std::bitset diff = i ^ (i + 1); + for (logical_qubit_id j = 0; j < n_qubits - 1; j++) { + if (diff[j]) sim.X(controls[j]); + } + } + } + + + + // Tests the multi-exp on a uniform superposition + TEST_METHOD(MultiExpWithHTest) { + const int num_qubits = 32; + auto qubit_prep = [](SparseSimulator& sim) { + sim.H(0); + sim.H(1); + sim.H(2); + }; + auto qubit_clear = [](SparseSimulator& sim) { + sim.H(2); + sim.release(2); + sim.H(1); + sim.release(1); + sim.H(0); + sim.release(0); + }; + MultiExpTest(qubit_prep, qubit_clear); + } + + // Tests the MultiExp on all computational basis states of 3 qubits + TEST_METHOD(MultiExpBasisTest) { + const int num_qubits = 32; + auto qubit_prep = [](SparseSimulator& sim, int index) { + if ((index & 1) == 0) { sim.X(0); } + if ((index & 2) == 0) { sim.X(1); } + if ((index & 4) == 0) { sim.X(2); } + }; + auto qubit_clear = [](SparseSimulator& sim, int index) { + if ((index & 1) == 0) { sim.X(0); } + sim.release(0); + if ((index & 2) == 0) { sim.X(1); } + sim.release(1); + if ((index & 4) == 0) { sim.X(2); } + sim.release(2); + }; + for (int i = 0; i < 8; i++) { + MultiExpTest([=](SparseSimulator& sim) {qubit_prep(sim, i); }, [=](SparseSimulator& sim) {qubit_clear(sim, i); }); + } + } + + // Tests a SWAP gate on all 2-qubit computational basis states + TEST_METHOD(SWAPGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + std::vector qubits{ 0,1 }; + sim.SWAP({ qubits[0] }, qubits[1]); // 00 -> 00 + assert_amplitude_equality(sim.probe("00"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 0.0, 0.0); + sim.X(qubits[0]); + sim.SWAP( qubits[0] , qubits[1]); // 10 -> 01 + assert_amplitude_equality(sim.probe("00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 0.0, 0.0); + sim.SWAP( qubits[0] , qubits[1]); // 01 -> 10 + assert_amplitude_equality(sim.probe("0"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 0.0, 0.0); + sim.X(qubits[1]); + sim.SWAP( qubits[0] , qubits[1]); // 11 -> 11 + assert_amplitude_equality(sim.probe("00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("11"), 1.0, 0.0); + } + + // Tests multi-controlled swap on all computational basis states of 4 qubits + // (2 controls, 2 targets) + TEST_METHOD(CSWAPGateTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + std::vector target_qubits{ 0,1 }; + std::vector control_qubits{ 2,3 }; + // Lambda to test when controls should cause no swap + auto no_swap_test = [&](std::string controls) { + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 00 -> 00 + assert_amplitude_equality(sim.probe(controls+"00"), 1.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "11"), 0.0, 0.0); + sim.X(target_qubits[0]); + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 01 -> 01 + assert_amplitude_equality(sim.probe(controls + "00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "01"), 1.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "11"), 0.0, 0.0); + sim.X(target_qubits[1]); + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 11 -> 11 + assert_amplitude_equality(sim.probe(controls + "00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "10"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "11"), 1.0, 0.0); + sim.X(target_qubits[0]); + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 10 -> 10 + assert_amplitude_equality(sim.probe(controls + "00"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "01"), 0.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "10"), 1.0, 0.0); + assert_amplitude_equality(sim.probe(controls + "11"), 0.0, 0.0); + sim.X(target_qubits[1]); + }; + // Controls are 00, no swap + no_swap_test("00"); + sim.X(control_qubits[0]); + // Controls are 01, no swap + no_swap_test("01"); + sim.X(control_qubits[1]); + // Controls are 11, test for swap + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 00 -> 00 + assert_amplitude_equality(sim.probe("1100"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1101"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1110"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1111"), 0.0, 0.0); + sim.X(target_qubits[0]); + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 10 -> 01 + assert_amplitude_equality(sim.probe("1100"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1101"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1110"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1111"), 0.0, 0.0); + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]);// 01 -> 10 + assert_amplitude_equality(sim.probe("1100"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1101"), 1.0, 0.0); + assert_amplitude_equality(sim.probe("1110"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1111"), 0.0, 0.0); + sim.X(target_qubits[1]); + sim.CSWAP(control_qubits, target_qubits[0], target_qubits[1]); // 11 -> 11 + assert_amplitude_equality(sim.probe("1100"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1101"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1110"), 0.0, 0.0); + assert_amplitude_equality(sim.probe("1111"), 1.0, 0.0); + sim.X(target_qubits[1]); + sim.X(target_qubits[0]); + sim.X(control_qubits[0]); + // Controls are 10, test for no swap + no_swap_test("10"); + } + + // Tests measurement probabilistically + // Based on the expected measurement probabilities for a Pauli Y + // rotation + // It samples a lot of measurements, and based on the + // current variance (of a binomial distrbution): + // - if it's very close to the expected distribution, + // it considers this a success + // - if it's very far from the expected distribution, + // it throws an exception + // - if it's in between, it runs more samples + // While this run-time is undetermined, the threshold + // for an exception shrinks with the number of tests + TEST_METHOD(MTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + const int n_tests = 5000; + const double log_false_positive_threshold = 0.1; + const double log_false_negative_threshold = 100.0; + for (double angle = 0.0; angle < M_PI / 2.0; angle += 0.1) { + SparseSimulator sim = SparseSimulator(num_qubits); + sim.set_random_seed(12345); + double expected_ratio = sin(angle / 2.0) * sin(angle / 2.0); + double ratio = 0.0; + unsigned long total_tests = 0; + unsigned long ones = 0; + double std_dev = 0.0; + double log_prob = 0.0; + logical_qubit_id qubit = 0; + do { + for (int i = 0; i < n_tests; i++) { + sim.R(Gates::Basis::PauliY, angle, qubit); + if (sim.M(qubit)) { + ones++; + sim.X(qubit); + } + } + total_tests += n_tests; + ratio = (double)ones / (double)total_tests; + double abs_diff = abs(expected_ratio - ratio); + // Based on Chernoff bounds + log_prob = abs_diff * abs_diff * expected_ratio * (double)total_tests; + std_dev = sqrt(expected_ratio * (1.0 - expected_ratio)) / (double)total_tests; + // Using variance of the binomial distribution + if (log_false_positive_threshold >= log_prob) { + break; + } + } while (log_false_negative_threshold >= log_prob); + if (log_false_negative_threshold < log_prob) { + throw std::runtime_error("Statistically improbable measurement results"); + } + } + + } + + // Tests an assortment of assertions to both pass and to throw exceptions + TEST_METHOD(AssertTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + using namespace Gates; + std::vector basis{ Basis::PauliZ, Basis::PauliZ, Basis::PauliZ }; + std::vector qubits{ 0,1,2 }; + sim.Assert(basis, qubits, false); + sim.update_state(); + // These require forcing the simulator to update the state for it to actually throw the exception + Assert::ExpectException([&]() {sim.Assert(basis, qubits, true); sim.update_state(); }); + basis = { Basis::PauliZ, Basis::PauliZ, Basis:: PauliI }; + sim.Assert(basis, qubits, false); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, true); sim.update_state(); }); + basis = { Basis::PauliX, Basis::PauliI, Basis::PauliI }; + Assert::ExpectException([&]() {sim.Assert(basis, qubits, false); sim.update_state(); }); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, true); sim.update_state(); }); + basis = { Basis::PauliY, Basis::PauliI, Basis::PauliI }; + Assert::ExpectException([&]() {sim.Assert(basis, qubits, false); sim.update_state(); }); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, true); sim.update_state(); }); + } + + // Tests an assortment of assertions on GHZ states + TEST_METHOD(AssertGHZTest) { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + SparseSimulator sim = SparseSimulator(num_qubits); + using namespace Gates; + std::vector basis(3, Basis::PauliX); + std::vector qubits{ 0,1,2 }; + sim.H(0); + sim.MCX({ 0 }, 1); + sim.MCX({ 0 }, 2); + sim.Assert(basis, qubits, false); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, true); }); + sim.Z(0); + sim.Assert(basis, qubits, true); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, false); }); + sim.S(0); + basis = { Basis::PauliY, Basis::PauliY, Basis::PauliY }; + sim.Assert(basis, qubits, false); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, true); }); + sim.Z(0); + sim.Assert(basis, qubits, true); + Assert::ExpectException([&]() {sim.Assert(basis, qubits, false); }); + sim.probe("0"); + } + + // Basic test of quantum teleportation + TEST_METHOD(TeleportationTest) + { + const logical_qubit_id num_qubits = 32; + const qubit_label_type zero(0); + for (double test_angle = 0; test_angle < 1.0; test_angle += 0.34) { + SparseSimulator sim = SparseSimulator(num_qubits); + sim.set_random_seed(12345); + std::vector qubits{ 0,1,2 }; + sim.H(qubits[1]); + sim.MCX({ qubits[1] }, qubits[2]); + + sim.R(Gates::Basis::PauliY, test_angle, 0); + + sim.MCX({ qubits[0] }, qubits[1]); + sim.H(qubits[0]); + bool result0 = sim.M(qubits[0]); + bool result1 = sim.M(qubits[1]); + if (result1) { + sim.X(qubits[2]); + sim.X(qubits[1]); + } + if (result0) { + sim.Z(qubits[2]); + sim.X(qubits[0]); + } + + amplitude teleported_qubit_0 = sim.probe("000"); + amplitude teleported_qubit_1 = sim.probe("100"); + Assert::AreEqual((float)cos(test_angle / 2.0), (float)teleported_qubit_0.real()); + Assert::AreEqual((float)0.0, (float)teleported_qubit_0.imag()); + Assert::AreEqual((float)sin(test_angle / 2.0), (float)teleported_qubit_1.real()); + Assert::AreEqual((float)0.0, (float)teleported_qubit_1.imag()); + } + } + + + // Tests that H gates properly cancel when executed + TEST_METHOD(HCancellationTest) + { + const int n_qubits = 128; + SparseSimulator sim = SparseSimulator(n_qubits); + sim.set_random_seed(12345); + std::vector qubits(n_qubits); + std::generate(qubits.begin(), qubits.end(), [] { static int i{ 0 }; return i++; }); + size_t buckets = 0; + // Will cause a huge memory problem if there is no cancellation + const int n_samples = 16; + for (int i = 0; i < n_qubits; i += n_samples) { + for (int ii = 0; ii < n_samples; ii++) { + sim.H(qubits[i + ii]); + } + sim.update_state(); + for (int ii = n_samples - 1; ii >= 0; ii--) { + sim.H(qubits[i + ii]); + } + sim.update_state(); + } + } + + // Checks that X and Z gates commute with H + TEST_METHOD(HXZCommutationTest) + { + const int n_qubits = 16; + SparseSimulator sim = SparseSimulator(n_qubits); + sim.set_random_seed(12345); + std::vector qubits(n_qubits); + std::generate(qubits.begin(), qubits.end(), [] { static int i{ 0 }; return i++; }); + for (int i = 0; i < n_qubits; i++) { + sim.H(qubits[i]); + } + // Here it will actually just commute the X and Z through the H in the queue + // without actually executing anything + std::bitset one_state = 0; + for (int i = 0; i < n_qubits - 1; i += 2) { + sim.Z(qubits[i]); + sim.X(qubits[i + 1]); + one_state.set(i); + } + for (int i = n_qubits - 1; i >= 0; i--) { + sim.H(qubits[i]); + } + for (__int64 i = 0; i < pow(2, n_qubits); i++) { + amplitude state = sim.probe(std::bitset(i).to_string()); + if (i == one_state.to_ulong()) { + assert_amplitude_equality(state, 1.0, 0.0); + } + else { + assert_amplitude_equality(state, 0.0, 0.0); + } + } + } + + }; + +} diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj new file mode 100644 index 00000000000..4cbcb322744 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj @@ -0,0 +1,179 @@ + + + + + Debug + Win32 + + + Release + Win32 + + + Debug + x64 + + + Release + x64 + + + + 16.0 + {E1663845-BBA8-41EA-AC31-50A092EEA728} + Win32Proj + SparseSimulatorTests + 10.0 + NativeUnitTestProject + + + + DynamicLibrary + true + v142 + Unicode + false + + + DynamicLibrary + false + v142 + true + Unicode + false + + + DynamicLibrary + true + v142 + Unicode + false + + + DynamicLibrary + false + v142 + true + Unicode + false + + + + + + + + + + + + + + + + + + + + + true + + + true + + + false + + + false + + + + NotUsing + Level3 + true + $(VCInstallDir)UnitTest\include;%(AdditionalIncludeDirectories) + WIN32;_DEBUG;%(PreprocessorDefinitions) + true + pch.h + stdcpp17 + + + Windows + $(VCInstallDir)UnitTest\lib;%(AdditionalLibraryDirectories) + + + + + NotUsing + Level3 + true + $(VCInstallDir)UnitTest\include;%(AdditionalIncludeDirectories) + _DEBUG;%(PreprocessorDefinitions) + true + pch.h + stdcpp17 + + + Windows + $(VCInstallDir)UnitTest\lib;%(AdditionalLibraryDirectories) + + + + + NotUsing + Level3 + true + true + true + $(VCInstallDir)UnitTest\include;%(AdditionalIncludeDirectories) + WIN32;NDEBUG;%(PreprocessorDefinitions) + true + pch.h + + + Windows + true + true + $(VCInstallDir)UnitTest\lib;%(AdditionalLibraryDirectories) + + + + + NotUsing + Level3 + true + true + true + $(VCInstallDir)UnitTest\include;%(AdditionalIncludeDirectories) + NDEBUG;%(PreprocessorDefinitions) + true + pch.h + stdcpp17 + + + Windows + true + true + $(VCInstallDir)UnitTest\lib;%(AdditionalLibraryDirectories) + + + + + + Create + Create + Create + Create + + + + + + + + + + + + \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj.filters b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj.filters new file mode 100644 index 00000000000..5947f4294c2 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/SparseSimulatorTests.vcxproj.filters @@ -0,0 +1,39 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;c++;cppm;ixx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hh;hpp;hxx;h++;hm;inl;inc;ipp;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms + + + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + + + Header Files + + + Source Files + + + \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp new file mode 100644 index 00000000000..e2d8135a0df --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.cpp @@ -0,0 +1,160 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#include "pch.h" +#include "CppUnitTest.h" +#include "../Native/SparseSimulator.h" +#include +#include + +using namespace Microsoft::VisualStudio::CppUnitTestFramework; +using namespace Microsoft::Quantum::SPARSESIMULATOR; + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define TEST_TOLERANCE 0.00000001 + + +namespace SparseSimulatorTestHelpers +{ + + inline std::size_t make_mask(std::vector const& qs) + { + std::size_t mask = 0; + for (std::size_t q : qs) + mask = mask | (1ull << q); + return mask; + } + // power of square root of -1 + inline amplitude iExp(int power) + { + using namespace std::literals::complex_literals; + 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; + } + + inline bool isDiagonal(std::vector const& b) + { + for (auto x : b) + if (x == Gates::Basis::PauliX || x == Gates::Basis::PauliY) return false; + return true; + } + + inline static void removeIdentities(std::vector& b, std::vector& qs) + { + unsigned i = 0; + while (i != b.size()) + { + if (b[i] == Gates::Basis::PauliI) + { + b.erase(b.begin() + i); + qs.erase(qs.begin() + i); + } + else + ++i; + } + } + + // Taken from qsharp-runtime + void apply_exp( + std::vector& wfn, + std::vector b, + double phi, + std::vector qs) + { + removeIdentities(b, qs); + if (qs.size() == 0) { return; } + unsigned lowest = *std::min_element(qs.begin(), qs.end()); + + std::size_t offset = 1ull << lowest; + if (isDiagonal(b)) + { + std::size_t mask = make_mask(qs); + amplitude phase = std::exp(amplitude(0., -phi)); + +#pragma omp parallel for schedule(static) + for (std::intptr_t x = 0; x < static_cast(wfn.size()); x++) + wfn[x] *= (std::bitset<64>(x & mask).count() % 2 ? phase : std::conj(phase)); + } + else { + std::size_t xy_bits = 0; + std::size_t yz_bits = 0; + int y_count = 0; + for (unsigned i = 0; i < b.size(); ++i) + { + switch (b[i]) + { + case Gates::Basis::PauliX: + xy_bits |= (1ull << qs[i]); + break; + case Gates::Basis::PauliY: + xy_bits |= (1ull << qs[i]); + yz_bits |= (1ull << qs[i]); + ++y_count; + break; + case Gates::Basis::PauliZ: + yz_bits |= (1ull << qs[i]); + break; + case Gates::Basis::PauliI: + break; + default: + break; + } + } + + amplitude alpha = std::cos(phi); + amplitude beta = std::sin(phi) * iExp(3 * y_count + 1); + amplitude gamma = std::sin(phi) * iExp(y_count + 1); + + for (std::intptr_t x = 0; x < static_cast(wfn.size()); x++) + { + std::intptr_t t = x ^ xy_bits; + if (x < t) + { + bool parity = std::bitset<64>(x & yz_bits).count() % 2; + auto a = wfn[x]; + auto b = wfn[t]; + wfn[x] = alpha * a + (parity ? -beta : beta) * b; + wfn[t] = alpha * b + (parity ? -gamma : gamma) * a; + } + } + } + } + + // Assertions for equality of amplitude types + inline void assert_double_equality_with_tolerance(double value1, double value2) { + if (value1 > value2) { + value1 = std::max(value1 - TEST_TOLERANCE, value2); + } + else { + value1 = std::min(value1 + TEST_TOLERANCE, value2); + } + Assert::AreEqual(value1, value2); + } + + void assert_amplitude_equality(amplitude amp, double real, double imag) { + assert_double_equality_with_tolerance(real, amp.real()); + assert_double_equality_with_tolerance(imag, amp.imag()); + } + + void assert_amplitude_equality(amplitude expected_amp, amplitude actual_amp) { + assert_amplitude_equality(actual_amp, expected_amp.real(), expected_amp.imag()); + } + + +} \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp new file mode 100644 index 00000000000..8eb71c0f7b1 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/TestHelpers.hpp @@ -0,0 +1,37 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +#pragma once + +#include "pch.h" +#include "CppUnitTest.h" +#include "../Native/SparseSimulator.h" +#include +#include + +using namespace Microsoft::VisualStudio::CppUnitTestFramework; +using namespace Microsoft::Quantum::SPARSESIMULATOR; + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + + + + +namespace SparseSimulatorTestHelpers +{ + + void apply_exp( + std::vector& wfn, + std::vector b, + double phi, + std::vector qs); + + void assert_amplitude_equality(amplitude amp, double real, double imag); + + void assert_amplitude_equality(amplitude expected_amp, amplitude actual_amp); + + // Prepares some qubits, then checks whether various Pauli exponentials work + +} \ No newline at end of file diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp new file mode 100644 index 00000000000..2d641a7e3b9 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.cpp @@ -0,0 +1,8 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +// pch.cpp: source file corresponding to the pre-compiled header + +#include "pch.h" + +// When you are using pre-compiled headers, this source file is necessary for compilation to succeed. diff --git a/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h new file mode 100644 index 00000000000..2cd3c85c965 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/SparseSimulatorTests/pch.h @@ -0,0 +1,16 @@ +// 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. +// However, files listed here are ALL re-compiled if any one of them is updated between builds. +// Do not add files here that you will be updating frequently as this negates the performance advantage. + +#ifndef PCH_H +#define PCH_H + +// add headers that you want to pre-compile here + + +#endif //PCH_H diff --git a/src/Simulation/Simulators/SparseSimulator/build.cake b/src/Simulation/Simulators/SparseSimulator/build.cake new file mode 100644 index 00000000000..9b923a0f297 --- /dev/null +++ b/src/Simulation/Simulators/SparseSimulator/build.cake @@ -0,0 +1,72 @@ +#addin nuget:?package=Cake.CMake&version=1.2.0 + +var release_configuration = Argument("configuration", "Release"); +var debug_configuration = Argument("configuration", "Debug"); +var target = Argument("target", "Build"); + +////////////////////////////////////////////////////////////////////// +// TASKS +////////////////////////////////////////////////////////////////////// + +Task("Clean") + .WithCriteria(c => HasArgument("rebuild")) + .Does(() => +{ + CleanDirectory($"./{release_configuration}"); + CleanDirectory($"./{debug_configuration}"); +}); + +Task("CMake" +) .IsDependentOn("Clean") + .Does(() => +{ + CMake(new CMakeSettings { + OutputPath = Directory("Native/build"), + SourcePath = Directory("Native/"), + }); + + CMakeBuild(new CMakeBuildSettings { + BinaryPath = Directory("Native/build"), + Configuration = release_configuration + }); +}); + +Task("Build") + .IsDependentOn("CMake") + .Does(() => +{ + DotNetCoreBuild(".", new DotNetCoreBuildSettings + { + Configuration = release_configuration, + }); +}); + +Task("CMake_debug" +) .IsDependentOn("Clean") + .Does(() => +{ + CMake(new CMakeSettings { + OutputPath = Directory("Native/build"), + SourcePath = Directory("Native/"), + }); + + CMakeBuild(new CMakeBuildSettings { + BinaryPath = Directory("Native/build"), + Configuration = debug_configuration + }); +}); + +Task("Build_debug") + .IsDependentOn("CMake_debug") + .Does(() => +{ + DotNetCoreBuild(".", new DotNetCoreBuildSettings + { + Configuration = debug_configuration, + }); +}); + + +RunTarget(target); +RunTarget(Argument("target", "Build_debug")); +