diff --git a/.clang-format b/.clang-format index a49304f2..f496c47a 100644 --- a/.clang-format +++ b/.clang-format @@ -7,3 +7,5 @@ AccessModifierOffset: -4 BinPackArguments: false BinPackParameters: false ColumnLimit: 100 +DerivePointerAlignment: false +PointerAlignment: Left diff --git a/.github/workflows/ci_ubuntu.yml b/.github/workflows/ci_ubuntu.yml index 2727155c..3f1a341a 100644 --- a/.github/workflows/ci_ubuntu.yml +++ b/.github/workflows/ci_ubuntu.yml @@ -65,6 +65,11 @@ jobs: - name: Test in Ubuntu run: | OMP_PROC_BIND=false ninja test -C build -j $(nproc) + + - name: Test if stub exists + run: | + echo -e "from scaluq import StateVector\nfrom scaluq.gate import I" > /tmp/stub_sample.py + mypy /tmp/stub_sample.py nvcc-build: name: NVCC build diff --git a/README.md b/README.md index 8074e3ba..31e5755b 100644 --- a/README.md +++ b/README.md @@ -81,7 +81,7 @@ int main() { circuit.add_gate(scaluq::gate::X(0)); circuit.add_gate(scaluq::gate::CNot(0, 1)); circuit.add_gate(scaluq::gate::Y(1)); - circuit.add_gate(scaluq::gate::RX(1, M_PI / 2)); + circuit.add_gate(scaluq::gate::RX(1, std::numbers::pi / 2)); circuit.update_quantum_state(state); scaluq::Operator observable(n_qubits); diff --git a/exe/main.cpp b/exe/main.cpp index 7354deb2..68b52618 100644 --- a/exe/main.cpp +++ b/exe/main.cpp @@ -1,20 +1,40 @@ -#include - -#include #include #include #include #include #include "../scaluq/all.hpp" -#include "../scaluq/util/utility.hpp" using namespace scaluq; using namespace std; void run() { - std::uint64_t n_qubits = 5; - auto state = StateVector::Haar_random_state(n_qubits); + auto y_gate = gate::Y(2); + std::cout << y_gate->to_string() << "\n\n"; + + auto cx_gate = gate::CX(0, 2); + std::cout << cx_gate << "\n\n"; + + auto swap_gate = gate::Swap(2, 3, {4, 6}); + std::cout << swap_gate << "\n\n"; + + auto rx_gate = gate::RX(2, 0.5); + std::cout << rx_gate << "\n\n"; + + auto prob_gate = gate::Probablistic({0.1, 0.1, 0.8}, {cx_gate, y_gate, swap_gate}); + std::cout << prob_gate << "\n\n"; + + auto prob_prob_gate = gate::Probablistic({0.5, 0.5}, {cx_gate, prob_gate}); + std::cout << prob_prob_gate << "\n\n"; + + auto prx_gate = gate::ParamRX(2); + std::cout << prx_gate << "\n\n"; + + auto pry_gate = gate::ParamRY(2, 2.5, {1, 3}); + std::cout << pry_gate << "\n\n"; + + auto pprob_gate = gate::ParamProbablistic({0.7, 0.3}, {prx_gate, pry_gate}); + std::cout << pprob_gate << std::endl; } int main() { diff --git a/exe/main.py b/exe/main.py index 349abafb..e67893b9 100644 --- a/exe/main.py +++ b/exe/main.py @@ -1,15 +1,11 @@ from scaluq import * +# from scaluq.gate import S def main(): state = StateVector(2) - gate = SWAP(0, 1) - gate = SWAPGate(gate) - mat = gate.get_matrix() - print(mat) - t = type(mat[0][0]) - print(t) + swap_gate = gate.Swap(0, 1) + print(state) + print(swap_gate) if __name__ == "__main__": - initialize(InitializationSettings().set_num_threads(8)) main() - finalize() diff --git a/pyproject.toml b/pyproject.toml index d36eff54..939ee578 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,7 @@ homepage = "http://www.scaluq.org" [project.optional-dependencies] dev = [ + "mypy == 1.11.2", "scikit-build == 0.17.6", "typing_extensions == 4.12.0", "numpy == 1.26.0", @@ -56,6 +57,7 @@ dev = [ ] ci = [ + "mypy == 1.11.2", "scikit-build == 0.17.6", "typing_extensions == 4.12.0", "numpy == 1.26.0", diff --git a/python/binding.cpp b/python/binding.cpp index 486d87ae..dcd732ae 100644 --- a/python/binding.cpp +++ b/python/binding.cpp @@ -13,10 +13,13 @@ #include #include -#include - namespace nb = nanobind; using namespace nb::literals; + +#define SCALUQ_USE_NANOBIND + +#include + using namespace scaluq; using namespace std::string_literals; @@ -28,7 +31,7 @@ struct type_caster> { NB_TYPE_CASTER(Kokkos::complex, const_name("complex")) template - bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept { + bool from_python(handle src, uint8_t flags, cleanup_list* cleanup) noexcept { (void)flags; (void)cleanup; @@ -60,7 +63,7 @@ struct type_caster> { } template - static handle from_cpp(T2 &&value, rv_policy policy, cleanup_list *cleanup) noexcept { + static handle from_cpp(T2&& value, rv_policy policy, cleanup_list* cleanup) noexcept { (void)policy; (void)cleanup; @@ -76,862 +79,30 @@ void cleanup() { } NB_MODULE(scaluq_core, m) { - m.def("finalize", - &finalize, - "Terminate the Kokkos execution environment. Release the resources.\n\n.. note:: " - "Finalization fails if there exists `StateVector` allocated. You must use " - "`StateVector` only inside inner scopes than the usage of `finalize` or delete all of " - "existing `StateVector`.\n\n.. note:: This is " - "automatically called when the program exits. If you call this manually, you cannot use " - "most of scaluq's functions until the program exits."); - m.def("is_finalized", &is_initialized, "Return true if `finalize()` is already called."); - - nb::class_(m, - "StateVector", - "Vector representation of quantum state.\n\n.. note:: Qubit index is " - "start from 0. If the amplitudes of $\\ket{b_{n-1}\\dots b_0}$ is " - "$b_i$, the state is $\\sum_i b_i 2^i$.") - .def(nb::init(), - "Construct state vector with specified qubits, initialized with computational " - "basis $\\ket{0\\dots0}$.") - .def(nb::init(), "Constructing state vector by copying other state.") - .def_static( - "Haar_random_state", - [](std::uint64_t n_qubits, std::optional seed) { - return StateVector::Haar_random_state(n_qubits, - seed.value_or(std::random_device{}())); - }, - "n_qubits"_a, - "seed"_a = std::nullopt, - "Constructing state vector with Haar random state. If seed is not specified, the value " - "from random device is used.") - .def("set_amplitude_at_index", - &StateVector::set_amplitude_at_index, - "Manually set amplitude at one index.") - .def("get_amplitude_at_index", - &StateVector::get_amplitude_at_index, - "Get amplitude at one index.\n\n.. note:: If you want to get all amplitudes, you " - "should " - "use `StateVector::amplitudes()`.") - .def("set_zero_state", - &StateVector::set_zero_state, - "Initialize with computational basis $\\ket{00\\dots0}$.") - .def("set_zero_norm_state", - &StateVector::set_zero_norm_state, - "Initialize with 0 (null vector).") - .def("set_computational_basis", - &StateVector::set_computational_basis, - "Initialize with computational basis \\ket{\\mathrm{basis}}.") - .def("amplitudes", &StateVector::amplitudes, "Get all amplitudes with as `list[complex]`.") - .def("n_qubits", &StateVector::n_qubits, "Get num of qubits.") - .def("dim", &StateVector::dim, "Get dimension of the vector ($=2^\\mathrm{n\\_qubits}$).") - .def("get_squared_norm", - &StateVector::get_squared_norm, - "Get squared norm of the state. $\\braket{\\psi|\\psi}$.") - .def("normalize", - &StateVector::normalize, - "Normalize state (let $\\braket{\\psi|\\psi} = 1$ by multiplying coef).") - .def("get_zero_probability", - &StateVector::get_zero_probability, - "Get the probability to observe $\\ket{0}$ at specified index.") - .def("get_marginal_probability", - &StateVector::get_marginal_probability, - "Get the marginal probability to observe as specified. Specify the result as n-length " - "list. `0` and `1` represent the qubit is observed and get the value. `2` represents " - "the qubit is not observed.") - .def("get_entropy", &StateVector::get_entropy, "Get the entropy of the vector.") - .def("add_state_vector", - &StateVector::add_state_vector, - "Add other state vector and make superposition. $\\ket{\\mathrm{this}} " - "\\leftarrow " - "\\ket{\\mathrm{this}} + \\ket{\\mathrm{state}}$.") - .def("add_state_vector_with_coef", - &StateVector::add_state_vector_with_coef, - "add other state vector with multiplying the coef and make superposition. " - "$\\ket{\\mathrm{this}}\\leftarrow\\ket{\\mathrm{this}}+\\mathrm{coef}" - "\\ket{\\mathrm{" - "state}}$.") - .def("multiply_coef", - &StateVector::multiply_coef, - "Multiply coef. " - "$\\ket{\\mathrm{this}}\\leftarrow\\mathrm{coef}\\ket{\\mathrm{this}}$.") - .def( - "sampling", - [](const StateVector &state, - std::uint64_t sampling_count, - std::optional seed) { - return state.sampling(sampling_count, seed.value_or(std::random_device{}())); - }, - "sampling_count"_a, - "seed"_a = std::nullopt, - "Sampling specified times. Result is `list[int]` with the `sampling_count` length.") - .def("to_string", &StateVector::to_string, "Information as `str`.") - .def("load", &StateVector::load, "Load amplitudes of `list[int]` with `dim` length.") - .def("__str__", &StateVector::to_string, "Information as `str`.") - .def_ro_static("UNMEASURED", - &StateVector::UNMEASURED, - "Constant used for `StateVector::get_marginal_probability` to express the " - "the qubit is not measured."); - - nb::class_( - m, - "StateVectorBatched", - "Batched vector representation of quantum state.\n\n.. note:: Qubit index is start from 0. " - "If the amplitudes of $\\ket{b_{n-1}\\dots b_0}$ is $b_i$, the state is $\\sum_i b_i " - "2^i$.") - .def(nb::init(), - "Construct batched state vector with specified batch size and qubits.") - .def(nb::init(), - "Constructing batched state vector by copying other batched state.") - .def("n_qubits", &StateVectorBatched::n_qubits, "Get num of qubits.") - .def("dim", - &StateVectorBatched::dim, - "Get dimension of the vector ($=2^\\mathrm{n\\_qubits}$).") - .def("batch_size", &StateVectorBatched::batch_size, "Get batch size.") - .def("set_state_vector", - nb::overload_cast(&StateVectorBatched::set_state_vector), - "Set the state vector for all batches.") - .def("set_state_vector", - nb::overload_cast( - &StateVectorBatched::set_state_vector), - "Set the state vector for a specific batch.") - .def("get_state_vector", - &StateVectorBatched::get_state_vector, - "Get the state vector for a specific batch.") - .def("set_zero_state", - &StateVectorBatched::set_zero_state, - "Initialize all batches with computational basis $\\ket{00\\dots0}$.") - .def("set_zero_norm_state", - &StateVectorBatched::set_zero_norm_state, - "Initialize with 0 (null vector).") - .def("set_computational_basis", - &StateVectorBatched::set_computational_basis, - "Initialize with computational basis \\ket{\\mathrm{basis}}.") - .def( - "sampling", - [](const StateVectorBatched &states, - std::uint64_t sampling_count, - std::optional seed) { - return states.sampling(sampling_count, seed.value_or(std::random_device{}())); - }, - "sampling_count"_a, - "seed"_a = std::nullopt, - "Sampling specified times. Result is `list[list[int]]` with the `sampling_count` " - "length.") - .def_static( - "Haar_random_states", - [](std::uint64_t batch_size, - std::uint64_t n_qubits, - bool set_same_state, - std::optional seed) { - return StateVectorBatched::Haar_random_states( - batch_size, n_qubits, set_same_state, seed.value_or(std::random_device{}())); - }, - "batch_size"_a, - "n_qubits"_a, - "set_same_state"_a, - "seed"_a = std::nullopt, - "Construct batched state vectors with Haar random states. If seed is not " - "specified, the value from random device is used.") - .def("amplitudes", - &StateVectorBatched::amplitudes, - "Get all amplitudes with as `list[list[complex]]`.") - .def("get_squared_norm", - &StateVectorBatched::get_squared_norm, - "Get squared norm of each state in the batch. $\\braket{\\psi|\\psi}$.") - .def("normalize", - &StateVectorBatched::normalize, - "Normalize each state in the batch (let $\\braket{\\psi|\\psi} = 1$ by " - "multiplying coef).") - .def("get_zero_probability", - &StateVectorBatched::get_zero_probability, - "Get the probability to observe $\\ket{0}$ at specified index for each state in " - "the batch.") - .def("get_marginal_probability", - &StateVectorBatched::get_marginal_probability, - "Get the marginal probability to observe as specified for each state in the batch. " - "Specify the result as n-length list. `0` and `1` represent the qubit is observed " - "and get the value. `2` represents the qubit is not observed.") - .def("get_entropy", - &StateVectorBatched::get_entropy, - "Get the entropy of each state in the batch.") - .def("add_state_vector", - &StateVectorBatched::add_state_vector, - "Add other batched state vectors and make superposition. $\\ket{\\mathrm{this}} " - "\\leftarrow \\ket{\\mathrm{this}} + \\ket{\\mathrm{states}}$.") - .def("add_state_vector_with_coef", - &StateVectorBatched::add_state_vector_with_coef, - "Add other batched state vectors with multiplying the coef and make superposition. " - "$\\ket{\\mathrm{this}}\\leftarrow\\ket{\\mathrm{this}}+\\mathrm{coef}" - "\\ket{\\mathrm{states}}$.") - .def("load", - &StateVectorBatched::load, - "Load batched amplitudes from `list[list[complex]]`.") - .def("copy", &StateVectorBatched::copy, "Create a copy of the batched state vector.") - .def("to_string", &StateVectorBatched::to_string, "Information as `str`.") - .def("__str__", &StateVectorBatched::to_string, "Information as `str`."); - - nb::enum_(m, "GateType", "Enum of Gate Type.") - .value("I", GateType::I) - .value("GlobalPhase", GateType::GlobalPhase) - .value("X", GateType::X) - .value("Y", GateType::Y) - .value("Z", GateType::Z) - .value("H", GateType::H) - .value("S", GateType::S) - .value("Sdag", GateType::Sdag) - .value("T", GateType::T) - .value("Tdag", GateType::Tdag) - .value("SqrtX", GateType::SqrtX) - .value("SqrtXdag", GateType::SqrtXdag) - .value("SqrtY", GateType::SqrtY) - .value("SqrtYdag", GateType::SqrtYdag) - .value("P0", GateType::P0) - .value("P1", GateType::P1) - .value("RX", GateType::RX) - .value("RY", GateType::RY) - .value("RZ", GateType::RZ) - .value("U1", GateType::U1) - .value("U2", GateType::U2) - .value("U3", GateType::U3) - .value("OneTargetMatrix", GateType::OneTargetMatrix) - .value("Swap", GateType::Swap) - .value("TwoTargetMatrix", GateType::TwoTargetMatrix) - .value("Pauli", GateType::Pauli) - .value("PauliRotation", GateType::PauliRotation); - -#define DEF_GATE_BASE(GATE_TYPE, DESCRIPTION) \ - nb::class_(m, #GATE_TYPE, DESCRIPTION) \ - .def("gate_type", &GATE_TYPE::gate_type, "Get gate type as `GateType` enum.") \ - .def( \ - "get_target_qubit_list", \ - [](const GATE_TYPE &gate) { return gate->get_target_qubit_list(); }, \ - "Get target qubits as `list[int]`. **Control qubits is not included.**") \ - .def( \ - "get_control_qubit_list", \ - [](const GATE_TYPE &gate) { return gate->get_control_qubit_list(); }, \ - "Get control qubits as `list[int]`.") \ - .def( \ - "get_operand_qubit_list", \ - [](const GATE_TYPE &gate) { return gate->get_operand_qubit_list(); }, \ - "Get target and control qubits as `list[int]`.") \ - .def( \ - "get_target_qubit_mask", \ - [](const GATE_TYPE &gate) { return gate->get_target_qubit_mask(); }, \ - "Get target qubits as mask. **Control qubits is not included.**") \ - .def( \ - "get_control_qubit_mask", \ - [](const GATE_TYPE &gate) { return gate->get_control_qubit_mask(); }, \ - "Get control qubits as mask.") \ - .def( \ - "get_operand_qubit_mask", \ - [](const GATE_TYPE &gate) { return gate->get_operand_qubit_mask(); }, \ - "Get target and control qubits as mask.") \ - .def( \ - "get_inverse", \ - [](const GATE_TYPE &gate) { return gate->get_inverse(); }, \ - "Generate inverse gate as `Gate` type. If not exists, return None.") \ - .def( \ - "update_quantum_state", \ - [](const GATE_TYPE &gate, StateVector &state_vector) { \ - gate->update_quantum_state(state_vector); \ - }, \ - "Apply gate to `state_vector`. `state_vector` in args is directly updated.") \ - .def( \ - "get_matrix", \ - [](const GATE_TYPE &gate) { return gate->get_matrix(); }, \ - "Get matrix representation of the gate.") - -#define DEF_GATE(GATE_TYPE, DESCRIPTION) \ - DEF_GATE_BASE( \ - GATE_TYPE, \ - DESCRIPTION \ - "\n\n.. note:: Upcast is required to use gate-general functions (ex: add to Circuit).") \ - .def(nb::init()) - - DEF_GATE_BASE(Gate, - "General class of QuantumGate.\n\n.. note:: Downcast to requred to use " - "gate-specific functions.") - .def(nb::init(), "Upcast from `IGate`.") - .def(nb::init(), "Upcast from `GlobalPhaseGate`.") - .def(nb::init(), "Upcast from `XGate`.") - .def(nb::init(), "Upcast from `YGate`.") - .def(nb::init(), "Upcast from `ZGate`.") - .def(nb::init(), "Upcast from `HGate`.") - .def(nb::init(), "Upcast from `SGate`.") - .def(nb::init(), "Upcast from `SdagGate`.") - .def(nb::init(), "Upcast from `TGate`.") - .def(nb::init(), "Upcast from `TdagGate`.") - .def(nb::init(), "Upcast from `SqrtXGate`.") - .def(nb::init(), "Upcast from `SqrtXdagGate`.") - .def(nb::init(), "Upcast from `SqrtYGate`.") - .def(nb::init(), "Upcast from `SqrtYdagGate`.") - .def(nb::init(), "Upcast from `P0Gate`.") - .def(nb::init(), "Upcast from `P1Gate`.") - .def(nb::init(), "Upcast from `RXGate`.") - .def(nb::init(), "Upcast from `RYGate`.") - .def(nb::init(), "Upcast from `RZGate`.") - .def(nb::init(), "Upcast from `U1Gate`.") - .def(nb::init(), "Upcast from `U2Gate`.") - .def(nb::init(), "Upcast from `U3Gate`.") - .def(nb::init(), "Upcast from `OneTargetMatrixGate`.") - .def(nb::init(), "Upcast from `SwapGate`.") - .def(nb::init(), "Upcast from `TwoTargetMatrixGate`.") - .def(nb::init(), "Upcast from `PauliGate`.") - .def(nb::init(), "Upcast from `PauliRotationGate`."); - - DEF_GATE(IGate, "Specific class of Pauli-I gate."); - DEF_GATE(GlobalPhaseGate, - "Specific class of gate, which rotate global phase, represented as " - "$e^{i\\mathrm{phase}}I$.") - .def( - "phase", - [](const GlobalPhaseGate &gate) { return gate->phase(); }, - "Get `phase` property"); - DEF_GATE(XGate, "Specific class of Pauli-X gate."); - DEF_GATE(YGate, "Specific class of Pauli-Y gate."); - DEF_GATE(ZGate, "Specific class of Pauli-Z gate."); - DEF_GATE(HGate, "Specific class of Hadamard gate."); - DEF_GATE(SGate, - "Specific class of S gate, represented as $\\begin{bmatrix}\n1 & 0\\\\\n0 & " - "i\n\\end{bmatrix}$."); - DEF_GATE(SdagGate, "Specific class of inverse of S gate."); - DEF_GATE(TGate, - "Specific class of T gate, represented as $\\begin{bmatrix}\n1 & 0\\\\\n0 & " - "e^{i\\pi/4}\n\\end{bmatrix}$."); - DEF_GATE(TdagGate, "Specific class of inverse of T gate."); - DEF_GATE(SqrtXGate, - "Specific class of sqrt(X) gate, represented as $\\begin{bmatrix}\n1+i & 1-i\\\\\n1-i " - "& 1+i\n\\end{bmatrix}$."); - DEF_GATE(SqrtXdagGate, "Specific class of inverse of sqrt(X) gate."); - DEF_GATE(SqrtYGate, - "Specific class of sqrt(Y) gate, represented as $\\begin{bmatrix}\n1+i & -1-i " - "\\\\\n1+i & 1+i\n\\end{bmatrix}$."); - DEF_GATE(SqrtYdagGate, "Specific class of inverse of sqrt(Y) gate."); - DEF_GATE( - P0Gate, - "Specific class of projection gate to $\\ket{0}$.\n\n.. note:: This gate is not unitary."); - DEF_GATE( - P1Gate, - "Specific class of projection gate to $\\ket{1}$.\n\n.. note:: This gate is not unitary."); - -#define DEF_ROTATION_GATE(GATE_TYPE, DESCRIPTION) \ - DEF_GATE(GATE_TYPE, DESCRIPTION) \ - .def( \ - "angle", [](const GATE_TYPE &gate) { return gate->angle(); }, "Get `angle` property.") - - DEF_ROTATION_GATE( - RXGate, - "Specific class of X rotation gate, represented as $e^{-i\\frac{\\mathrm{angle}}{2}X}$."); - DEF_ROTATION_GATE( - RYGate, - "Specific class of Y rotation gate, represented as $e^{-i\\frac{\\mathrm{angle}}{2}Y}$."); - DEF_ROTATION_GATE( - RZGate, - "Specific class of Z rotation gate, represented as $e^{-i\\frac{\\mathrm{angle}}{2}Z}$."); - - DEF_GATE(U1Gate, - "Specific class of IBMQ's U1 Gate, which is a rotation abount Z-axis, " - "represented as " - "$\\begin{bmatrix}\n1 & 0\\\\\n0 & e^{i\\lambda}\n\\end{bmatrix}$.") - .def( - "lambda_", [](const U1Gate &gate) { return gate->lambda(); }, "Get `lambda` property."); - DEF_GATE(U2Gate, - "Specific class of IBMQ's U2 Gate, which is a rotation about X+Z-axis, " - "represented as " - "$\\frac{1}{\\sqrt{2}} \\begin{bmatrix}1 & -e^{-i\\lambda}\\\\\n" - "e^{i\\phi} & e^{i(\\phi+\\lambda)}\n\\end{bmatrix}$.") - .def( - "phi", [](const U2Gate &gate) { return gate->phi(); }, "Get `phi` property.") - .def( - "lambda_", [](const U2Gate &gate) { return gate->lambda(); }, "Get `lambda` property."); - DEF_GATE(U3Gate, - "Specific class of IBMQ's U3 Gate, which is a rotation abount 3 axis, " - "represented as " - "$\\begin{bmatrix}\n\\cos \\frac{\\theta}{2} & " - "-e^{i\\lambda}\\sin\\frac{\\theta}{2}\\\\\n" - "e^{i\\phi}\\sin\\frac{\\theta}{2} & " - "e^{i(\\phi+\\lambda)}\\cos\\frac{\\theta}{2}\n\\end{bmatrix}$.") - .def( - "theta", [](const U3Gate &gate) { return gate->theta(); }, "Get `theta` property.") - .def( - "phi", [](const U3Gate &gate) { return gate->phi(); }, "Get `phi` property.") - .def( - "lambda_", [](const U3Gate &gate) { return gate->lambda(); }, "Get `lambda` property."); - - DEF_GATE(OneTargetMatrixGate, "Specific class of single-qubit dense matrix gate.") - .def("matrix", [](const OneTargetMatrixGate &gate) { return gate->matrix(); }); - - DEF_GATE(PauliGate, - "Specific class of multi-qubit pauli gate, which applies single-qubit Pauli " - "gate to " - "each of qubit."); - DEF_GATE(PauliRotationGate, - "Specific class of multi-qubit pauli-rotation gate, represented as " - "$e^{-i\\frac{\\mathrm{angle}}{2}P}$."); + internal::bind_types_hpp(m); - DEF_GATE(ProbablisticGate, - "Specific class of probablistic gate. The gate to apply is picked from a cirtain " - "distribution.") - .def( - "gate_list", - [](const ProbablisticGate &gate) { return gate->gate_list(); }, - nb::rv_policy::reference) - .def( - "distribution", - [](const ProbablisticGate &gate) { return gate->distribution(); }, - nb::rv_policy::reference); + internal::bind_state_state_vector_hpp(m); + internal::bind_state_state_vector_batched_hpp(m); auto mgate = m.def_submodule("gate", "Define gates."); -#define DEF_GATE_FACTORY(GATE_NAME) \ - mgate.def( \ - #GATE_NAME, &gate::GATE_NAME, "Generate general Gate class instance of " #GATE_NAME ".") + internal::bind_gate_gate_hpp(m); + internal::bind_gate_gate_standard_hpp(m); + internal::bind_gate_gate_matrix_hpp(m); + internal::bind_gate_gate_pauli_hpp(m); + internal::bind_gate_gate_factory_hpp(mgate); + internal::bind_gate_param_gate_hpp(m); + internal::bind_gate_param_gate_standard_hpp(m); + internal::bind_gate_param_gate_pauli_hpp(m); + internal::bind_gate_param_gate_probablistic_hpp(m); + internal::bind_gate_param_gate_factory(mgate); + // internal::bind_gate_merge_gate_hpp(m); - DEF_GATE_FACTORY(I); - DEF_GATE_FACTORY(GlobalPhase); - DEF_GATE_FACTORY(X); - DEF_GATE_FACTORY(Y); - DEF_GATE_FACTORY(Z); - DEF_GATE_FACTORY(H); - DEF_GATE_FACTORY(S); - DEF_GATE_FACTORY(Sdag); - DEF_GATE_FACTORY(T); - DEF_GATE_FACTORY(Tdag); - DEF_GATE_FACTORY(SqrtX); - DEF_GATE_FACTORY(SqrtXdag); - DEF_GATE_FACTORY(SqrtY); - DEF_GATE_FACTORY(SqrtYdag); - DEF_GATE_FACTORY(P0); - DEF_GATE_FACTORY(P1); - DEF_GATE_FACTORY(RX); - DEF_GATE_FACTORY(RY); - DEF_GATE_FACTORY(RZ); - DEF_GATE_FACTORY(U1); - DEF_GATE_FACTORY(U2); - DEF_GATE_FACTORY(U3); - DEF_GATE_FACTORY(OneTargetMatrix); - mgate.def( - "CX", - &gate::CX, - "Generate general Gate class instance of CX.\n\n.. note:: CX is a specialization of X."); - mgate.def("CNot", - &gate::CX, - "Generate general Gate class instance of CNot.\n\n.. note:: CNot is an alias of CX."); - mgate.def( - "CZ", - &gate::CZ, - "Generate general Gate class instance of CZ.\n\n.. note:: CZ is a specialization of Z."); - mgate.def( - "CCX", - &gate::CCX, - "Generate general Gate class instance of CXX.\n\n.. note:: CX is a specialization of X."); - mgate.def( - "CCNot", - &gate::CCX, - "Generate general Gate class instance of CCNot.\n\n.. note:: CCNot is an alias of CCX."); - mgate.def("Toffoli", - &gate::CCX, - "Generate general Gate class instance of Toffoli.\n\n.. note:: Toffoli is an alias " - "of CCX."); - DEF_GATE_FACTORY(Swap); - DEF_GATE_FACTORY(TwoTargetMatrix); - DEF_GATE_FACTORY(Pauli); - DEF_GATE_FACTORY(PauliRotation); - mgate.def("DenseMatrix", - &gate::DenseMatrix, - "Generate general Gate class instance of DenseMatrix. IGate, OneTargetMatrixGate or " - "TwoTargetMatrixGate correspond to len(target) is created. The case len(target) >= 3 " - "is currently not supported."); - DEF_GATE_FACTORY(Probablistic); + internal::bind_circuit_circuit_hpp(m); - nb::enum_(m, "ParamGateType", "Enum of ParamGate Type.") - .value("PRX", ParamGateType::PRX) - .value("PRY", ParamGateType::PRY) - .value("PRZ", ParamGateType::PRZ) - .value("PPauliRotation", ParamGateType::PPauliRotation); - - m.def( - "merge_gate", &merge_gate, "Merge two gates. return value is (merged gate, global phase)."); - -#define DEF_PGATE_BASE(PGATE_TYPE, DESCRIPTION) \ - nb::class_(m, #PGATE_TYPE, DESCRIPTION) \ - .def("param_gate_type", \ - &PGATE_TYPE::param_gate_type, \ - "Get parametric gate type as `ParamGateType` enum.") \ - .def( \ - "get_target_qubit_list", \ - [](const PGATE_TYPE &gate) { return gate->get_target_qubit_list(); }, \ - "Get target qubits as `list[int]`. **Control qubits is not included.**") \ - .def( \ - "get_control_qubit_list", \ - [](const PGATE_TYPE &gate) { return gate->get_control_qubit_list(); }, \ - "Get control qubits as `list[int]`.") \ - .def( \ - "get_operand_qubit_list", \ - [](const PGATE_TYPE &gate) { return gate->get_operand_qubit_list(); }, \ - "Get target and control qubits as `list[int]`.") \ - .def( \ - "get_target_qubit_mask", \ - [](const PGATE_TYPE &gate) { return gate->get_target_qubit_mask(); }, \ - "Get target qubits as mask. **Control qubits is not included.**") \ - .def( \ - "get_control_qubit_mask", \ - [](const PGATE_TYPE &gate) { return gate->get_control_qubit_mask(); }, \ - "Get control qubits as mask.") \ - .def( \ - "get_operand_qubit_mask", \ - [](const PGATE_TYPE &gate) { return gate->get_operand_qubit_mask(); }, \ - "Get target and control qubits as mask.") \ - .def( \ - "get_inverse", \ - [](const PGATE_TYPE ¶m_gate) { return param_gate->get_inverse(); }, \ - "Generate inverse parametric-gate as `ParamGate` type. If not exists, return None.") \ - .def( \ - "update_quantum_state", \ - [](const PGATE_TYPE ¶m_gate, StateVector &state_vector, double param) { \ - param_gate->update_quantum_state(state_vector, param); \ - }, \ - "Apply gate to `state_vector` with holding the parameter. `state_vector` in args is " \ - "directly updated.") \ - .def( \ - "get_matrix", \ - [](const PGATE_TYPE &gate, double param) { return gate->get_matrix(param); }, \ - "Get matrix representation of the gate with holding the parameter.") - -#define DEF_PGATE(PGATE_TYPE, DESCRIPTION) \ - DEF_PGATE_BASE( \ - PGATE_TYPE, \ - DESCRIPTION \ - "\n\n.. note:: Upcast is required to use gate-general functions (ex: add to Circuit).") \ - .def(nb::init()) - - DEF_PGATE_BASE( - ParamGate, - "General class of parametric quantum gate.\n\n.. note:: Downcast to requred to use " - "gate-specific functions.") - .def(nb::init(), "Upcast from `PRXGate`.") - .def(nb::init(), "Upcast from `PRYGate`.") - .def(nb::init(), "Upcast from `PRZGate`.") - .def(nb::init(), "Upcast from `PPauliRotationGate`."); - - DEF_PGATE(PRXGate, - "Specific class of parametric X rotation gate, represented as " - "$e^{-i\\frac{\\mathrm{angle}}{2}X}$. `angle` is given as `param * pcoef`."); - DEF_PGATE(PRYGate, - "Specific class of parametric Y rotation gate, represented as " - "$e^{-i\\frac{\\mathrm{angle}}{2}Y}$. `angle` is given as `param * pcoef`."); - DEF_PGATE(PRZGate, - "Specific class of parametric Z rotation gate, represented as " - "$e^{-i\\frac{\\mathrm{angle}}{2}Z}$. `angle` is given as `param * pcoef`."); - - DEF_PGATE(PPauliRotationGate, - "Specific class of parametric multi-qubit pauli-rotation gate, represented as " - "$e^{-i\\frac{\\mathrm{angle}}{2}P}$. `angle` is given as `param * pcoef`."); - - DEF_PGATE(PProbablisticGate, - "Specific class of parametric probablistic gate. The gate to apply is picked from a " - "cirtain " - "distribution.") - .def( - "gate_list", - [](const PProbablisticGate &gate) { return gate->gate_list(); }, - nb::rv_policy::reference) - .def( - "distribution", - [](const PProbablisticGate &gate) { return gate->distribution(); }, - nb::rv_policy::reference); - - mgate.def("PRX", - &gate::PRX, - "Generate general ParamGate class instance of PRX.", - "target"_a, - "coef"_a = 1., - "controls"_a = std::vector{}); - mgate.def("PRY", - &gate::PRY, - "Generate general ParamGate class instance of PRY.", - "target"_a, - "coef"_a = 1., - "controls"_a = std::vector{}); - mgate.def("PRZ", - &gate::PRZ, - "Generate general ParamGate class instance of PRZ.", - "target"_a, - "coef"_a = 1., - "controls"_a = std::vector{}); - mgate.def("PPauliRotation", - &gate::PPauliRotation, - "Generate general ParamGate class instance of PPauliRotation.", - "pauli"_a, - "coef"_a = 1., - "controls"_a = std::vector{}); - mgate.def("PProbablistic", - &gate::PProbablistic, - "Generate general ParamGate class instance of PProbablistic."); - mgate.def( - "PProbablistic", - [](const std::vector>> &prob_gate_list) { - std::vector distribution; - std::vector> gate_list; - distribution.reserve(prob_gate_list.size()); - gate_list.reserve(prob_gate_list.size()); - for (const auto &[prob, gate] : prob_gate_list) { - distribution.push_back(prob); - gate_list.push_back(gate); - } - return gate::PProbablistic(distribution, gate_list); - }, - "Generate general ParamGate class instance of PProbablistic."); - - nb::class_(m, "Circuit", "Quantum circuit represented as gate array") - .def(nb::init(), "Initialize empty circuit of specified qubits.") - .def("n_qubits", &Circuit::n_qubits, "Get property of `n_qubits`.") - .def("gate_list", - &Circuit::gate_list, - "Get property of `gate_list`.", - nb::rv_policy::reference) - .def("gate_count", &Circuit::gate_count, "Get property of `gate_count`.") - .def("key_set", &Circuit::key_set, "Get set of keys of parameters.") - .def("get", &Circuit::get, "Get reference of i-th gate.") - .def("get_key", - &Circuit::get_key, - "Get parameter key of i-th gate. If it is not parametric, return None.") - .def("calculate_depth", &Circuit::calculate_depth, "Get depth of circuit.") - .def("add_gate", - nb::overload_cast(&Circuit::add_gate), - "Add gate. Given gate is copied.") - .def("add_param_gate", - nb::overload_cast(&Circuit::add_param_gate), - "Add parametric gate with specifing key. Given param_gate is copied.") - .def("add_circuit", - nb::overload_cast(&Circuit::add_circuit), - "Add all gates in specified circuit. Given gates are copied.") - .def("update_quantum_state", - &Circuit::update_quantum_state, - "Apply gate to the StateVector. StateVector in args is directly updated. If the " - "circuit contains parametric gate, you have to give real value of parameter as " - "dict[str, float] in 2nd arg.") - .def( - "update_quantum_state", - [&](const Circuit &circuit, StateVector &state, nb::kwargs kwargs) { - std::map parameters; - for (auto &&[key, param] : kwargs) { - parameters[nb::cast(key)] = nb::cast(param); - } - circuit.update_quantum_state(state, parameters); - }, - "Apply gate to the StateVector. StateVector in args is directly updated. If the " - "circuit contains parametric gate, you have to give real value of parameter as " - "\"name=value\" format in kwargs.") - .def( - "update_quantum_state", - [](const Circuit &circuit, StateVector &state) { circuit.update_quantum_state(state); }) - .def("copy", &Circuit::copy, "Copy circuit. All the gates inside is copied.") - .def("get_inverse", - &Circuit::get_inverse, - "Get inverse of circuit. ALl the gates are newly created."); - - nb::enum_(m, "PauliID") - .value("I", PauliOperator::I) - .value("X", PauliOperator::X) - .value("Y", PauliOperator::Y) - .value("Z", PauliOperator::Z) - .export_values(); - - nb::class_( - m, "PauliOperatorData", "Internal data structure for PauliOperator.") - .def(nb::init(), "coef"_a = 1., "Initialize data with coefficient.") - .def(nb::init(), - "pauli_string"_a, - "coef"_a = 1., - "Initialize data with pauli string.") - .def(nb::init &, - const std::vector &, - Complex>(), - "target_qubit_list"_a, - "pauli_id_list"_a, - "coef"_a = 1., - "Initialize data with target qubits and pauli ids.") - .def(nb::init &, Complex>(), - "pauli_id_par_qubit"_a, - "coef"_a = 1., - "Initialize data with pauli ids per qubit.") - .def(nb::init &, const std::vector &, Complex>(), - "bit_flip_mask"_a, - "phase_flip_mask"_a, - "coef"_a = 1., - "Initialize data with bit flip and phase flip masks.") - .def(nb::init(), - "data"_a, - "Initialize pauli operator from Data object.") - .def("add_single_pauli", - &PauliOperator::Data::add_single_pauli, - "target_qubit"_a, - "pauli_id"_a, - "Add a single pauli operation to the data.") - .def("get_coef", - &PauliOperator::Data::get_coef, - "Get the coefficient of the Pauli operator.") - .def("set_coef", - &PauliOperator::Data::set_coef, - "c"_a, - "Set the coefficient of the Pauli operator.") - .def("get_target_qubit_list", - &PauliOperator::Data::get_target_qubit_list, - "Get the list of target qubits.") - .def("get_pauli_id_list", - &PauliOperator::Data::get_pauli_id_list, - "Get the list of Pauli IDs.") - .def("get_XZ_mask_representation", - &PauliOperator::Data::get_XZ_mask_representation, - "Get the X and Z mask representation as a tuple of vectors."); - - nb::class_( - m, - "PauliOperator", - "Pauli operator as coef and tensor product of single pauli for each qubit.") - .def(nb::init(), "coef"_a = 1., "Initialize operator which just multiplying coef.") - .def(nb::init &, - const std::vector &, - Complex>(), - "target_qubit_list"_a, - "pauli_id_list"_a, - "coef"_a = 1., - "Initialize pauli operator. For each `i`, single pauli correspond to " - "`pauli_id_list[i]` is applied to `target_qubit_list`-th qubit.") - .def(nb::init(), - "pauli_string"_a, - "coef"_a = 1., - "Initialize pauli operator. If `pauli_string` is `\"X0Y2\"`, Pauli-X is applied to " - "0-th qubit and Pauli-Y is applied to 2-th qubit. In `pauli_string`, spaces are " - "ignored.") - .def(nb::init &, Complex>(), - "pauli_id_par_qubit"_a, - "coef"_a = 1., - "Initialize pauli operator. For each `i`, single pauli correspond to " - "`paul_id_per_qubit` is applied to `i`-th qubit.") - .def( - "__init__", - [](PauliOperator *t, - nb::int_ bit_flip_mask_py, - nb::int_ phase_flip_mask_py, - Complex coef) { - internal::BitVector bit_flip_mask(0), phase_flip_mask(0); - const nb::int_ mask(~0ULL); - auto &bit_flip_raw = bit_flip_mask.data_raw(); - assert(bit_flip_raw.empty()); - while (bit_flip_mask_py > nb::int_(0)) { - bit_flip_raw.push_back((std::uint64_t)nb::int_(bit_flip_mask_py & mask)); - bit_flip_mask_py >>= nb::int_(64); - } - auto &phase_flip_raw = phase_flip_mask.data_raw(); - assert(phase_flip_raw.empty()); - while (phase_flip_mask_py > nb::int_(0)) { - phase_flip_raw.push_back((std::uint64_t)nb::int_(phase_flip_mask_py & mask)); - phase_flip_mask_py >>= nb::int_(64); - } - new (t) PauliOperator(bit_flip_mask, phase_flip_mask, coef); - }, - "bit_flip_mask"_a, - "phase_flip_mask"_a, - "coef"_a = 1., - "Initialize pauli operator. For each `i`, single pauli applied to `i`-th qubit is got " - "from `i-th` bit of `bit_flip_mask` and `phase_flip_mask` as follows.\n\n.. " - "csv-table::\n\n \"bit_flip\",\"phase_flip\",\"pauli\"\n \"0\",\"0\",\"I\"\n " - "\"0\",\"1\",\"Z\"\n \"1\",\"0\",\"X\"\n \"1\",\"1\",\"Y\"") - .def("get_coef", &PauliOperator::get_coef, "Get property `coef`.") - .def("get_target_qubit_list", - &PauliOperator::get_target_qubit_list, - "Get qubits to be applied pauli.") - .def("get_pauli_id_list", - &PauliOperator::get_pauli_id_list, - "Get pauli id to be applied. The order is correspond to the result of " - "`get_target_qubit_list`") - .def( - "get_XZ_mask_representation", - [](const PauliOperator &pauli) { - const auto [x_mask, z_mask] = pauli.get_XZ_mask_representation(); - nb::int_ x_mask_py(0); - for (std::uint64_t i = 0; i < x_mask.size(); ++i) { - x_mask_py |= nb::int_(x_mask[i]) << nb::int_(i); - } - nb::int_ z_mask_py(0); - for (std::uint64_t i = 0; i < z_mask.size(); ++i) { - z_mask_py |= nb::int_(z_mask[i]) << nb::int_(i); - } - return std::make_tuple(x_mask_py, z_mask_py); - }, - "Get single-pauli property as binary integer representation. See description of " - "`__init__(bit_flip_mask_py: int, phase_flip_mask_py: int, coef: float=1.)` for " - "details.") - .def("get_pauli_string", - &PauliOperator::get_pauli_string, - "Get single-pauli property as string representation. See description of " - "`__init__(pauli_string: str, coef: float=1.)` for details.") - .def("get_dagger", &PauliOperator::get_dagger, "Get adjoint operator.") - .def("get_qubit_count", - &PauliOperator::get_qubit_count, - "Get num of qubits to applied with, when count from 0-th qubit. Subset of $[0, " - "\\mathrm{qubit_count})$ is the target.") - .def("apply_to_state", &PauliOperator::apply_to_state, "Apply pauli to state vector.") - .def("get_expectation_value", - &PauliOperator::get_expectation_value, - "Get expectation value of measuring state vector. $\\bra{\\psi}P\\ket{\\psi}$.") - .def("get_transition_amplitude", - &PauliOperator::get_transition_amplitude, - "Get transition amplitude of measuring state vector. $\\bra{\\chi}P\\ket{\\psi}$.") - .def(nb::self * nb::self) - .def(nb::self * Complex()); + internal::bind_operator_pauli_operator_hpp(m); + internal::bind_operator_operator_hpp(m); - nb::class_(m, "Operator", "General quantum operator class.") - .def(nb::init(), - "qubit_count"_a, - "Initialize operator with specified number of qubits.") - .def("is_hermitian", &Operator::is_hermitian, "Check if the operator is Hermitian.") - .def("n_qubits", &Operator::n_qubits, "Get the number of qubits the operator acts on.") - .def("terms", &Operator::terms, "Get the list of Pauli terms that make up the operator.") - .def("to_string", &Operator::to_string, "Get string representation of the operator.") - .def("add_operator", - nb::overload_cast(&Operator::add_operator), - "Add a Pauli operator to this operator.") - .def( - "add_random_operator", - [](Operator &op, std::uint64_t operator_count, std::optional seed) { - return op.add_random_operator(operator_count, - seed.value_or(std::random_device{}())); - }, - "operator_count"_a, - "seed"_a = std::nullopt, - "Add a specified number of random Pauli operators to this operator. An optional seed " - "can be provided for reproducibility.") - .def("optimize", &Operator::optimize, "Optimize the operator by combining like terms.") - .def("get_dagger", - &Operator::get_dagger, - "Get the adjoint (Hermitian conjugate) of the operator.") - .def("apply_to_state", &Operator::apply_to_state, "Apply the operator to a state vector.") - .def("get_expectation_value", - &Operator::get_expectation_value, - "Get the expectation value of the operator with respect to a state vector.") - .def("get_transition_amplitude", - &Operator::get_transition_amplitude, - "Get the transition amplitude of the operator between two state vectors.") - .def(nb::self *= Complex()) - .def(nb::self * Complex()) - .def(+nb::self) - .def(-nb::self) - .def(nb::self += nb::self) - .def(nb::self + nb::self) - .def(nb::self -= nb::self) - .def(nb::self - nb::self) - .def(nb::self * nb::self) - .def(nb::self *= nb::self) - .def(nb::self += PauliOperator()) - .def(nb::self + PauliOperator()) - .def(nb::self -= PauliOperator()) - .def(nb::self - PauliOperator()) - .def(nb::self *= PauliOperator()) - .def(nb::self * PauliOperator()); initialize(); std::atexit(&cleanup); } diff --git a/scaluq/CMakeLists.txt b/scaluq/CMakeLists.txt index a7de166d..274fd106 100644 --- a/scaluq/CMakeLists.txt +++ b/scaluq/CMakeLists.txt @@ -5,8 +5,8 @@ target_sources(scaluq PRIVATE gate/update_ops_dense_matrix.cpp gate/update_ops_sparse_matrix.cpp gate/update_ops_standard.cpp - gate/update_ops_pauli.cpp # gate/merge_gate.cpp + operator/apply_pauli.cpp operator/pauli_operator.cpp operator/operator.cpp state/state_vector.cpp diff --git a/scaluq/all.hpp b/scaluq/all.hpp index 4330187f..f79a69a4 100644 --- a/scaluq/all.hpp +++ b/scaluq/all.hpp @@ -8,6 +8,7 @@ #include "gate/param_gate.hpp" #include "gate/param_gate_factory.hpp" #include "gate/update_ops.hpp" +#include "operator/apply_pauli.hpp" #include "operator/operator.hpp" #include "operator/pauli_operator.hpp" #include "state/state_vector.hpp" diff --git a/scaluq/circuit/circuit.cpp b/scaluq/circuit/circuit.cpp index a349e877..4cf3badb 100644 --- a/scaluq/circuit/circuit.cpp +++ b/scaluq/circuit/circuit.cpp @@ -7,11 +7,11 @@ std::uint64_t Circuit::calculate_depth() const { std::vector filled_step(_n_qubits, 0ULL); for (const auto& gate : _gate_list) { std::vector control_qubits = - gate.index() == 0 ? std::get<0>(gate)->get_control_qubit_list() - : std::get<1>(gate).first->get_control_qubit_list(); + gate.index() == 0 ? std::get<0>(gate)->control_qubit_list() + : std::get<1>(gate).first->control_qubit_list(); std::vector target_qubits = - gate.index() == 0 ? std::get<0>(gate)->get_target_qubit_list() - : std::get<1>(gate).first->get_target_qubit_list(); + gate.index() == 0 ? std::get<0>(gate)->target_qubit_list() + : std::get<1>(gate).first->target_qubit_list(); std::uint64_t max_step_amount_target_qubits = 0; for (std::uint64_t control : control_qubits) { if (max_step_amount_target_qubits < filled_step[control]) { @@ -123,8 +123,8 @@ Circuit Circuit::get_inverse() const { } void Circuit::check_gate_is_valid(const Gate& gate) const { - auto targets = gate->get_target_qubit_list(); - auto controls = gate->get_control_qubit_list(); + auto targets = gate->target_qubit_list(); + auto controls = gate->control_qubit_list(); bool valid = true; if (!targets.empty()) valid &= *std::max_element(targets.begin(), targets.end()) < _n_qubits; if (!controls.empty()) valid &= *std::max_element(controls.begin(), controls.end()) < _n_qubits; @@ -134,8 +134,8 @@ void Circuit::check_gate_is_valid(const Gate& gate) const { } void Circuit::check_gate_is_valid(const ParamGate& param_gate) const { - auto targets = param_gate->get_target_qubit_list(); - auto controls = param_gate->get_control_qubit_list(); + auto targets = param_gate->target_qubit_list(); + auto controls = param_gate->control_qubit_list(); bool valid = true; if (!targets.empty()) valid &= *std::max_element(targets.begin(), targets.end()) < _n_qubits; if (!controls.empty()) valid &= *std::max_element(controls.begin(), controls.end()) < _n_qubits; diff --git a/scaluq/circuit/circuit.hpp b/scaluq/circuit/circuit.hpp index 656b87fe..5d173d5a 100644 --- a/scaluq/circuit/circuit.hpp +++ b/scaluq/circuit/circuit.hpp @@ -15,7 +15,7 @@ class Circuit { [[nodiscard]] inline std::uint64_t n_qubits() const { return _n_qubits; } [[nodiscard]] inline const std::vector& gate_list() const { return _gate_list; } - [[nodiscard]] inline std::uint64_t gate_count() { return _gate_list.size(); } + [[nodiscard]] inline std::uint64_t n_gates() { return _gate_list.size(); } [[nodiscard]] inline const std::set key_set() const { std::set key_set; for (auto&& gate : _gate_list) { @@ -23,13 +23,13 @@ class Circuit { } return key_set; } - [[nodiscard]] inline const GateWithKey& get(std::uint64_t idx) const { + [[nodiscard]] inline const GateWithKey& get_gate_at(std::uint64_t idx) const { if (idx >= _gate_list.size()) { - throw std::runtime_error("Circuit::get(std::uint64_t): index out of bounds"); + throw std::runtime_error("Circuit::get_gate_at(std::uint64_t): index out of bounds"); } return _gate_list[idx]; } - [[nodiscard]] inline std::optional get_key(std::uint64_t idx) { + [[nodiscard]] inline std::optional get_param_key_at(std::uint64_t idx) { if (idx >= _gate_list.size()) { throw std::runtime_error( "Circuit::get_parameter_key(std::uint64_t): index out of bounds"); @@ -62,4 +62,58 @@ class Circuit { void check_gate_is_valid(const Gate& gate) const; void check_gate_is_valid(const ParamGate& gate) const; }; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_circuit_circuit_hpp(nb::module_& m) { + nb::class_(m, "Circuit", "Quantum circuit represented as gate array") + .def(nb::init(), "Initialize empty circuit of specified qubits.") + .def("n_qubits", &Circuit::n_qubits, "Get property of `n_qubits`.") + .def("gate_list", + &Circuit::gate_list, + "Get property of `gate_list`.", + nb::rv_policy::reference) + .def("n_gates", &Circuit::n_gates, "Get property of `n_gates`.") + .def("key_set", &Circuit::key_set, "Get set of keys of parameters.") + .def("get_gate_at", &Circuit::get_gate_at, "Get reference of i-th gate.") + .def("get_param_key_at", + &Circuit::get_param_key_at, + "Get parameter key of i-th gate. If it is not parametric, return None.") + .def("calculate_depth", &Circuit::calculate_depth, "Get depth of circuit.") + .def("add_gate", + nb::overload_cast(&Circuit::add_gate), + "Add gate. Given gate is copied.") + .def("add_param_gate", + nb::overload_cast(&Circuit::add_param_gate), + "Add parametric gate with specifing key. Given param_gate is copied.") + .def("add_circuit", + nb::overload_cast(&Circuit::add_circuit), + "Add all gates in specified circuit. Given gates are copied.") + .def("update_quantum_state", + &Circuit::update_quantum_state, + "Apply gate to the StateVector. StateVector in args is directly updated. If the " + "circuit contains parametric gate, you have to give real value of parameter as " + "dict[str, float] in 2nd arg.") + .def( + "update_quantum_state", + [&](const Circuit& circuit, StateVector& state, nb::kwargs kwargs) { + std::map parameters; + for (auto&& [key, param] : kwargs) { + parameters[nb::cast(key)] = nb::cast(param); + } + circuit.update_quantum_state(state, parameters); + }, + "Apply gate to the StateVector. StateVector in args is directly updated. If the " + "circuit contains parametric gate, you have to give real value of parameter as " + "\"name=value\" format in kwargs.") + .def( + "update_quantum_state", + [](const Circuit& circuit, StateVector& state) { circuit.update_quantum_state(state); }) + .def("copy", &Circuit::copy, "Copy circuit. All the gates inside is copied.") + .def("get_inverse", + &Circuit::get_inverse, + "Get inverse of circuit. All the gates are newly created."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/constant.hpp b/scaluq/constant.hpp index 590bd9b6..cc415f5b 100644 --- a/scaluq/constant.hpp +++ b/scaluq/constant.hpp @@ -6,17 +6,17 @@ #include "types.hpp" namespace scaluq { -//! PI value -KOKKOS_INLINE_FUNCTION -double PI() { return 3.141592653589793; } -//! square root of 2 -KOKKOS_INLINE_FUNCTION -double SQRT2() { return 1.4142135623730950; } +namespace internal { + +#define DEF_MATH_CONSTANT(TRAIT, VALUE) \ + template \ + inline constexpr auto TRAIT##_v = std::enable_if_t, T>(VALUE); \ + inline constexpr auto TRAIT = TRAIT##_v //! inverse square root of 2 KOKKOS_INLINE_FUNCTION -double INVERSE_SQRT2() { return 0.707106781186547; } +double INVERSE_SQRT2() { return Kokkos::numbers::sqrt2 / 2; } //! cosine pi/8 KOKKOS_INLINE_FUNCTION @@ -28,71 +28,69 @@ double SINPI8() { return 0.382683432365090; } //! identity matrix KOKKOS_INLINE_FUNCTION -matrix_2_2 I_GATE() { return {1, 0, 0, 1}; } +Matrix2x2 I_GATE() { return {1, 0, 0, 1}; } //! Pauli matrix X KOKKOS_INLINE_FUNCTION -matrix_2_2 X_GATE() { return {0, 1, 1, 0}; } +Matrix2x2 X_GATE() { return {0, 1, 1, 0}; } //! Pauli matrix Y KOKKOS_INLINE_FUNCTION -matrix_2_2 Y_GATE() { return {0, Complex(0, -1), Complex(0, 1), 0}; } +Matrix2x2 Y_GATE() { return {0, Complex(0, -1), Complex(0, 1), 0}; } //! Pauli matrix Z KOKKOS_INLINE_FUNCTION -matrix_2_2 Z_GATE() { return {1, 0, 0, -1}; } +Matrix2x2 Z_GATE() { return {1, 0, 0, -1}; } //! list of Pauli matrix I,X,Y,Z -// std::array PAULI_MATRIX = {I_GATE, X_GATE, Y_GATE, Z_GATE}; +// std::array PAULI_MATRIX = {I_GATE, X_GATE, +// Y_GATE, Z_GATE}; //! S-gate KOKKOS_INLINE_FUNCTION -matrix_2_2 S_GATE_MATRIX() { return {1, 0, 0, Complex(0, 1)}; } +Matrix2x2 S_GATE_MATRIX() { return {1, 0, 0, Complex(0, 1)}; } //! Sdag-gate KOKKOS_INLINE_FUNCTION -matrix_2_2 S_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(0, -1)}; } +Matrix2x2 S_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(0, -1)}; } //! T-gate KOKKOS_INLINE_FUNCTION -matrix_2_2 T_GATE_MATRIX() { - return {COSPI8() - Complex(0, SINPI8()), 0., 0., COSPI8() + Complex(0, SINPI8()) * SINPI8()}; -} +Matrix2x2 T_GATE_MATRIX() { return {1, 0, 0, Complex(INVERSE_SQRT2(), INVERSE_SQRT2())}; } //! Tdag-gate KOKKOS_INLINE_FUNCTION -matrix_2_2 T_DAG_GATE_MATRIX() { - return {COSPI8() + Complex(0, SINPI8()), 0., 0., COSPI8() - Complex(0, SINPI8())}; -} +Matrix2x2 T_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(INVERSE_SQRT2(), -INVERSE_SQRT2())}; } //! Hadamard gate KOKKOS_INLINE_FUNCTION -matrix_2_2 HADAMARD_MATRIX() { +Matrix2x2 HADAMARD_MATRIX() { return {INVERSE_SQRT2(), INVERSE_SQRT2(), INVERSE_SQRT2(), -INVERSE_SQRT2()}; } //! square root of X gate KOKKOS_INLINE_FUNCTION -matrix_2_2 SQRT_X_GATE_MATRIX() { +Matrix2x2 SQRT_X_GATE_MATRIX() { return {Complex(0.5, 0.5), Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(0.5, 0.5)}; } //! square root of Y gate KOKKOS_INLINE_FUNCTION -matrix_2_2 SQRT_Y_GATE_MATRIX() { +Matrix2x2 SQRT_Y_GATE_MATRIX() { return {Complex(0.5, 0.5), Complex(-0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5)}; } //! square root dagger of X gate KOKKOS_INLINE_FUNCTION -matrix_2_2 SQRT_X_DAG_GATE_MATRIX() { +Matrix2x2 SQRT_X_DAG_GATE_MATRIX() { return {Complex(0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5), Complex(0.5, -0.5)}; } //! square root dagger of Y gate KOKKOS_INLINE_FUNCTION -matrix_2_2 SQRT_Y_DAG_GATE_MATRIX() { +Matrix2x2 SQRT_Y_DAG_GATE_MATRIX() { return {Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(-0.5, 0.5), Complex(0.5, -0.5)}; } //! Projection to 0 KOKKOS_INLINE_FUNCTION -matrix_2_2 PROJ_0_MATRIX() { return {1, 0, 0, 0}; } +Matrix2x2 PROJ_0_MATRIX() { return {1, 0, 0, 0}; } //! Projection to 1 KOKKOS_INLINE_FUNCTION -matrix_2_2 PROJ_1_MATRIX() { return {0, 0, 0, 1}; } +Matrix2x2 PROJ_1_MATRIX() { return {0, 0, 0, 1}; } //! complex values for exp(j * i*pi/4 ) KOKKOS_INLINE_FUNCTION -array_4 PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; } +Kokkos::Array PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; } //! complex values for exp(-j * i*pi/4 ) KOKKOS_INLINE_FUNCTION -array_4 PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; } +Kokkos::Array PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; } +} // namespace internal } // namespace scaluq diff --git a/scaluq/gate/gate.hpp b/scaluq/gate/gate.hpp index 9a726d1c..6b8d197c 100644 --- a/scaluq/gate/gate.hpp +++ b/scaluq/gate/gate.hpp @@ -36,13 +36,11 @@ class U1GateImpl; class U2GateImpl; class U3GateImpl; class OneTargetMatrixGateImpl; -class CXGateImpl; -class CZGateImpl; -class CCXGateImpl; class SwapGateImpl; class TwoTargetMatrixGateImpl; class PauliGateImpl; class PauliRotationGateImpl; +class ProbablisticGateImpl; class SparseMatrixGateImpl; class DenseMatrixGateImpl; @@ -76,57 +74,81 @@ enum class GateType { U2, U3, OneTargetMatrix, - CX, - CZ, - CCX, Swap, TwoTargetMatrix, Pauli, PauliRotation, SparseMatrix, - DenseMatrix + DenseMatrix, + Probablistic }; template constexpr GateType get_gate_type() { - if constexpr (std::is_same_v) return GateType::Unknown; - if constexpr (std::is_same_v) return GateType::I; - if constexpr (std::is_same_v) return GateType::GlobalPhase; - if constexpr (std::is_same_v) return GateType::X; - if constexpr (std::is_same_v) return GateType::Y; - if constexpr (std::is_same_v) return GateType::Z; - if constexpr (std::is_same_v) return GateType::H; - if constexpr (std::is_same_v) return GateType::S; - if constexpr (std::is_same_v) return GateType::Sdag; - if constexpr (std::is_same_v) return GateType::T; - if constexpr (std::is_same_v) return GateType::Tdag; - if constexpr (std::is_same_v) return GateType::SqrtX; - if constexpr (std::is_same_v) return GateType::SqrtXdag; - if constexpr (std::is_same_v) return GateType::SqrtY; - if constexpr (std::is_same_v) return GateType::SqrtYdag; - if constexpr (std::is_same_v) return GateType::P0; - if constexpr (std::is_same_v) return GateType::P1; - if constexpr (std::is_same_v) return GateType::RX; - if constexpr (std::is_same_v) return GateType::RY; - if constexpr (std::is_same_v) return GateType::RZ; - if constexpr (std::is_same_v) return GateType::U1; - if constexpr (std::is_same_v) return GateType::U2; - if constexpr (std::is_same_v) return GateType::U3; - if constexpr (std::is_same_v) + using TWithoutConst = std::remove_cv_t; + if constexpr (std::is_same_v) + return GateType::Unknown; + else if constexpr (std::is_same_v) + return GateType::I; + else if constexpr (std::is_same_v) + return GateType::GlobalPhase; + else if constexpr (std::is_same_v) + return GateType::X; + else if constexpr (std::is_same_v) + return GateType::Y; + else if constexpr (std::is_same_v) + return GateType::Z; + else if constexpr (std::is_same_v) + return GateType::H; + else if constexpr (std::is_same_v) + return GateType::S; + else if constexpr (std::is_same_v) + return GateType::Sdag; + else if constexpr (std::is_same_v) + return GateType::T; + else if constexpr (std::is_same_v) + return GateType::Tdag; + else if constexpr (std::is_same_v) + return GateType::SqrtX; + else if constexpr (std::is_same_v) + return GateType::SqrtXdag; + else if constexpr (std::is_same_v) + return GateType::SqrtY; + else if constexpr (std::is_same_v) + return GateType::SqrtYdag; + else if constexpr (std::is_same_v) + return GateType::P0; + else if constexpr (std::is_same_v) + return GateType::P1; + else if constexpr (std::is_same_v) + return GateType::RX; + else if constexpr (std::is_same_v) + return GateType::RY; + else if constexpr (std::is_same_v) + return GateType::RZ; + else if constexpr (std::is_same_v) + return GateType::U1; + else if constexpr (std::is_same_v) + return GateType::U2; + else if constexpr (std::is_same_v) + return GateType::U3; + else if constexpr (std::is_same_v) return GateType::OneTargetMatrix; - if constexpr (std::is_same_v) return GateType::CX; - if constexpr (std::is_same_v) return GateType::CZ; - if constexpr (std::is_same_v) return GateType::CCX; - if constexpr (std::is_same_v) return GateType::Swap; - if constexpr (std::is_same_v) + else if constexpr (std::is_same_v) + return GateType::Swap; + else if constexpr (std::is_same_v) return GateType::TwoTargetMatrix; - if constexpr (std::is_same_v) return GateType::Pauli; - if constexpr (std::is_same_v) + else if constexpr (std::is_same_v) + return GateType::Pauli; + else if constexpr (std::is_same_v) return GateType::PauliRotation; - if constexpr (std::is_same_v) return GateType::SparseMatrix; - if constexpr (std::is_same_v) return GateType::DenseMatrix; - static_assert("unknown GateImpl"); - return GateType::Unknown; + else if constexpr (std::is_same_v) + return GateType::Probablistic; + else if constexpr (std::is_same_v) + return GateType::SparseMatrix; + else if constexpr (std::is_same_v) + return GateType::DenseMatrix; + static_assert(internal::lazy_false_v, "unknown GateImpl"); } namespace internal { @@ -142,37 +164,53 @@ class GateBase : public std::enable_shared_from_this { } } + std::string get_qubit_info_as_string(const std::string& indent) const { + std::ostringstream ss; + auto targets = target_qubit_list(); + auto controls = control_qubit_list(); + ss << indent << " Target Qubits: {"; + for (std::uint32_t i = 0; i < targets.size(); ++i) + ss << targets[i] << (i == targets.size() - 1 ? "" : ", "); + ss << "}\n"; + ss << indent << " Control Qubits: {"; + for (std::uint32_t i = 0; i < controls.size(); ++i) + ss << controls[i] << (i == controls.size() - 1 ? "" : ", "); + ss << "}"; + return ss.str(); + } + public: GateBase(std::uint64_t target_mask, std::uint64_t control_mask) : _target_mask(target_mask), _control_mask(control_mask) { if (_target_mask & _control_mask) [[unlikely]] { throw std::runtime_error( "Error: Gate::Gate(std::uint64_t target_mask, std::uint64_t control_mask) : Target " - "and " - "control qubits must not overlap."); + "and control qubits must not overlap."); } } virtual ~GateBase() = default; - [[nodiscard]] virtual std::vector get_target_qubit_list() const { + [[nodiscard]] virtual std::vector target_qubit_list() const { return mask_to_vector(_target_mask); } - [[nodiscard]] virtual std::vector get_control_qubit_list() const { + [[nodiscard]] virtual std::vector control_qubit_list() const { return mask_to_vector(_control_mask); } - [[nodiscard]] virtual std::vector get_operand_qubit_list() const { + [[nodiscard]] virtual std::vector operand_qubit_list() const { return mask_to_vector(_target_mask | _control_mask); } - [[nodiscard]] virtual std::uint64_t get_target_qubit_mask() const { return _target_mask; } - [[nodiscard]] virtual std::uint64_t get_control_qubit_mask() const { return _control_mask; } - [[nodiscard]] virtual std::uint64_t get_operand_qubit_mask() const { + [[nodiscard]] virtual std::uint64_t target_qubit_mask() const { return _target_mask; } + [[nodiscard]] virtual std::uint64_t control_qubit_mask() const { return _control_mask; } + [[nodiscard]] virtual std::uint64_t operand_qubit_mask() const { return _target_mask | _control_mask; } [[nodiscard]] virtual Gate get_inverse() const = 0; - [[nodiscard]] virtual ComplexMatrix get_matrix() const = 0; + [[nodiscard]] virtual internal::ComplexMatrix get_matrix() const = 0; virtual void update_quantum_state(StateVector& state_vector) const = 0; + + [[nodiscard]] virtual std::string to_string(const std::string& indent = "") const = 0; }; template @@ -232,7 +270,113 @@ class GatePtr { } return _gate_ptr.get(); } + + friend std::ostream& operator<<(std::ostream& os, GatePtr gate) { + os << gate->to_string(); + return os; + } }; } // namespace internal +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +#define DEF_GATE_BASE(GATE_TYPE, DESCRIPTION) \ + nb::class_(m, #GATE_TYPE, DESCRIPTION) \ + .def("gate_type", &GATE_TYPE::gate_type, "Get gate type as `GateType` enum.") \ + .def( \ + "target_qubit_list", \ + [](const GATE_TYPE& gate) { return gate->target_qubit_list(); }, \ + "Get target qubits as `list[int]`. **Control qubits is not included.**") \ + .def( \ + "control_qubit_list", \ + [](const GATE_TYPE& gate) { return gate->control_qubit_list(); }, \ + "Get control qubits as `list[int]`.") \ + .def( \ + "operand_qubit_list", \ + [](const GATE_TYPE& gate) { return gate->operand_qubit_list(); }, \ + "Get target and control qubits as `list[int]`.") \ + .def( \ + "target_qubit_mask", \ + [](const GATE_TYPE& gate) { return gate->target_qubit_mask(); }, \ + "Get target qubits as mask. **Control qubits is not included.**") \ + .def( \ + "control_qubit_mask", \ + [](const GATE_TYPE& gate) { return gate->control_qubit_mask(); }, \ + "Get control qubits as mask.") \ + .def( \ + "operand_qubit_mask", \ + [](const GATE_TYPE& gate) { return gate->operand_qubit_mask(); }, \ + "Get target and control qubits as mask.") \ + .def( \ + "get_inverse", \ + [](const GATE_TYPE& gate) { return gate->get_inverse(); }, \ + "Generate inverse gate as `Gate` type. If not exists, return None.") \ + .def( \ + "update_quantum_state", \ + [](const GATE_TYPE& gate, StateVector& state_vector) { \ + gate->update_quantum_state(state_vector); \ + }, \ + "Apply gate to `state_vector`. `state_vector` in args is directly updated.") \ + .def( \ + "get_matrix", \ + [](const GATE_TYPE& gate) { return gate->get_matrix(); }, \ + "Get matrix representation of the gate.") \ + .def( \ + "to_string", \ + [](const GATE_TYPE& gate) { return gate->to_string(""); }, \ + "Get string representation of the gate.") \ + .def( \ + "__str__", \ + [](const GATE_TYPE& gate) { return gate->to_string(""); }, \ + "Get string representation of the gate.") + +nb::class_ gate_base_def; + +#define DEF_GATE(GATE_TYPE, DESCRIPTION) \ + ::scaluq::internal::gate_base_def.def(nb::init(), "Upcast from `" #GATE_TYPE "`."); \ + DEF_GATE_BASE( \ + GATE_TYPE, \ + DESCRIPTION \ + "\n\n.. note:: Upcast is required to use gate-general functions (ex: add to Circuit).") \ + .def(nb::init()) + +void bind_gate_gate_hpp(nb::module_& m) { + nb::enum_(m, "GateType", "Enum of Gate Type.") + .value("I", GateType::I) + .value("GlobalPhase", GateType::GlobalPhase) + .value("X", GateType::X) + .value("Y", GateType::Y) + .value("Z", GateType::Z) + .value("H", GateType::H) + .value("S", GateType::S) + .value("Sdag", GateType::Sdag) + .value("T", GateType::T) + .value("Tdag", GateType::Tdag) + .value("SqrtX", GateType::SqrtX) + .value("SqrtXdag", GateType::SqrtXdag) + .value("SqrtY", GateType::SqrtY) + .value("SqrtYdag", GateType::SqrtYdag) + .value("P0", GateType::P0) + .value("P1", GateType::P1) + .value("RX", GateType::RX) + .value("RY", GateType::RY) + .value("RZ", GateType::RZ) + .value("U1", GateType::U1) + .value("U2", GateType::U2) + .value("U3", GateType::U3) + .value("OneTargetMatrix", GateType::OneTargetMatrix) + .value("Swap", GateType::Swap) + .value("TwoTargetMatrix", GateType::TwoTargetMatrix) + .value("Pauli", GateType::Pauli) + .value("PauliRotation", GateType::PauliRotation); + + gate_base_def = + DEF_GATE_BASE(Gate, + "General class of QuantumGate.\n\n.. note:: Downcast to requred to use " + "gate-specific functions.") + .def(nb::init(), "Just copy shallowly."); +} +} // namespace internal +#endif + } // namespace scaluq diff --git a/scaluq/gate/gate_factory.hpp b/scaluq/gate/gate_factory.hpp index 6ce749b0..1c1b4275 100644 --- a/scaluq/gate/gate_factory.hpp +++ b/scaluq/gate/gate_factory.hpp @@ -150,9 +150,8 @@ inline Gate TwoTargetMatrix(std::uint64_t target1, return internal::GateFactory::create_gate( internal::vector_to_mask({target1, target2}), internal::vector_to_mask(controls), matrix); } -// まだ inline Gate Pauli(const PauliOperator& pauli, const std::vector& controls = {}) { - auto tar = pauli.get_target_qubit_list(); + auto tar = pauli.target_qubit_list(); return internal::GateFactory::create_gate( internal::vector_to_mask(controls), pauli); } @@ -171,10 +170,8 @@ inline Gate DenseMatrix(const std::vector& targets, if (static_cast(matrix.rows()) != dim || static_cast(matrix.cols()) != dim) { throw std::runtime_error( - "gate::DenseMatrix(const std::vector&, const ComplexMatrix&): matrix " - "size " - "must " - "be 2^{n_qubits} x 2^{n_qubits}."); + "gate::DenseMatrix(const std::vector&, const internal::ComplexMatrix&): " + "matrix size must be 2^{n_qubits} x 2^{n_qubits}."); } if (targets.size() == 0) return I(); return internal::GateFactory::create_gate( @@ -192,4 +189,192 @@ inline Gate Probablistic(const std::vector& distribution, gate_list); } } // namespace gate + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_gate_factory_hpp(nb::module_& mgate) { + mgate.def("I", &gate::I, "Generate general Gate class instance of I."); + mgate.def("GlobalPhase", + &gate::GlobalPhase, + "Generate general Gate class instance of GlobalPhase.", + "phase"_a, + "controls"_a = std::vector{}); + mgate.def("X", + &gate::X, + "Generate general Gate class instance of X.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("Y", + &gate::Y, + "Generate general Gate class instance of Y.", + "taget"_a, + "controls"_a = std::vector{}); + mgate.def("Z", + &gate::Z, + "Generate general Gate class instance of Z.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("H", + &gate::H, + "Generate general Gate class instance of H.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("S", + &gate::S, + "Generate general Gate class instance of S.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("Sdag", + &gate::Sdag, + "Generate general Gate class instance of Sdag.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("T", + &gate::T, + "Generate general Gate class instance of T.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("Tdag", + &gate::Tdag, + "Generate general Gate class instance of Tdag.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("SqrtX", + &gate::SqrtX, + "Generate general Gate class instance of SqrtX.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("SqrtXdag", + &gate::SqrtXdag, + "Generate general Gate class instance of SqrtXdag.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("SqrtY", + &gate::SqrtY, + "Generate general Gate class instance of SqrtY.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("SqrtYdag", + &gate::SqrtYdag, + "Generate general Gate class instance of SqrtYdag.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("P0", + &gate::P0, + "Generate general Gate class instance of P0.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("P1", + &gate::P1, + "Generate general Gate class instance of P1.", + "target"_a, + "controls"_a = std::vector{}); + mgate.def("RX", + &gate::RX, + "Generate general Gate class instance of RX.", + "target"_a, + "angle"_a, + "controls"_a = std::vector{}); + mgate.def("RY", + &gate::RY, + "Generate general Gate class instance of RY.", + "target"_a, + "angle"_a, + "controls"_a = std::vector{}); + mgate.def("RZ", + &gate::RZ, + "Generate general Gate class instance of RZ.", + "target"_a, + "angle"_a, + "controls"_a = std::vector{}); + mgate.def("U1", + &gate::U1, + "Generate general Gate class instance of U1.", + "target"_a, + "lambda_"_a, + "controls"_a = std::vector{}); + mgate.def("U2", + &gate::U2, + "Generate general Gate class instance of U2.", + "target"_a, + "phi"_a, + "lambda_"_a, + "controls"_a = std::vector{}); + mgate.def("U3", + &gate::U3, + "Generate general Gate class instance of U3.", + "target"_a, + "theta"_a, + "phi"_a, + "lambda_"_a, + "controls"_a = std::vector{}); + mgate.def("Swap", + &gate::Swap, + "Generate general Gate class instance of Swap.", + "target1"_a, + "target2"_a, + "controls"_a = std::vector{}); + mgate.def( + "CX", + &gate::CX, + "Generate general Gate class instance of CX.\n\n.. note:: CX is a specialization of X."); + mgate.def("CNot", + &gate::CX, + "Generate general Gate class instance of CNot.\n\n.. note:: CNot is an alias of CX."); + mgate.def( + "CZ", + &gate::CZ, + "Generate general Gate class instance of CZ.\n\n.. note:: CZ is a specialization of Z."); + mgate.def( + "CCX", + &gate::CCX, + "Generate general Gate class instance of CXX.\n\n.. note:: CX is a specialization of X."); + mgate.def( + "CCNot", + &gate::CCX, + "Generate general Gate class instance of CCNot.\n\n.. note:: CCNot is an alias of CCX."); + mgate.def("Toffoli", + &gate::CCX, + "Generate general Gate class instance of Toffoli.\n\n.. note:: Toffoli is an alias " + "of CCX."); + mgate.def("OneTargetMatrix", + &gate::OneTargetMatrix, + "Generate general Gate class instance of OneTargetMatrix.", + "target"_a, + "matrix"_a, + "controls"_a = std::vector{}); + mgate.def("TwoTargetMatrix", + &gate::TwoTargetMatrix, + "Generate general Gate class instance of TwoTargetMatrix.", + "target1"_a, + "target2"_a, + "matrix"_a, + "controls"_a = std::vector{}); + mgate.def("DenseMatrix", + &gate::DenseMatrix, + "Generate general Gate class instance of DenseMatrix. IGate, OneTargetMatrixGate or " + "TwoTargetMatrixGate correspond to len(target) is created. The case len(target) >= 3 " + "is currently not supported.", + "targets"_a, + "matrix"_a, + "controls"_a = std::vector{}); + mgate.def("Pauli", + &gate::Pauli, + "Generate general Gate class instance of Pauli.", + "pauli"_a, + "controls"_a = std::vector{}); + mgate.def("PauliRotation", + &gate::PauliRotation, + "Generate general Gate class instance of PauliRotation.", + "pauli"_a, + "angle"_a, + "controls"_a = std::vector{}); + mgate.def("Probablistic", + &gate::Probablistic, + "Generate general Gate class instance of Probablistic.", + "distribution"_a, + "gate_list"_a); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/gate_matrix.hpp b/scaluq/gate/gate_matrix.hpp index e9468325..77548ac3 100644 --- a/scaluq/gate/gate_matrix.hpp +++ b/scaluq/gate/gate_matrix.hpp @@ -13,36 +13,35 @@ namespace scaluq { namespace internal { class OneTargetMatrixGateImpl : public GateBase { - matrix_2_2 _matrix; + Matrix2x2 _matrix; public: OneTargetMatrixGateImpl(std::uint64_t target_mask, std::uint64_t control_mask, const std::array, 2>& matrix) : GateBase(target_mask, control_mask) { - _matrix.val[0][0] = matrix[0][0]; - _matrix.val[0][1] = matrix[0][1]; - _matrix.val[1][0] = matrix[1][0]; - _matrix.val[1][1] = matrix[1][1]; + _matrix[0][0] = matrix[0][0]; + _matrix[0][1] = matrix[0][1]; + _matrix[1][0] = matrix[1][0]; + _matrix[1][1] = matrix[1][1]; } std::array, 2> matrix() const { - return {_matrix.val[0][0], _matrix.val[0][1], _matrix.val[1][0], _matrix.val[1][1]}; + return {_matrix[0][0], _matrix[0][1], _matrix[1][0], _matrix[1][1]}; } Gate get_inverse() const override { return std::make_shared( _target_mask, _control_mask, - std::array, 2>{Kokkos::conj(_matrix.val[0][0]), - Kokkos::conj(_matrix.val[1][0]), - Kokkos::conj(_matrix.val[0][1]), - Kokkos::conj(_matrix.val[1][1])}); - } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); - mat << this->_matrix.val[0][0], this->_matrix.val[0][1], this->_matrix.val[1][0], - this->_matrix.val[1][1]; + std::array, 2>{Kokkos::conj(_matrix[0][0]), + Kokkos::conj(_matrix[1][0]), + Kokkos::conj(_matrix[0][1]), + Kokkos::conj(_matrix[1][1])}); + } + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); + mat << this->_matrix[0][0], this->_matrix[0][1], this->_matrix[1][0], this->_matrix[1][1]; return mat; } @@ -50,10 +49,17 @@ class OneTargetMatrixGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); one_target_dense_matrix_gate(_target_mask, _control_mask, _matrix, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: OneTargetMatrix\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class TwoTargetMatrixGateImpl : public GateBase { - matrix_4_4 _matrix; + Matrix4x4 _matrix; public: TwoTargetMatrixGateImpl(std::uint64_t target_mask, @@ -62,7 +68,7 @@ class TwoTargetMatrixGateImpl : public GateBase { : GateBase(target_mask, control_mask) { for (std::uint64_t i : std::views::iota(0, 4)) { for (std::uint64_t j : std::views::iota(0, 4)) { - _matrix.val[i][j] = matrix[i][j]; + _matrix[i][j] = matrix[i][j]; } } } @@ -71,7 +77,7 @@ class TwoTargetMatrixGateImpl : public GateBase { std::array, 4> matrix; for (std::uint64_t i : std::views::iota(0, 4)) { for (std::uint64_t j : std::views::iota(0, 4)) { - matrix[i][j] = _matrix.val[i][j]; + matrix[i][j] = _matrix[i][j]; } } return matrix; @@ -81,20 +87,18 @@ class TwoTargetMatrixGateImpl : public GateBase { std::array, 4> matrix_dag; for (std::uint64_t i : std::views::iota(0, 4)) { for (std::uint64_t j : std::views::iota(0, 4)) { - matrix_dag[i][j] = Kokkos::conj(_matrix.val[j][i]); + matrix_dag[i][j] = Kokkos::conj(_matrix[j][i]); } } return std::make_shared( _target_mask, _control_mask, matrix_dag); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(4, 4); - mat << this->_matrix.val[0][0], this->_matrix.val[0][1], this->_matrix.val[0][2], - this->_matrix.val[0][3], this->_matrix.val[1][0], this->_matrix.val[1][1], - this->_matrix.val[1][2], this->_matrix.val[1][3], this->_matrix.val[2][0], - this->_matrix.val[2][1], this->_matrix.val[2][2], this->_matrix.val[2][3], - this->_matrix.val[3][0], this->_matrix.val[3][1], this->_matrix.val[3][2], - this->_matrix.val[3][3]; + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(4, 4); + mat << this->_matrix[0][0], this->_matrix[0][1], this->_matrix[0][2], this->_matrix[0][3], + this->_matrix[1][0], this->_matrix[1][1], this->_matrix[1][2], this->_matrix[1][3], + this->_matrix[2][0], this->_matrix[2][1], this->_matrix[2][2], this->_matrix[2][3], + this->_matrix[3][0], this->_matrix[3][1], this->_matrix[3][2], this->_matrix[3][3]; return mat; } @@ -102,6 +106,13 @@ class TwoTargetMatrixGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); two_target_dense_matrix_gate(_target_mask, _control_mask, _matrix, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: TwoTargetMatrix\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class SparseMatrixGateImpl : public GateBase { @@ -183,4 +194,15 @@ using OneTargetMatrixGate = internal::GatePtr using TwoTargetMatrixGate = internal::GatePtr; using SparseMatrixGate = internal::GatePtr; using DenseMatrixGate = internal::GatePtr; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_gate_matrix_hpp(nb::module_& m) { + DEF_GATE(OneTargetMatrixGate, "Specific class of one-qubit dense matrix gate.") + .def("matrix", [](const OneTargetMatrixGate& gate) { return gate->matrix(); }); + DEF_GATE(TwoTargetMatrixGate, "Specific class of two-qubit dense matrix gate.") + .def("matrix", [](const TwoTargetMatrixGate& gate) { return gate->matrix(); }); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/gate_pauli.hpp b/scaluq/gate/gate_pauli.hpp index ee864ae7..59767412 100644 --- a/scaluq/gate/gate_pauli.hpp +++ b/scaluq/gate/gate_pauli.hpp @@ -2,6 +2,7 @@ #include +#include "../operator/apply_pauli.hpp" #include "../operator/pauli_operator.hpp" #include "../util/utility.hpp" #include "gate.hpp" @@ -13,17 +14,29 @@ class PauliGateImpl : public GateBase { public: PauliGateImpl(std::uint64_t control_mask, const PauliOperator& pauli) - : GateBase(vector_to_mask(pauli.get_target_qubit_list()), control_mask), - _pauli(pauli) {} + : GateBase(vector_to_mask(pauli.target_qubit_list()), control_mask), _pauli(pauli) {} PauliOperator pauli() const { return _pauli; }; - std::vector get_pauli_id_list() const { return _pauli.get_pauli_id_list(); } + std::vector pauli_id_list() const { return _pauli.pauli_id_list(); } Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { return this->_pauli.get_matrix(); } + internal::ComplexMatrix get_matrix() const override { return this->_pauli.get_matrix(); } void update_quantum_state(StateVector& state_vector) const override { - pauli_gate(_control_mask, _pauli, state_vector); + auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation(); + apply_pauli(_control_mask, bit_flip_mask, phase_flip_mask, _pauli.coef(), state_vector); + } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + auto controls = control_qubit_list(); + ss << indent << "Gate Type: Pauli\n"; + ss << indent << " Control Qubits: {"; + for (std::uint32_t i = 0; i < controls.size(); ++i) + ss << controls[i] << (i == controls.size() - 1 ? "" : ", "); + ss << "}\n"; + ss << indent << " Pauli Operator: \"" << _pauli.get_pauli_string() << "\""; + return ss.str(); } }; @@ -33,33 +46,62 @@ class PauliRotationGateImpl : public GateBase { public: PauliRotationGateImpl(std::uint64_t control_mask, const PauliOperator& pauli, double angle) - : GateBase(vector_to_mask(pauli.get_target_qubit_list()), control_mask), + : GateBase(vector_to_mask(pauli.target_qubit_list()), control_mask), _pauli(pauli), _angle(angle) {} PauliOperator pauli() const { return _pauli; } - std::vector get_pauli_id_list() const { return _pauli.get_pauli_id_list(); } + std::vector pauli_id_list() const { return _pauli.pauli_id_list(); } double angle() const { return _angle; } Gate get_inverse() const override { return std::make_shared(_control_mask, _pauli, -_angle); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat = this->_pauli.get_matrix_ignoring_coef(); - Complex true_angle = _angle * _pauli.get_coef(); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat = this->_pauli.get_matrix_ignoring_coef(); + Complex true_angle = _angle * _pauli.coef(); StdComplex imag_unit(0, 1); mat = (StdComplex)Kokkos::cos(-true_angle / 2) * - ComplexMatrix::Identity(mat.rows(), mat.cols()) + + internal::ComplexMatrix::Identity(mat.rows(), mat.cols()) + imag_unit * (StdComplex)Kokkos::sin(-true_angle / 2) * mat; return mat; } void update_quantum_state(StateVector& state_vector) const override { - pauli_rotation_gate(_control_mask, _pauli, _angle, state_vector); + auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation(); + apply_pauli_rotation( + _control_mask, bit_flip_mask, phase_flip_mask, _pauli.coef(), _angle, state_vector); + } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + auto controls = control_qubit_list(); + ss << indent << "Gate Type: PauliRotation\n"; + ss << indent << " Angle: " << _angle << "\n"; + ss << indent << " Control Qubits: {"; + for (std::uint32_t i = 0; i < controls.size(); ++i) + ss << controls[i] << (i == controls.size() - 1 ? "" : ", "); + ss << "}\n"; + ss << indent << " Pauli Operator: \"" << _pauli.get_pauli_string() << "\""; + return ss.str(); } }; } // namespace internal using PauliGate = internal::GatePtr; using PauliRotationGate = internal::GatePtr; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_gate_pauli_hpp(nb::module_& m) { + DEF_GATE(PauliGate, + "Specific class of multi-qubit pauli gate, which applies single-qubit Pauli " + "gate to " + "each of qubit."); + DEF_GATE(PauliRotationGate, + "Specific class of multi-qubit pauli-rotation gate, represented as " + "$e^{-i\\frac{\\mathrm{angle}}{2}P}$."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/gate_probablistic.hpp b/scaluq/gate/gate_probablistic.hpp index 9be3946c..bf0a826d 100644 --- a/scaluq/gate/gate_probablistic.hpp +++ b/scaluq/gate/gate_probablistic.hpp @@ -31,34 +31,34 @@ class ProbablisticGateImpl : public GateBase { const std::vector& gate_list() const { return _gate_list; } const std::vector& distribution() const { return _distribution; } - std::vector get_target_qubit_list() const override { + std::vector target_qubit_list() const override { throw std::runtime_error( - "ProbablisticGateImpl::get_target_qubit_list(): This function must not be used in " + "ProbablisticGateImpl::target_qubit_list(): This function must not be used in " "ProbablisticGateImpl."); } - std::vector get_control_qubit_list() const override { + std::vector control_qubit_list() const override { throw std::runtime_error( - "ProbablisticGateImpl::get_control_qubit_list(): This function must not be used in " + "ProbablisticGateImpl::control_qubit_list(): This function must not be used in " "ProbablisticGateImpl."); } - std::vector get_operand_qubit_list() const override { + std::vector operand_qubit_list() const override { throw std::runtime_error( - "ProbablisticGateImpl::get_operand_qubit_list(): This function must not be used in " + "ProbablisticGateImpl::operand_qubit_list(): This function must not be used in " "ProbablisticGateImpl."); } - std::uint64_t get_target_qubit_mask() const override { + std::uint64_t target_qubit_mask() const override { throw std::runtime_error( - "ProbablisticGateImpl::get_target_qubit_mask(): This function must not be used in " + "ProbablisticGateImpl::target_qubit_mask(): This function must not be used in " "ProbablisticGateImpl."); } - std::uint64_t get_control_qubit_mask() const override { + std::uint64_t control_qubit_mask() const override { throw std::runtime_error( - "ProbablisticGateImpl::get_control_qubit_mask(): This function must not be used in " + "ProbablisticGateImpl::control_qubit_mask(): This function must not be used in " "ProbablisticGateImpl."); } - std::uint64_t get_operand_qubit_mask() const override { + std::uint64_t operand_qubit_mask() const override { throw std::runtime_error( - "ProbablisticGateImpl::get_operand_qubit_mask(): This function must not be used in " + "ProbablisticGateImpl::operand_qubit_mask(): This function must not be used in " "ProbablisticGateImpl."); } @@ -70,7 +70,7 @@ class ProbablisticGateImpl : public GateBase { }); return std::make_shared(_distribution, inv_gate_list); } - ComplexMatrix get_matrix() const override { + internal::ComplexMatrix get_matrix() const override { throw std::runtime_error( "ProbablisticGateImpl::get_matrix(): This function must not be used in " "ProbablisticGateImpl."); @@ -85,8 +85,38 @@ class ProbablisticGateImpl : public GateBase { if (i >= _gate_list.size()) i = _gate_list.size() - 1; _gate_list[i]->update_quantum_state(state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + const auto dist = distribution(); + ss << indent << "Gate Type: Probablistic\n"; + for (std::size_t i = 0; i < dist.size(); ++i) { + ss << indent << " --------------------\n"; + ss << indent << " Probability: " << dist[i] << "\n"; + ss << gate_list()[i]->to_string(indent + " ") << (i == dist.size() - 1 ? "" : "\n"); + } + return ss.str(); + } }; } // namespace internal using ProbablisticGate = internal::GatePtr; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_gate_probablistic(nb::module_& m) { + DEF_GATE(ProbablisticGate, + "Specific class of probablistic gate. The gate to apply is picked from a cirtain " + "distribution.") + .def( + "gate_list", + [](const ProbablisticGate& gate) { return gate->gate_list(); }, + nb::rv_policy::reference) + .def( + "distribution", + [](const ProbablisticGate& gate) { return gate->distribution(); }, + nb::rv_policy::reference); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/gate_standard.hpp b/scaluq/gate/gate_standard.hpp index e1046155..c98c4828 100644 --- a/scaluq/gate/gate_standard.hpp +++ b/scaluq/gate/gate_standard.hpp @@ -11,11 +11,20 @@ class IGateImpl : public GateBase { IGateImpl() : GateBase(0, 0) {} Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { return ComplexMatrix::Identity(1, 1); } + internal::ComplexMatrix get_matrix() const override { + return internal::ComplexMatrix::Identity(1, 1); + } void update_quantum_state(StateVector& state_vector) const override { i_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: I\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class GlobalPhaseGateImpl : public GateBase { @@ -31,14 +40,22 @@ class GlobalPhaseGateImpl : public GateBase { Gate get_inverse() const override { return std::make_shared(_control_mask, -_phase); } - ComplexMatrix get_matrix() const override { - return ComplexMatrix::Identity(1, 1) * std::exp(std::complex(0, _phase)); + internal::ComplexMatrix get_matrix() const override { + return internal::ComplexMatrix::Identity(1, 1) * std::exp(std::complex(0, _phase)); } void update_quantum_state(StateVector& state_vector) const override { check_qubit_mask_within_bounds(state_vector); global_phase_gate(_target_mask, _control_mask, _phase, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: GlobalPhase\n"; + ss << indent << " Phase: " << _phase << "\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class RotationGateBase : public GateBase { @@ -57,8 +74,8 @@ class XGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0, 1, 1, 0; return mat; } @@ -67,6 +84,13 @@ class XGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); x_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: X\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class YGateImpl : public GateBase { @@ -74,8 +98,8 @@ class YGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0, -1i, 1i, 0; return mat; } @@ -84,6 +108,13 @@ class YGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); y_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: Y\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class ZGateImpl : public GateBase { @@ -91,8 +122,8 @@ class ZGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, -1; return mat; } @@ -101,6 +132,13 @@ class ZGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); z_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: Z\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class HGateImpl : public GateBase { @@ -108,8 +146,8 @@ class HGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 1, 1, -1; mat /= std::sqrt(2); return mat; @@ -119,6 +157,13 @@ class HGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); h_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: H\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class SGateImpl; @@ -135,8 +180,8 @@ class SGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override; - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, 1i; return mat; } @@ -145,6 +190,13 @@ class SGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); s_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: S\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class SdagGateImpl : public GateBase { @@ -154,8 +206,8 @@ class SdagGateImpl : public GateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, -1i; return mat; } @@ -164,6 +216,13 @@ class SdagGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); sdag_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: Sdag\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; // for resolving dependency issues inline Gate SGateImpl::get_inverse() const { @@ -175,8 +234,8 @@ class TGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override; - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, (1. + 1i) / std::sqrt(2); return mat; } @@ -185,6 +244,13 @@ class TGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); t_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: T\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class TdagGateImpl : public GateBase { @@ -194,8 +260,8 @@ class TdagGateImpl : public GateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, (1. - 1.i) / std::sqrt(2); return mat; } @@ -204,6 +270,13 @@ class TdagGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); tdag_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: Tdag\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; // for resolving dependency issues inline Gate TGateImpl::get_inverse() const { @@ -215,8 +288,8 @@ class SqrtXGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override; - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0.5 + 0.5i, 0.5 - 0.5i, 0.5 - 0.5i, 0.5 + 0.5i; return mat; } @@ -225,6 +298,13 @@ class SqrtXGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); sqrtx_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: SqrtX\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class SqrtXdagGateImpl : public GateBase { @@ -234,8 +314,8 @@ class SqrtXdagGateImpl : public GateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0.5 - 0.5i, 0.5 + 0.5i, 0.5 + 0.5i, 0.5 - 0.5i; return mat; } @@ -244,6 +324,13 @@ class SqrtXdagGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); sqrtxdag_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: SqrtXdag\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; // for resolving dependency issues inline Gate SqrtXGateImpl::get_inverse() const { @@ -255,8 +342,8 @@ class SqrtYGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override; - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0.5 + 0.5i, -0.5 - 0.5i, 0.5 + 0.5i, 0.5 + 0.5i; return mat; } @@ -265,6 +352,13 @@ class SqrtYGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); sqrty_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: SqrtY\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class SqrtYdagGateImpl : public GateBase { @@ -274,8 +368,8 @@ class SqrtYdagGateImpl : public GateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0.5 - 0.5i, 0.5 - 0.5i, -0.5 + 0.5i, 0.5 - 0.5i; return mat; } @@ -284,6 +378,13 @@ class SqrtYdagGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); sqrtydag_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: SqrtYdag\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; // for resolving dependency issues inline Gate SqrtYGateImpl::get_inverse() const { @@ -297,8 +398,8 @@ class P0GateImpl : public GateBase { Gate get_inverse() const override { throw std::runtime_error("P0::get_inverse: Projection gate doesn't have inverse gate"); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, 0; return mat; } @@ -307,6 +408,13 @@ class P0GateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); p0_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: P0\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class P1GateImpl : public GateBase { @@ -316,8 +424,8 @@ class P1GateImpl : public GateBase { Gate get_inverse() const override { throw std::runtime_error("P1::get_inverse: Projection gate doesn't have inverse gate"); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 0, 0, 0, 1; return mat; } @@ -326,6 +434,13 @@ class P1GateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); p1_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: P1\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class RXGateImpl : public RotationGateBase { @@ -335,8 +450,8 @@ class RXGateImpl : public RotationGateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask, -_angle); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << std::cos(_angle / 2), -1i * std::sin(_angle / 2), -1i * std::sin(_angle / 2), std::cos(_angle / 2); return mat; @@ -346,6 +461,14 @@ class RXGateImpl : public RotationGateBase { check_qubit_mask_within_bounds(state_vector); rx_gate(_target_mask, _control_mask, _angle, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: RX\n"; + ss << indent << " Angle: " << this->_angle << "\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class RYGateImpl : public RotationGateBase { @@ -355,8 +478,8 @@ class RYGateImpl : public RotationGateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask, -_angle); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << std::cos(_angle / 2), -std::sin(_angle / 2), std::sin(_angle / 2), std::cos(_angle / 2); return mat; @@ -366,6 +489,14 @@ class RYGateImpl : public RotationGateBase { check_qubit_mask_within_bounds(state_vector); ry_gate(_target_mask, _control_mask, _angle, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: RY\n"; + ss << indent << " Angle: " << this->_angle << "\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class RZGateImpl : public RotationGateBase { @@ -375,8 +506,8 @@ class RZGateImpl : public RotationGateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask, -_angle); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << std::exp(-0.5i * _angle), 0, 0, std::exp(0.5i * _angle); return mat; } @@ -385,6 +516,14 @@ class RZGateImpl : public RotationGateBase { check_qubit_mask_within_bounds(state_vector); rz_gate(_target_mask, _control_mask, _angle, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: RZ\n"; + ss << indent << " Angle: " << this->_angle << "\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class U1GateImpl : public GateBase { @@ -399,8 +538,8 @@ class U1GateImpl : public GateBase { Gate get_inverse() const override { return std::make_shared(_target_mask, _control_mask, -_lambda); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << 1, 0, 0, std::exp(1i * _lambda); return mat; } @@ -409,6 +548,13 @@ class U1GateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); u1_gate(_target_mask, _control_mask, _lambda, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: U1\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class U2GateImpl : public GateBase { double _phi, _lambda; @@ -421,14 +567,17 @@ class U2GateImpl : public GateBase { double lambda() const { return _lambda; } Gate get_inverse() const override { - return std::make_shared( - _target_mask, _control_mask, -_lambda - PI(), -_phi + PI()); - } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); - mat << std::cos(PI() / 4.), -std::exp(1i * _lambda) * std::sin(PI() / 4.), - std::exp(1i * _phi) * std::sin(PI() / 4.), - std::exp(1i * _phi) * std::exp(1i * _lambda) * std::cos(PI() / 4.); + return std::make_shared(_target_mask, + _control_mask, + -_lambda - Kokkos::numbers::pi, + -_phi + Kokkos::numbers::pi); + } + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); + mat << std::cos(Kokkos::numbers::pi / 4.), + -std::exp(1i * _lambda) * std::sin(Kokkos::numbers::pi / 4.), + std::exp(1i * _phi) * std::sin(Kokkos::numbers::pi / 4.), + std::exp(1i * _phi) * std::exp(1i * _lambda) * std::cos(Kokkos::numbers::pi / 4.); return mat; } @@ -436,6 +585,13 @@ class U2GateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); u2_gate(_target_mask, _control_mask, _phi, _lambda, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: U2\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class U3GateImpl : public GateBase { @@ -457,8 +613,8 @@ class U3GateImpl : public GateBase { return std::make_shared( _target_mask, _control_mask, -_theta, -_lambda, -_phi); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat(2, 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat(2, 2); mat << std::cos(_theta / 2.), -std::exp(1i * _lambda) * std::sin(_theta / 2.), std::exp(1i * _phi) * std::sin(_theta / 2.), std::exp(1i * _phi) * std::exp(1i * _lambda) * std::cos(_theta / 2.); @@ -469,6 +625,13 @@ class U3GateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); u3_gate(_target_mask, _control_mask, _theta, _phi, _lambda, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: U3\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; class SwapGateImpl : public GateBase { @@ -476,8 +639,8 @@ class SwapGateImpl : public GateBase { using GateBase::GateBase; Gate get_inverse() const override { return shared_from_this(); } - ComplexMatrix get_matrix() const override { - ComplexMatrix mat = ComplexMatrix::Identity(1 << 2, 1 << 2); + internal::ComplexMatrix get_matrix() const override { + internal::ComplexMatrix mat = internal::ComplexMatrix::Identity(1 << 2, 1 << 2); mat << 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1; return mat; } @@ -486,6 +649,13 @@ class SwapGateImpl : public GateBase { check_qubit_mask_within_bounds(state_vector); swap_gate(_target_mask, _control_mask, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: Swap\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; } // namespace internal @@ -515,4 +685,90 @@ using U2Gate = internal::GatePtr; using U3Gate = internal::GatePtr; using SwapGate = internal::GatePtr; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_gate_standard_hpp(nb::module_& m) { + DEF_GATE(IGate, "Specific class of Pauli-I gate."); + DEF_GATE(GlobalPhaseGate, + "Specific class of gate, which rotate global phase, represented as " + "$e^{i\\mathrm{phase}}I$.") + .def( + "phase", + [](const GlobalPhaseGate& gate) { return gate->phase(); }, + "Get `phase` property"); + DEF_GATE(XGate, "Specific class of Pauli-X gate."); + DEF_GATE(YGate, "Specific class of Pauli-Y gate."); + DEF_GATE(ZGate, "Specific class of Pauli-Z gate."); + DEF_GATE(HGate, "Specific class of Hadamard gate."); + DEF_GATE(SGate, + "Specific class of S gate, represented as $\\begin{bmatrix}\n1 & 0\\\\\n0 & " + "i\n\\end{bmatrix}$."); + DEF_GATE(SdagGate, "Specific class of inverse of S gate."); + DEF_GATE(TGate, + "Specific class of T gate, represented as $\\begin{bmatrix}\n1 & 0\\\\\n0 & " + "e^{i\\pi/4}\n\\end{bmatrix}$."); + DEF_GATE(TdagGate, "Specific class of inverse of T gate."); + DEF_GATE(SqrtXGate, + "Specific class of sqrt(X) gate, represented as $\\begin{bmatrix}\n1+i & 1-i\\\\\n1-i " + "& 1+i\n\\end{bmatrix}$."); + DEF_GATE(SqrtXdagGate, "Specific class of inverse of sqrt(X) gate."); + DEF_GATE(SqrtYGate, + "Specific class of sqrt(Y) gate, represented as $\\begin{bmatrix}\n1+i & -1-i " + "\\\\\n1+i & 1+i\n\\end{bmatrix}$."); + DEF_GATE(SqrtYdagGate, "Specific class of inverse of sqrt(Y) gate."); + DEF_GATE( + P0Gate, + "Specific class of projection gate to $\\ket{0}$.\n\n.. note:: This gate is not unitary."); + DEF_GATE( + P1Gate, + "Specific class of projection gate to $\\ket{1}$.\n\n.. note:: This gate is not unitary."); + +#define DEF_ROTATION_GATE(GATE_TYPE, DESCRIPTION) \ + DEF_GATE(GATE_TYPE, DESCRIPTION) \ + .def( \ + "angle", [](const GATE_TYPE& gate) { return gate->angle(); }, "Get `angle` property.") + + DEF_ROTATION_GATE( + RXGate, + "Specific class of X rotation gate, represented as $e^{-i\\frac{\\mathrm{angle}}{2}X}$."); + DEF_ROTATION_GATE( + RYGate, + "Specific class of Y rotation gate, represented as $e^{-i\\frac{\\mathrm{angle}}{2}Y}$."); + DEF_ROTATION_GATE( + RZGate, + "Specific class of Z rotation gate, represented as $e^{-i\\frac{\\mathrm{angle}}{2}Z}$."); + + DEF_GATE(U1Gate, + "Specific class of IBMQ's U1 Gate, which is a rotation abount Z-axis, " + "represented as " + "$\\begin{bmatrix}\n1 & 0\\\\\n0 & e^{i\\lambda}\n\\end{bmatrix}$.") + .def( + "lambda_", [](const U1Gate& gate) { return gate->lambda(); }, "Get `lambda` property."); + DEF_GATE(U2Gate, + "Specific class of IBMQ's U2 Gate, which is a rotation about X+Z-axis, " + "represented as " + "$\\frac{1}{\\sqrt{2}} \\begin{bmatrix}1 & -e^{-i\\lambda}\\\\\n" + "e^{i\\phi} & e^{i(\\phi+\\lambda)}\n\\end{bmatrix}$.") + .def( + "phi", [](const U2Gate& gate) { return gate->phi(); }, "Get `phi` property.") + .def( + "lambda_", [](const U2Gate& gate) { return gate->lambda(); }, "Get `lambda` property."); + DEF_GATE(U3Gate, + "Specific class of IBMQ's U3 Gate, which is a rotation abount 3 axis, " + "represented as " + "$\\begin{bmatrix}\n\\cos \\frac{\\theta}{2} & " + "-e^{i\\lambda}\\sin\\frac{\\theta}{2}\\\\\n" + "e^{i\\phi}\\sin\\frac{\\theta}{2} & " + "e^{i(\\phi+\\lambda)}\\cos\\frac{\\theta}{2}\n\\end{bmatrix}$.") + .def( + "theta", [](const U3Gate& gate) { return gate->theta(); }, "Get `theta` property.") + .def( + "phi", [](const U3Gate& gate) { return gate->phi(); }, "Get `phi` property.") + .def( + "lambda_", [](const U3Gate& gate) { return gate->lambda(); }, "Get `lambda` property."); + DEF_GATE(SwapGate, "Specific class of two-qubit swap gate."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/merge_gate.cpp b/scaluq/gate/merge_gate.cpp index 844984fc..158886b9 100644 --- a/scaluq/gate/merge_gate.cpp +++ b/scaluq/gate/merge_gate.cpp @@ -32,21 +32,21 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { assert(!pauli_id1 || pauli_id1 != 0); assert(!pauli_id2 || pauli_id2 != 0); if (pauli_id1 && pauli_id2) { - std::uint64_t target1 = gate1->get_target_qubit_list()[0]; - std::uint64_t target2 = gate2->get_target_qubit_list()[0]; + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; if (target1 == target2) { if (pauli_id1 == pauli_id2) return {gate::I(), 0.}; if (pauli_id1 == 1) { - if (pauli_id2 == 2) return {gate::Z(target1), -PI() / 2}; - if (pauli_id2 == 3) return {gate::Y(target1), PI() / 2}; + if (pauli_id2 == 2) return {gate::Z(target1), -Kokkos::numbers::pi / 2}; + if (pauli_id2 == 3) return {gate::Y(target1), Kokkos::numbers::pi / 2}; } if (pauli_id1 == 2) { - if (pauli_id2 == 3) return {gate::X(target1), -PI() / 2}; - if (pauli_id2 == 1) return {gate::Z(target1), PI() / 2}; + if (pauli_id2 == 3) return {gate::X(target1), -Kokkos::numbers::pi / 2}; + if (pauli_id2 == 1) return {gate::Z(target1), Kokkos::numbers::pi / 2}; } if (pauli_id1 == 3) { - if (pauli_id2 == 1) return {gate::Y(target1), -PI() / 2}; - if (pauli_id2 == 2) return {gate::X(target1), PI() / 2}; + if (pauli_id2 == 1) return {gate::Y(target1), -Kokkos::numbers::pi / 2}; + if (pauli_id2 == 2) return {gate::X(target1), Kokkos::numbers::pi / 2}; } } } @@ -54,11 +54,11 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { (pauli_id2 || gate2.gate_type() == GateType::Pauli)) { auto pauli1 = gate_type1 == GateType::Pauli ? PauliGate(gate1)->pauli() - : PauliOperator(std::vector{gate1->get_target_qubit_list()[0]}, + : PauliOperator(std::vector{gate1->target_qubit_list()[0]}, std::vector{pauli_id1.value()}); auto pauli2 = gate_type2 == GateType::Pauli ? PauliGate(gate2)->pauli() - : PauliOperator(std::vector{gate2->get_target_qubit_list()[0]}, + : PauliOperator(std::vector{gate2->target_qubit_list()[0]}, std::vector{pauli_id2.value()}); return {gate::Pauli(pauli2 * pauli1), 0.}; } @@ -87,8 +87,8 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { auto oct_phase1 = get_oct_phase(gate_type1); auto oct_phase2 = get_oct_phase(gate_type2); if (oct_phase1 && oct_phase2) { - std::uint64_t target1 = gate1->get_target_qubit_list()[0]; - std::uint64_t target2 = gate2->get_target_qubit_list()[0]; + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; if (target1 == target2) { auto g = oct_phase_gate(oct_phase1.value() + oct_phase2.value(), target1); if (g) return {g.value(), 0.}; @@ -96,14 +96,14 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { } if ((oct_phase1 || gate_type1 == GateType::RZ || gate_type1 == GateType::U1) && (oct_phase2 || gate_type2 == GateType::RZ || gate_type2 == GateType::U1)) { - std::uint64_t target1 = gate1->get_target_qubit_list()[0]; - std::uint64_t target2 = gate2->get_target_qubit_list()[0]; + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; if (target1 == target2) { - double phase1 = oct_phase1 ? oct_phase1.value() * PI() / 4 + double phase1 = oct_phase1 ? oct_phase1.value() * Kokkos::numbers::pi / 4 : gate_type1 == GateType::RZ ? RZGate(gate1)->angle() : U1Gate(gate1)->lambda(); double global_phase1 = gate_type1 == GateType::RZ ? -RZGate(gate1)->angle() / 2 : 0.; - double phase2 = oct_phase2 ? oct_phase2.value() * PI() / 4 + double phase2 = oct_phase2 ? oct_phase2.value() * Kokkos::numbers::pi / 4 : gate_type2 == GateType::RZ ? RZGate(gate2)->angle() : U1Gate(gate2)->lambda(); double global_phase2 = gate_type2 == GateType::RZ ? -RZGate(gate2)->angle() / 2 : 0.; @@ -114,17 +114,17 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { // Special case: RX auto get_rx_angle = [&](Gate gate, GateType gate_type) -> std::optional { if (gate_type == GateType::I) return 0.; - if (gate_type == GateType::X) return PI(); - if (gate_type == GateType::SqrtX) return PI() / 2; - if (gate_type == GateType::SqrtXdag) return -PI() / 2; + if (gate_type == GateType::X) return Kokkos::numbers::pi; + if (gate_type == GateType::SqrtX) return Kokkos::numbers::pi / 2; + if (gate_type == GateType::SqrtXdag) return -Kokkos::numbers::pi / 2; if (gate_type == GateType::RX) return RXGate(gate)->angle(); return std::nullopt; }; auto rx_param1 = get_rx_angle(gate1, gate_type1); auto rx_param2 = get_rx_angle(gate2, gate_type2); if (rx_param1 && rx_param2) { - std::uint64_t target1 = gate1->get_target_qubit_list()[0]; - std::uint64_t target2 = gate2->get_target_qubit_list()[0]; + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; double global_phase1 = gate_type1 == GateType::RX ? 0. : rx_param1.value() / 2; double global_phase2 = gate_type2 == GateType::RX ? 0. : rx_param2.value() / 2; if (target1 == target2) { @@ -136,17 +136,17 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { // Special case: RY auto get_ry_angle = [&](Gate gate, GateType gate_type) -> std::optional { if (gate_type == GateType::I) return 0.; - if (gate_type == GateType::Y) return PI(); - if (gate_type == GateType::SqrtY) return PI() / 2; - if (gate_type == GateType::SqrtYdag) return -PI() / 2; + if (gate_type == GateType::Y) return Kokkos::numbers::pi; + if (gate_type == GateType::SqrtY) return Kokkos::numbers::pi / 2; + if (gate_type == GateType::SqrtYdag) return -Kokkos::numbers::pi / 2; if (gate_type == GateType::RY) return RYGate(gate)->angle(); return std::nullopt; }; auto ry_param1 = get_ry_angle(gate1, gate_type1); auto ry_param2 = get_ry_angle(gate2, gate_type2); if (ry_param1 && ry_param2) { - std::uint64_t target1 = gate1->get_target_qubit_list()[0]; - std::uint64_t target2 = gate2->get_target_qubit_list()[0]; + std::uint64_t target1 = gate1->target_qubit_list()[0]; + std::uint64_t target2 = gate2->target_qubit_list()[0]; double global_phase1 = gate_type1 == GateType::RY ? 0. : ry_param1.value() / 2; double global_phase2 = gate_type2 == GateType::RY ? 0. : ry_param2.value() / 2; if (target1 == target2) { @@ -188,10 +188,10 @@ std::pair merge_gate(const Gate& gate1, const Gate& gate2) { } // General case - auto gate1_targets = gate1->get_target_qubit_list(); - std::ranges::copy(gate1->get_control_qubit_list(), std::back_inserter(gate1_targets)); - auto gate2_targets = gate2->get_target_qubit_list(); - std::ranges::copy(gate2->get_control_qubit_list(), std::back_inserter(gate2_targets)); + auto gate1_targets = gate1->target_qubit_list(); + std::ranges::copy(gate1->control_qubit_list(), std::back_inserter(gate1_targets)); + auto gate2_targets = gate2->target_qubit_list(); + std::ranges::copy(gate2->control_qubit_list(), std::back_inserter(gate2_targets)); std::vector merged_targets(gate1_targets.size() + gate2_targets.size()); std::ranges::copy(gate1_targets, merged_targets.begin()); std::ranges::copy(gate2_targets, merged_targets.begin() + gate1_targets.size()); diff --git a/scaluq/gate/merge_gate.hpp b/scaluq/gate/merge_gate.hpp index c92281dc..55130d1f 100644 --- a/scaluq/gate/merge_gate.hpp +++ b/scaluq/gate/merge_gate.hpp @@ -4,4 +4,13 @@ namespace scaluq { std::pair merge_gate(const Gate& gate1, const Gate& gate2); + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_merge_gate_hpp(nb::module_& m) { + m.def( + "merge_gate", &merge_gate, "Merge two gates. return value is (merged gate, global phase)."); } +} // namespace internal +#endif +} // namespace scaluq diff --git a/scaluq/gate/param_gate.hpp b/scaluq/gate/param_gate.hpp index e0ff1cb8..d1e86050 100644 --- a/scaluq/gate/param_gate.hpp +++ b/scaluq/gate/param_gate.hpp @@ -12,28 +12,44 @@ class ParamGateBase; template concept ParamGateImpl = std::derived_from; -class PRXGateImpl; -class PRYGateImpl; -class PRZGateImpl; -class PPauliRotationGateImpl; +class ParamRXGateImpl; +class ParamRYGateImpl; +class ParamRZGateImpl; +class ParamPauliRotationGateImpl; +class ParamProbablisticGateImpl; template class ParamGatePtr; } // namespace internal using ParamGate = internal::ParamGatePtr; -enum class ParamGateType { Unknown, PRX, PRY, PRZ, PPauliRotation }; +enum class ParamGateType { + Unknown, + ParamRX, + ParamRY, + ParamRZ, + ParamPauliRotation, + ParamProbablistic, + Error +}; template constexpr ParamGateType get_param_gate_type() { - if constexpr (std::is_same_v) return ParamGateType::Unknown; - if constexpr (std::is_same_v) return ParamGateType::PRX; - if constexpr (std::is_same_v) return ParamGateType::PRY; - if constexpr (std::is_same_v) return ParamGateType::PRZ; - if constexpr (std::is_same_v) - return ParamGateType::PPauliRotation; - static_assert("unknown ParamGateImpl"); - return ParamGateType::Unknown; + using TWithoutConst = std::remove_cv_t; + if constexpr (std::is_same_v) + return ParamGateType::Unknown; + else if constexpr (std::is_same_v) + return ParamGateType::ParamRX; + else if constexpr (std::is_same_v) + return ParamGateType::ParamRY; + else if constexpr (std::is_same_v) + return ParamGateType::ParamRZ; + else if constexpr (std::is_same_v) + return ParamGateType::ParamPauliRotation; + else if constexpr (std::is_same_v) + return ParamGateType::ParamProbablistic; + else + static_assert(internal::lazy_false_v, "unknown GateImpl"); } namespace internal { @@ -50,39 +66,56 @@ class ParamGateBase { } } + std::string get_qubit_info_as_string(const std::string& indent) const { + std::ostringstream ss; + auto targets = target_qubit_list(); + auto controls = control_qubit_list(); + ss << indent << " Parameter Coefficient: " << _pcoef << "\n"; + ss << indent << " Target Qubits: {"; + for (std::uint32_t i = 0; i < targets.size(); ++i) + ss << targets[i] << (i == targets.size() - 1 ? "" : ", "); + ss << "}\n"; + ss << indent << " Control Qubits: {"; + for (std::uint32_t i = 0; i < controls.size(); ++i) + ss << controls[i] << (i == controls.size() - 1 ? "" : ", "); + ss << "}"; + return ss.str(); + } + public: - ParamGateBase(std::uint64_t target_mask, std::uint64_t control_mask, double pcoef = 1.) - : _target_mask(target_mask), _control_mask(control_mask), _pcoef(pcoef) { + ParamGateBase(std::uint64_t target_mask, std::uint64_t control_mask, double param_coef = 1.) + : _target_mask(target_mask), _control_mask(control_mask), _pcoef(param_coef) { if (_target_mask & _control_mask) [[unlikely]] { throw std::runtime_error( "Error: ParamGate::ParamGate(std::uint64_t target_mask, std::uint64_t " - "control_mask) : Target and " - "control qubits must not overlap."); + "control_mask) : Target and control qubits must not overlap."); } } virtual ~ParamGateBase() = default; - [[nodiscard]] double pcoef() { return _pcoef; } + [[nodiscard]] double param_coef() const { return _pcoef; } - [[nodiscard]] virtual std::vector get_target_qubit_list() const { + [[nodiscard]] virtual std::vector target_qubit_list() const { return mask_to_vector(_target_mask); } - [[nodiscard]] virtual std::vector get_control_qubit_list() const { + [[nodiscard]] virtual std::vector control_qubit_list() const { return mask_to_vector(_control_mask); } - [[nodiscard]] virtual std::vector get_operand_qubit_list() const { + [[nodiscard]] virtual std::vector operand_qubit_list() const { return mask_to_vector(_target_mask | _control_mask); } - [[nodiscard]] virtual std::uint64_t get_target_qubit_mask() const { return _target_mask; } - [[nodiscard]] virtual std::uint64_t get_control_qubit_mask() const { return _control_mask; } - [[nodiscard]] virtual std::uint64_t get_operand_qubit_mask() const { + [[nodiscard]] virtual std::uint64_t target_qubit_mask() const { return _target_mask; } + [[nodiscard]] virtual std::uint64_t control_qubit_mask() const { return _control_mask; } + [[nodiscard]] virtual std::uint64_t operand_qubit_mask() const { return _target_mask | _control_mask; } [[nodiscard]] virtual ParamGate get_inverse() const = 0; - [[nodiscard]] virtual ComplexMatrix get_matrix(double param) const = 0; + [[nodiscard]] virtual internal::ComplexMatrix get_matrix(double param) const = 0; virtual void update_quantum_state(StateVector& state_vector, double param) const = 0; + + [[nodiscard]] virtual std::string to_string(const std::string& indent = "") const = 0; }; template @@ -99,7 +132,7 @@ class ParamGatePtr { ParamGatePtr() : _param_gate_ptr(nullptr), _param_gate_type(get_param_gate_type()) {} ParamGatePtr(const ParamGatePtr& param_gate) = default; template - ParamGatePtr(const std::shared_ptr& param_gate_ptr) { + ParamGatePtr(const std::shared_ptr& param_gate_ptr) { if constexpr (std::is_same_v) { _param_gate_type = get_param_gate_type(); _param_gate_ptr = param_gate_ptr; @@ -142,7 +175,89 @@ class ParamGatePtr { } return _param_gate_ptr.get(); } + + friend std::ostream& operator<<(std::ostream& os, ParamGatePtr gate) { + os << gate->to_string(); + return os; + } }; } // namespace internal +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +#define DEF_PARAM_GATE_BASE(PARAM_GATE_TYPE, DESCRIPTION) \ + nb::class_(m, #PARAM_GATE_TYPE, DESCRIPTION) \ + .def("param_gate_type", \ + &PARAM_GATE_TYPE::param_gate_type, \ + "Get parametric gate type as `ParamGateType` enum.") \ + .def( \ + "param_coef", \ + [](const PARAM_GATE_TYPE& gate) { return gate->param_coef(); }, \ + "Get coefficient of parameter.") \ + .def( \ + "target_qubit_list", \ + [](const PARAM_GATE_TYPE& gate) { return gate->target_qubit_list(); }, \ + "Get target qubits as `list[int]`. **Control qubits is not included.**") \ + .def( \ + "control_qubit_list", \ + [](const PARAM_GATE_TYPE& gate) { return gate->control_qubit_list(); }, \ + "Get control qubits as `list[int]`.") \ + .def( \ + "operand_qubit_list", \ + [](const PARAM_GATE_TYPE& gate) { return gate->operand_qubit_list(); }, \ + "Get target and control qubits as `list[int]`.") \ + .def( \ + "target_qubit_mask", \ + [](const PARAM_GATE_TYPE& gate) { return gate->target_qubit_mask(); }, \ + "Get target qubits as mask. **Control qubits is not included.**") \ + .def( \ + "control_qubit_mask", \ + [](const PARAM_GATE_TYPE& gate) { return gate->control_qubit_mask(); }, \ + "Get control qubits as mask.") \ + .def( \ + "operand_qubit_mask", \ + [](const PARAM_GATE_TYPE& gate) { return gate->operand_qubit_mask(); }, \ + "Get target and control qubits as mask.") \ + .def( \ + "get_inverse", \ + [](const PARAM_GATE_TYPE& param_gate) { return param_gate->get_inverse(); }, \ + "Generate inverse parametric-gate as `ParamGate` type. If not exists, return None.") \ + .def( \ + "update_quantum_state", \ + [](const PARAM_GATE_TYPE& param_gate, StateVector& state_vector, double param) { \ + param_gate->update_quantum_state(state_vector, param); \ + }, \ + "Apply gate to `state_vector` with holding the parameter. `state_vector` in args is " \ + "directly updated.") \ + .def( \ + "get_matrix", \ + [](const PARAM_GATE_TYPE& gate, double param) { return gate->get_matrix(param); }, \ + "Get matrix representation of the gate with holding the parameter.") + +nb::class_ param_gate_base_def; + +#define DEF_PARAM_GATE(PARAM_GATE_TYPE, DESCRIPTION) \ + ::scaluq::internal::param_gate_base_def.def(nb::init(), \ + "Upcast from `" #PARAM_GATE_TYPE "`."); \ + DEF_PARAM_GATE_BASE( \ + PARAM_GATE_TYPE, \ + DESCRIPTION \ + "\n\n.. note:: Upcast is required to use gate-general functions (ex: add to Circuit).") \ + .def(nb::init()) + +void bind_gate_param_gate_hpp(nb::module_& m) { + nb::enum_(m, "ParamGateType", "Enum of ParamGate Type.") + .value("ParamRX", ParamGateType::ParamRX) + .value("ParamRY", ParamGateType::ParamRY) + .value("ParamRZ", ParamGateType::ParamRZ) + .value("ParamPauliRotation", ParamGateType::ParamPauliRotation); + + param_gate_base_def = DEF_PARAM_GATE_BASE( + ParamGate, + "General class of parametric quantum gate.\n\n.. note:: Downcast to requred to use " + "gate-specific functions."); +} + +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/param_gate_factory.hpp b/scaluq/gate/param_gate_factory.hpp index cfbd3e2f..bc7a7ec6 100644 --- a/scaluq/gate/param_gate_factory.hpp +++ b/scaluq/gate/param_gate_factory.hpp @@ -15,35 +15,82 @@ class ParamGateFactory { }; } // namespace internal namespace gate { -inline ParamGate PRX(std::uint64_t target, - double pcoef = 1., - const std::vector& controls = {}) { - return internal::ParamGateFactory::create_gate( - internal::vector_to_mask({target}), internal::vector_to_mask(controls), pcoef); +inline ParamGate ParamRX(std::uint64_t target, + double param_coef = 1., + const std::vector& controls = {}) { + return internal::ParamGateFactory::create_gate( + internal::vector_to_mask({target}), internal::vector_to_mask(controls), param_coef); } -inline ParamGate PRY(std::uint64_t target, - double pcoef = 1., - const std::vector& controls = {}) { - return internal::ParamGateFactory::create_gate( - internal::vector_to_mask({target}), internal::vector_to_mask(controls), pcoef); +inline ParamGate ParamRY(std::uint64_t target, + double param_coef = 1., + const std::vector& controls = {}) { + return internal::ParamGateFactory::create_gate( + internal::vector_to_mask({target}), internal::vector_to_mask(controls), param_coef); } -inline ParamGate PRZ(std::uint64_t target, - double pcoef = 1., - const std::vector& controls = {}) { - return internal::ParamGateFactory::create_gate( - internal::vector_to_mask({target}), internal::vector_to_mask(controls), pcoef); +inline ParamGate ParamRZ(std::uint64_t target, + double param_coef = 1., + const std::vector& controls = {}) { + return internal::ParamGateFactory::create_gate( + internal::vector_to_mask({target}), internal::vector_to_mask(controls), param_coef); } -// まだ -inline ParamGate PPauliRotation(const PauliOperator& pauli, - double pcoef = 1., - const std::vector& controls = {}) { - return internal::ParamGateFactory::create_gate( - internal::vector_to_mask(controls), pauli, pcoef); +inline ParamGate ParamPauliRotation(const PauliOperator& pauli, + double param_coef = 1., + const std::vector& controls = {}) { + return internal::ParamGateFactory::create_gate( + internal::vector_to_mask(controls), pauli, param_coef); } -inline ParamGate PProbablistic(const std::vector& distribution, - const std::vector>& gate_list) { - return internal::ParamGateFactory::create_gate(distribution, - gate_list); +inline ParamGate ParamProbablistic(const std::vector& distribution, + const std::vector>& gate_list) { + return internal::ParamGateFactory::create_gate( + distribution, gate_list); } } // namespace gate + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_param_gate_factory(nb::module_& mgate) { + mgate.def("ParamRX", + &gate::ParamRX, + "Generate general ParamGate class instance of ParamRX.", + "target"_a, + "coef"_a = 1., + "controls"_a = std::vector{}); + mgate.def("ParamRY", + &gate::ParamRY, + "Generate general ParamGate class instance of ParamRY.", + "target"_a, + "coef"_a = 1., + "controls"_a = std::vector{}); + mgate.def("ParamRZ", + &gate::ParamRZ, + "Generate general ParamGate class instance of ParamRZ.", + "target"_a, + "coef"_a = 1., + "controls"_a = std::vector{}); + mgate.def("ParamPauliRotation", + &gate::ParamPauliRotation, + "Generate general ParamGate class instance of ParamPauliRotation.", + "pauli"_a, + "coef"_a = 1., + "controls"_a = std::vector{}); + mgate.def("ParamProbablistic", + &gate::ParamProbablistic, + "Generate general ParamGate class instance of ParamProbablistic."); + mgate.def( + "ParamProbablistic", + [](const std::vector>>& prob_gate_list) { + std::vector distribution; + std::vector> gate_list; + distribution.reserve(prob_gate_list.size()); + gate_list.reserve(prob_gate_list.size()); + for (const auto& [prob, gate] : prob_gate_list) { + distribution.push_back(prob); + gate_list.push_back(gate); + } + return gate::ParamProbablistic(distribution, gate_list); + }, + "Generate general ParamGate class instance of ParamProbablistic."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/param_gate_pauli.hpp b/scaluq/gate/param_gate_pauli.hpp index 0ba0a736..f0dffdf7 100644 --- a/scaluq/gate/param_gate_pauli.hpp +++ b/scaluq/gate/param_gate_pauli.hpp @@ -2,43 +2,73 @@ #include +#include "../operator/apply_pauli.hpp" #include "../operator/pauli_operator.hpp" #include "../util/utility.hpp" #include "param_gate.hpp" namespace scaluq { namespace internal { -class PPauliRotationGateImpl : public ParamGateBase { +class ParamPauliRotationGateImpl : public ParamGateBase { const PauliOperator _pauli; public: - PPauliRotationGateImpl(std::uint64_t control_mask, - const PauliOperator& pauli, - double pcoef = 1.) - : ParamGateBase(vector_to_mask(pauli.get_target_qubit_list()), control_mask, pcoef), + ParamPauliRotationGateImpl(std::uint64_t control_mask, + const PauliOperator& pauli, + double param_coef = 1.) + : ParamGateBase(vector_to_mask(pauli.target_qubit_list()), control_mask, param_coef), _pauli(pauli) {} PauliOperator pauli() const { return _pauli; } - std::vector get_pauli_id_list() const { return _pauli.get_pauli_id_list(); } + std::vector pauli_id_list() const { return _pauli.pauli_id_list(); } ParamGate get_inverse() const override { - return std::make_shared(_control_mask, _pauli, -_pcoef); + return std::make_shared(_control_mask, _pauli, -_pcoef); } - ComplexMatrix get_matrix(double param) const override { + internal::ComplexMatrix get_matrix(double param) const override { double angle = _pcoef * param; - Complex true_angle = angle * this->_pauli.get_coef(); - ComplexMatrix mat = this->_pauli.get_matrix_ignoring_coef(); + Complex true_angle = angle * this->_pauli.coef(); + internal::ComplexMatrix mat = this->_pauli.get_matrix_ignoring_coef(); StdComplex imag_unit(0, 1); mat = (StdComplex)Kokkos::cos(-true_angle / 2) * - ComplexMatrix::Identity(mat.rows(), mat.cols()) + + internal::ComplexMatrix::Identity(mat.rows(), mat.cols()) + imag_unit * (StdComplex)Kokkos::sin(-true_angle / 2) * mat; return mat; } void update_quantum_state(StateVector& state_vector, double param) const override { - pauli_rotation_gate(_control_mask, _pauli, _pcoef * param, state_vector); + auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation(); + apply_pauli_rotation(_control_mask, + bit_flip_mask, + phase_flip_mask, + _pauli.coef(), + _pcoef * param, + state_vector); + } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + auto controls = control_qubit_list(); + ss << indent << "Gate Type: ParamPauliRotation\n"; + ss << indent << " Control Qubits: {"; + for (std::uint32_t i = 0; i < controls.size(); ++i) + ss << controls[i] << (i == controls.size() - 1 ? "" : ", "); + ss << "}\n"; + ss << indent << " Pauli Operator: \"" << _pauli.get_pauli_string() << "\""; + return ss.str(); } }; } // namespace internal -using PPauliRotationGate = internal::ParamGatePtr; +using ParamPauliRotationGate = internal::ParamGatePtr; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_param_gate_pauli_hpp(nb::module_& m) { + DEF_PARAM_GATE( + ParamPauliRotationGate, + "Specific class of parametric multi-qubit pauli-rotation gate, represented as " + "$e^{-i\\frac{\\mathrm{angle}}{2}P}$. `angle` is given as `param * param_coef`."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/param_gate_probablistic.hpp b/scaluq/gate/param_gate_probablistic.hpp index 063a1ea9..d47b0e2e 100644 --- a/scaluq/gate/param_gate_probablistic.hpp +++ b/scaluq/gate/param_gate_probablistic.hpp @@ -8,15 +8,15 @@ namespace scaluq { namespace internal { -class PProbablisticGateImpl : public ParamGateBase { +class ParamProbablisticGateImpl : public ParamGateBase { using EitherGate = std::variant; std::vector _distribution; std::vector _cumlative_distribution; std::vector _gate_list; public: - PProbablisticGateImpl(const std::vector& distribution, - const std::vector>& gate_list) + ParamProbablisticGateImpl(const std::vector& distribution, + const std::vector>& gate_list) : ParamGateBase(0, 0), _distribution(distribution), _gate_list(gate_list) { std::uint64_t n = distribution.size(); if (n == 0) { @@ -35,35 +35,35 @@ class PProbablisticGateImpl : public ParamGateBase { const std::vector>& gate_list() const { return _gate_list; } const std::vector& distribution() const { return _distribution; } - std::vector get_target_qubit_list() const override { + std::vector target_qubit_list() const override { throw std::runtime_error( - "PProbablisticGateImpl::get_target_qubit_list(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::target_qubit_list(): This function must not be used in " + "ParamProbablisticGateImpl."); } - std::vector get_control_qubit_list() const override { + std::vector control_qubit_list() const override { throw std::runtime_error( - "PProbablisticGateImpl::get_control_qubit_list(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::control_qubit_list(): This function must not be used in " + "ParamProbablisticGateImpl."); } - std::vector get_operand_qubit_list() const override { + std::vector operand_qubit_list() const override { throw std::runtime_error( - "PProbablisticGateImpl::get_operand_qubit_list(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::operand_qubit_list(): This function must not be used in " + "ParamProbablisticGateImpl."); } - std::uint64_t get_target_qubit_mask() const override { + std::uint64_t target_qubit_mask() const override { throw std::runtime_error( - "PProbablisticGateImpl::get_target_qubit_mask(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::target_qubit_mask(): This function must not be used in " + "ParamProbablisticGateImpl."); } - std::uint64_t get_control_qubit_mask() const override { + std::uint64_t control_qubit_mask() const override { throw std::runtime_error( - "PProbablisticGateImpl::get_control_qubit_mask(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::control_qubit_mask(): This function must not be used in " + "ParamProbablisticGateImpl."); } - std::uint64_t get_operand_qubit_mask() const override { + std::uint64_t operand_qubit_mask() const override { throw std::runtime_error( - "PProbablisticGateImpl::get_operand_qubit_mask(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::operand_qubit_mask(): This function must not be used in " + "ParamProbablisticGateImpl."); } ParamGate get_inverse() const override { @@ -73,12 +73,12 @@ class PProbablisticGateImpl : public ParamGateBase { _gate_list, std::back_inserter(inv_gate_list), [](const EitherGate& gate) { return std::visit([](const auto& g) { return EitherGate{g->get_inverse()}; }, gate); }); - return std::make_shared(_distribution, inv_gate_list); + return std::make_shared(_distribution, inv_gate_list); } - ComplexMatrix get_matrix(double) const override { + internal::ComplexMatrix get_matrix(double) const override { throw std::runtime_error( - "PProbablisticGateImpl::get_matrix(): This function must not be used in " - "PProbablisticGateImpl."); + "ParamProbablisticGateImpl::get_matrix(): This function must not be used in " + "ParamProbablisticGateImpl."); } void update_quantum_state(StateVector& state_vector, double param) const override { @@ -95,8 +95,44 @@ class PProbablisticGateImpl : public ParamGateBase { std::get<1>(gate)->update_quantum_state(state_vector, param); } } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + const auto dist = distribution(); + ss << indent << "Gate Type: Probablistic\n"; + for (std::size_t i = 0; i < dist.size(); ++i) { + ss << indent << " --------------------\n"; + ss << indent << " Probability: " << dist[i] << "\n"; + std::visit( + [&](auto&& arg) { + ss << arg->to_string(indent + " ") << (i == dist.size() - 1 ? "" : "\n"); + }, + gate_list()[i]); + } + return ss.str(); + } }; } // namespace internal -using PProbablisticGate = internal::ParamGatePtr; +using ParamProbablisticGate = internal::ParamGatePtr; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_param_gate_probablistic_hpp(nb::module_& m) { + DEF_PARAM_GATE( + ParamProbablisticGate, + "Specific class of parametric probablistic gate. The gate to apply is picked from a " + "cirtain " + "distribution.") + .def( + "gate_list", + [](const ParamProbablisticGate& gate) { return gate->gate_list(); }, + nb::rv_policy::reference) + .def( + "distribution", + [](const ParamProbablisticGate& gate) { return gate->distribution(); }, + nb::rv_policy::reference); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/param_gate_standard.hpp b/scaluq/gate/param_gate_standard.hpp index cf12afae..755a162c 100644 --- a/scaluq/gate/param_gate_standard.hpp +++ b/scaluq/gate/param_gate_standard.hpp @@ -8,16 +8,16 @@ namespace scaluq { namespace internal { -class PRXGateImpl : public ParamGateBase { +class ParamRXGateImpl : public ParamGateBase { public: using ParamGateBase::ParamGateBase; ParamGate get_inverse() const override { - return std::make_shared(_target_mask, _control_mask, -_pcoef); + return std::make_shared(_target_mask, _control_mask, -_pcoef); } - ComplexMatrix get_matrix(double param) const override { + internal::ComplexMatrix get_matrix(double param) const override { double angle = _pcoef * param; - ComplexMatrix mat(2, 2); + internal::ComplexMatrix mat(2, 2); mat << std::cos(angle / 2), -1i * std::sin(angle / 2), -1i * std::sin(angle / 2), std::cos(angle / 2); return mat; @@ -27,18 +27,25 @@ class PRXGateImpl : public ParamGateBase { check_qubit_mask_within_bounds(state_vector); rx_gate(_target_mask, _control_mask, _pcoef * param, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: ParamRX\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; -class PRYGateImpl : public ParamGateBase { +class ParamRYGateImpl : public ParamGateBase { public: using ParamGateBase::ParamGateBase; ParamGate get_inverse() const override { - return std::make_shared(_target_mask, _control_mask, -_pcoef); + return std::make_shared(_target_mask, _control_mask, -_pcoef); } - ComplexMatrix get_matrix(double param) const override { + internal::ComplexMatrix get_matrix(double param) const override { double angle = _pcoef * param; - ComplexMatrix mat(2, 2); + internal::ComplexMatrix mat(2, 2); mat << std::cos(angle / 2), -std::sin(angle / 2), std::sin(angle / 2), std::cos(angle / 2); return mat; } @@ -47,18 +54,25 @@ class PRYGateImpl : public ParamGateBase { check_qubit_mask_within_bounds(state_vector); ry_gate(_target_mask, _control_mask, _pcoef * param, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: ParamRY\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; -class PRZGateImpl : public ParamGateBase { +class ParamRZGateImpl : public ParamGateBase { public: using ParamGateBase::ParamGateBase; ParamGate get_inverse() const override { - return std::make_shared(_target_mask, _control_mask, -_pcoef); + return std::make_shared(_target_mask, _control_mask, -_pcoef); } - ComplexMatrix get_matrix(double param) const override { + internal::ComplexMatrix get_matrix(double param) const override { double angle = param * _pcoef; - ComplexMatrix mat(2, 2); + internal::ComplexMatrix mat(2, 2); mat << std::exp(-0.5i * angle), 0, 0, std::exp(0.5i * angle); return mat; } @@ -67,12 +81,37 @@ class PRZGateImpl : public ParamGateBase { check_qubit_mask_within_bounds(state_vector); rz_gate(_target_mask, _control_mask, _pcoef * param, state_vector); } + + std::string to_string(const std::string& indent) const override { + std::ostringstream ss; + ss << indent << "Gate Type: ParamRZ\n"; + ss << get_qubit_info_as_string(indent); + return ss.str(); + } }; } // namespace internal -using PRXGate = internal::ParamGatePtr; -using PRYGate = internal::ParamGatePtr; -using PRZGate = internal::ParamGatePtr; +using ParamRXGate = internal::ParamGatePtr; +using ParamRYGate = internal::ParamGatePtr; +using ParamRZGate = internal::ParamGatePtr; +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_gate_param_gate_standard_hpp(nb::module_& m) { + DEF_PARAM_GATE( + ParamRXGate, + "Specific class of parametric X rotation gate, represented as " + "$e^{-i\\frac{\\mathrm{angle}}{2}X}$. `angle` is given as `param * param_coef`."); + DEF_PARAM_GATE( + ParamRYGate, + "Specific class of parametric Y rotation gate, represented as " + "$e^{-i\\frac{\\mathrm{angle}}{2}Y}$. `angle` is given as `param * param_coef`."); + DEF_PARAM_GATE( + ParamRZGate, + "Specific class of parametric Z rotation gate, represented as " + "$e^{-i\\frac{\\mathrm{angle}}{2}Z}$. `angle` is given as `param * param_coef`."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/gate/update_ops.hpp b/scaluq/gate/update_ops.hpp index a3ceeba9..3af48e23 100644 --- a/scaluq/gate/update_ops.hpp +++ b/scaluq/gate/update_ops.hpp @@ -57,21 +57,21 @@ void rz_gate(std::uint64_t target_mask, double angle, StateVector& state); -matrix_2_2 get_IBMQ_matrix(double _theta, double _phi, double _lambda); +Matrix2x2 get_IBMQ_matrix(double _theta, double _phi, double _lambda); void one_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, - const matrix_2_2& matrix, + const Matrix2x2& matrix, StateVector& state); void two_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, - const matrix_4_4& matrix, + const Matrix4x4& matrix, StateVector& state); void one_target_diagonal_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, - const diagonal_matrix_2_2& diag, + const DiagonalMatrix2x2& diag, StateVector& state); void u1_gate(std::uint64_t target_mask, diff --git a/scaluq/gate/update_ops_dense_matrix.cpp b/scaluq/gate/update_ops_dense_matrix.cpp index b8d1b2de..cfabc059 100644 --- a/scaluq/gate/update_ops_dense_matrix.cpp +++ b/scaluq/gate/update_ops_dense_matrix.cpp @@ -9,7 +9,7 @@ namespace scaluq { namespace internal { void one_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, - const matrix_2_2& matrix, + const Matrix2x2& matrix, StateVector& state) { Kokkos::parallel_for( state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) { @@ -18,8 +18,8 @@ void one_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t basis_1 = basis_0 | target_mask; Complex val0 = state._raw[basis_0]; Complex val1 = state._raw[basis_1]; - Complex res0 = matrix.val[0][0] * val0 + matrix.val[0][1] * val1; - Complex res1 = matrix.val[1][0] * val0 + matrix.val[1][1] * val1; + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1; state._raw[basis_0] = res0; state._raw[basis_1] = res1; }); @@ -28,7 +28,7 @@ void one_target_dense_matrix_gate(std::uint64_t target_mask, void two_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, - const matrix_4_4& matrix, + const Matrix4x4& matrix, StateVector& state) { std::uint64_t lower_target_mask = -target_mask & target_mask; std::uint64_t upper_target_mask = target_mask ^ lower_target_mask; @@ -44,14 +44,14 @@ void two_target_dense_matrix_gate(std::uint64_t target_mask, Complex val1 = state._raw[basis_1]; Complex val2 = state._raw[basis_2]; Complex val3 = state._raw[basis_3]; - Complex res0 = matrix.val[0][0] * val0 + matrix.val[0][1] * val1 + - matrix.val[0][2] * val2 + matrix.val[0][3] * val3; - Complex res1 = matrix.val[1][0] * val0 + matrix.val[1][1] * val1 + - matrix.val[1][2] * val2 + matrix.val[1][3] * val3; - Complex res2 = matrix.val[2][0] * val0 + matrix.val[2][1] * val1 + - matrix.val[2][2] * val2 + matrix.val[2][3] * val3; - Complex res3 = matrix.val[3][0] * val0 + matrix.val[3][1] * val1 + - matrix.val[3][2] * val2 + matrix.val[3][3] * val3; + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1 + matrix[0][2] * val2 + + matrix[0][3] * val3; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1 + matrix[1][2] * val2 + + matrix[1][3] * val3; + Complex res2 = matrix[2][0] * val0 + matrix[2][1] * val1 + matrix[2][2] * val2 + + matrix[2][3] * val3; + Complex res3 = matrix[3][0] * val0 + matrix[3][1] * val1 + matrix[3][2] * val2 + + matrix[3][3] * val3; state._raw[basis_0] = res0; state._raw[basis_1] = res1; state._raw[basis_2] = res2; diff --git a/scaluq/gate/update_ops_pauli.cpp b/scaluq/gate/update_ops_pauli.cpp deleted file mode 100644 index 58695a21..00000000 --- a/scaluq/gate/update_ops_pauli.cpp +++ /dev/null @@ -1,72 +0,0 @@ -#include -#include - -#include "../constant.hpp" -#include "../operator/pauli_operator.hpp" -#include "../types.hpp" -#include "../util/utility.hpp" -#include "update_ops.hpp" - -namespace scaluq { -namespace internal { - -// まだ -void pauli_gate(std::uint64_t control_mask, const PauliOperator& pauli, StateVector& state) { - pauli.apply_to_state(state); -} - -// まだ -void pauli_rotation_gate(std::uint64_t control_mask, - const PauliOperator& pauli, - double angle, - StateVector& state) { - auto [bit_flip_mask_vector, phase_flip_mask_vector] = pauli.get_XZ_mask_representation(); - std::uint64_t bit_flip_mask = internal::BitVector(bit_flip_mask_vector).data_raw()[0]; - std::uint64_t phase_flip_mask = internal::BitVector(phase_flip_mask_vector).data_raw()[0]; - std::uint64_t global_phase_90_rot_count = std::popcount(bit_flip_mask & phase_flip_mask); - Complex true_angle = angle * pauli.get_coef(); - const Complex cosval = Kokkos::cos(-true_angle / 2); - const Complex sinval = Kokkos::sin(-true_angle / 2); - if (bit_flip_mask == 0) { - const Complex cval_min = cosval - Complex(0, 1) * sinval; - const Complex cval_pls = cosval + Complex(0, 1) * sinval; - Kokkos::parallel_for( - state.dim(), KOKKOS_LAMBDA(std::uint64_t state_idx) { - if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { - state._raw[state_idx] *= cval_min; - } else { - state._raw[state_idx] *= cval_pls; - } - }); - Kokkos::fence(); - return; - } else { - const std::uint64_t insert_idx = internal::BitVector(bit_flip_mask_vector).msb(); - Kokkos::parallel_for( - state.dim() >> 1, KOKKOS_LAMBDA(std::uint64_t state_idx) { - std::uint64_t basis_0 = internal::insert_zero_to_basis_index(state_idx, insert_idx); - std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; - - int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1; - int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) & 1; - - // fetch values - Complex cval_0 = state._raw[basis_0]; - Complex cval_1 = state._raw[basis_1]; - - // set values - state._raw[basis_0] = - cosval * cval_0 + - Complex(0, 1) * sinval * cval_1 * - PHASE_M90ROT().val[(global_phase_90_rot_count + bit_parity_0 * 2) % 4]; - state._raw[basis_1] = - cosval * cval_1 + - Complex(0, 1) * sinval * cval_0 * - PHASE_M90ROT().val[(global_phase_90_rot_count + bit_parity_1 * 2) % 4]; - }); - Kokkos::fence(); - } -} - -} // namespace internal -} // namespace scaluq diff --git a/scaluq/gate/update_ops_standard.cpp b/scaluq/gate/update_ops_standard.cpp index 123bd36c..1e64854b 100644 --- a/scaluq/gate/update_ops_standard.cpp +++ b/scaluq/gate/update_ops_standard.cpp @@ -118,7 +118,7 @@ void rx_gate(std::uint64_t target_mask, StateVector& state) { const double cosval = std::cos(angle / 2.); const double sinval = std::sin(angle / 2.); - matrix_2_2 matrix = {cosval, Complex(0, -sinval), Complex(0, -sinval), cosval}; + Matrix2x2 matrix = {cosval, Complex(0, -sinval), Complex(0, -sinval), cosval}; one_target_dense_matrix_gate(target_mask, control_mask, matrix, state); } @@ -128,20 +128,20 @@ void ry_gate(std::uint64_t target_mask, StateVector& state) { const double cosval = std::cos(angle / 2.); const double sinval = std::sin(angle / 2.); - matrix_2_2 matrix = {cosval, -sinval, sinval, cosval}; + Matrix2x2 matrix = {cosval, -sinval, sinval, cosval}; one_target_dense_matrix_gate(target_mask, control_mask, matrix, state); } void one_target_diagonal_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, - const diagonal_matrix_2_2& diag, + const DiagonalMatrix2x2& diag, StateVector& state) { Kokkos::parallel_for( state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) { std::uint64_t basis = insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; - state._raw[basis] *= diag.val[0]; - state._raw[basis | target_mask] *= diag.val[1]; + state._raw[basis] *= diag[0]; + state._raw[basis | target_mask] *= diag[1]; }); Kokkos::fence(); } @@ -152,11 +152,11 @@ void rz_gate(std::uint64_t target_mask, StateVector& state) { const double cosval = std::cos(angle / 2.); const double sinval = std::sin(angle / 2.); - diagonal_matrix_2_2 diag = {Complex(cosval, -sinval), Complex(cosval, sinval)}; + DiagonalMatrix2x2 diag = {Complex(cosval, -sinval), Complex(cosval, sinval)}; one_target_diagonal_matrix_gate(target_mask, control_mask, diag, state); } -matrix_2_2 get_IBMQ_matrix(double theta, double phi, double lambda) { +Matrix2x2 get_IBMQ_matrix(double theta, double phi, double lambda) { Complex exp_val1 = Kokkos::exp(Complex(0, phi)); Complex exp_val2 = Kokkos::exp(Complex(0, lambda)); Complex cos_val = Kokkos::cos(theta / 2.); @@ -186,7 +186,7 @@ void u2_gate(std::uint64_t target_mask, double lambda, StateVector& state) { one_target_dense_matrix_gate( - target_mask, control_mask, get_IBMQ_matrix(PI() / 2., phi, lambda), state); + target_mask, control_mask, get_IBMQ_matrix(Kokkos::numbers::pi / 2., phi, lambda), state); } void u3_gate(std::uint64_t target_mask, diff --git a/scaluq/operator/apply_pauli.cpp b/scaluq/operator/apply_pauli.cpp new file mode 100644 index 00000000..27549754 --- /dev/null +++ b/scaluq/operator/apply_pauli.cpp @@ -0,0 +1,101 @@ +#include "apply_pauli.hpp" + +#include + +#include "../constant.hpp" +#include "../types.hpp" +#include "../util/utility.hpp" + +namespace scaluq::internal { +void apply_pauli(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + StateVector& state_vector) { + if (bit_flip_mask == 0) { + Kokkos::parallel_for( + state_vector.dim() >> std::popcount(control_mask), KOKKOS_LAMBDA(std::uint64_t i) { + std::uint64_t state_idx = + insert_zero_at_mask_positions(i, control_mask) | control_mask; + if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { + state_vector._raw[state_idx] *= -coef; + } else { + state_vector._raw[state_idx] *= coef; + } + }); + Kokkos::fence(); + return; + } + std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; + std::uint64_t global_phase_90rot_count = std::popcount(bit_flip_mask & phase_flip_mask); + Complex global_phase = PHASE_M90ROT()[global_phase_90rot_count % 4]; + Kokkos::parallel_for( + state_vector.dim() >> (std::popcount(control_mask) + 1), KOKKOS_LAMBDA(std::uint64_t i) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) | control_mask; + std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; + Complex tmp1 = state_vector._raw[basis_0] * global_phase; + Complex tmp2 = state_vector._raw[basis_1] * global_phase; + if (Kokkos::popcount(basis_0 & phase_flip_mask) & 1) tmp2 = -tmp2; + if (Kokkos::popcount(basis_1 & phase_flip_mask) & 1) tmp1 = -tmp1; + state_vector._raw[basis_0] = tmp2 * coef; + state_vector._raw[basis_1] = tmp1 * coef; + }); + Kokkos::fence(); +} +void apply_pauli_rotation(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + double angle, + StateVector& state_vector) { + std::uint64_t global_phase_90_rot_count = std::popcount(bit_flip_mask & phase_flip_mask); + Complex true_angle = angle * coef; + const Complex cosval = Kokkos::cos(-true_angle / 2); + const Complex sinval = Kokkos::sin(-true_angle / 2); + if (bit_flip_mask == 0) { + const Complex cval_min = cosval - Complex(0, 1) * sinval; + const Complex cval_pls = cosval + Complex(0, 1) * sinval; + Kokkos::parallel_for( + state_vector.dim() >> std::popcount(control_mask), KOKKOS_LAMBDA(std::uint64_t i) { + std::uint64_t state_idx = + insert_zero_at_mask_positions(i, control_mask) | control_mask; + if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { + state_vector._raw[state_idx] *= cval_min; + } else { + state_vector._raw[state_idx] *= cval_pls; + } + }); + Kokkos::fence(); + return; + } else { + std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; + Kokkos::parallel_for( + state_vector.dim() >> (std::popcount(control_mask) + 1), + KOKKOS_LAMBDA(std::uint64_t i) { + std::uint64_t basis_0 = + internal::insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) | + control_mask; + std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; + + int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1; + int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) & 1; + + // fetch values + Complex cval_0 = state_vector._raw[basis_0]; + Complex cval_1 = state_vector._raw[basis_1]; + + // set values + state_vector._raw[basis_0] = + cosval * cval_0 + + Complex(0, 1) * sinval * cval_1 * + PHASE_M90ROT()[(global_phase_90_rot_count + bit_parity_0 * 2) % 4]; + state_vector._raw[basis_1] = + cosval * cval_1 + + Complex(0, 1) * sinval * cval_0 * + PHASE_M90ROT()[(global_phase_90_rot_count + bit_parity_1 * 2) % 4]; + }); + Kokkos::fence(); + } +} +} // namespace scaluq::internal diff --git a/scaluq/operator/apply_pauli.hpp b/scaluq/operator/apply_pauli.hpp new file mode 100644 index 00000000..ceb41ab5 --- /dev/null +++ b/scaluq/operator/apply_pauli.hpp @@ -0,0 +1,17 @@ +#pragma once + +#include "../state/state_vector.hpp" + +namespace scaluq::internal { +void apply_pauli(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + StateVector& state_vector); +void apply_pauli_rotation(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + double angle, + StateVector& state_vector); +} // namespace scaluq::internal diff --git a/scaluq/operator/operator.cpp b/scaluq/operator/operator.cpp index ee17d8d2..42883295 100644 --- a/scaluq/operator/operator.cpp +++ b/scaluq/operator/operator.cpp @@ -11,7 +11,7 @@ Operator::Operator(std::uint64_t n_qubits) : _n_qubits(n_qubits) {} std::string Operator::to_string() const { std::stringstream ss; for (auto itr = _terms.begin(); itr != _terms.end(); ++itr) { - ss << itr->get_coef() << " " << itr->get_pauli_string(); + ss << itr->coef() << " " << itr->get_pauli_string(); if (itr != prev(_terms.end())) { ss << " + "; } @@ -21,9 +21,9 @@ std::string Operator::to_string() const { void Operator::add_operator(const PauliOperator& mpt) { add_operator(PauliOperator{mpt}); } void Operator::add_operator(PauliOperator&& mpt) { - _is_hermitian &= mpt.get_coef().imag() == 0.; + _is_hermitian &= mpt.coef().imag() == 0.; if (![&] { - const auto& target_list = mpt.get_target_qubit_list(); + const auto& target_list = mpt.target_qubit_list(); if (target_list.empty()) return true; return *std::max_element(target_list.begin(), target_list.end()) < _n_qubits; }()) { @@ -47,9 +47,9 @@ void Operator::add_random_operator(std::uint64_t operator_count, std::uint64_t s } void Operator::optimize() { - std::map, Complex> pauli_and_coef; + std::map, Complex> pauli_and_coef; for (const auto& pauli : _terms) { - pauli_and_coef[pauli.get_XZ_mask_representation()] += pauli.get_coef(); + pauli_and_coef[pauli.get_XZ_mask_representation()] += pauli.coef(); } _terms.clear(); for (const auto& [mask, coef] : pauli_and_coef) { @@ -92,12 +92,12 @@ Complex Operator::get_expectation_value(const StateVector& state_vector) const { Kokkos::DefaultHostExecutionSpace(), terms_view, bmasks_host, - [](const PauliOperator& pauli) { return pauli._ptr->_bit_flip_mask.data_raw()[0]; }); + [](const PauliOperator& pauli) { return pauli._ptr->_bit_flip_mask; }); Kokkos::Experimental::transform( Kokkos::DefaultHostExecutionSpace(), terms_view, pmasks_host, - [](const PauliOperator& pauli) { return pauli._ptr->_phase_flip_mask.data_raw()[0]; }); + [](const PauliOperator& pauli) { return pauli._ptr->_phase_flip_mask; }); Kokkos::Experimental::transform(Kokkos::DefaultHostExecutionSpace(), terms_view, coefs_host, @@ -133,7 +133,7 @@ Complex Operator::get_expectation_value(const StateVector& state_vector) const { sizeof(std::uint64_t) * 8 - Kokkos::countl_zero(bit_flip_mask) - 1; std::uint64_t global_phase_90rot_count = Kokkos::popcount(bit_flip_mask & phase_flip_mask); - Complex global_phase = PHASE_90ROT().val[global_phase_90rot_count % 4]; + Complex global_phase = internal::PHASE_90ROT()[global_phase_90rot_count % 4]; std::uint64_t basis_0 = internal::insert_zero_to_basis_index(state_idx, pivot); std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; double tmp = @@ -166,11 +166,11 @@ Complex Operator::get_transition_amplitude(const StateVector& state_vector_bra, std::vector coefs_vector(nterms); std::transform( _terms.begin(), _terms.end(), bmasks_vector.begin(), [](const PauliOperator& pauli) { - return pauli._ptr->_bit_flip_mask.data_raw()[0]; + return pauli._ptr->_bit_flip_mask; }); std::transform( _terms.begin(), _terms.end(), pmasks_vector.begin(), [](const PauliOperator& pauli) { - return pauli._ptr->_phase_flip_mask.data_raw()[0]; + return pauli._ptr->_phase_flip_mask; }); std::transform( _terms.begin(), _terms.end(), coefs_vector.begin(), [](const PauliOperator& pauli) { @@ -204,7 +204,7 @@ Complex Operator::get_transition_amplitude(const StateVector& state_vector_bra, sizeof(std::uint64_t) * 8 - Kokkos::countl_zero(bit_flip_mask) - 1; std::uint64_t global_phase_90rot_count = Kokkos::popcount(bit_flip_mask & phase_flip_mask); - Complex global_phase = PHASE_90ROT().val[global_phase_90rot_count % 4]; + Complex global_phase = internal::PHASE_90ROT()[global_phase_90rot_count % 4]; std::uint64_t basis_0 = internal::insert_zero_to_basis_index(state_idx, pivot); std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; Complex tmp1 = Kokkos::conj(state_vector_bra._raw[basis_1]) * diff --git a/scaluq/operator/operator.hpp b/scaluq/operator/operator.hpp index cc71abf7..2c472dda 100644 --- a/scaluq/operator/operator.hpp +++ b/scaluq/operator/operator.hpp @@ -67,4 +67,59 @@ class Operator { std::uint64_t _n_qubits; bool _is_hermitian = true; }; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_operator_operator_hpp(nb::module_& m) { + nb::class_(m, "Operator", "General quantum operator class.") + .def(nb::init(), + "qubit_count"_a, + "Initialize operator with specified number of qubits.") + .def("is_hermitian", &Operator::is_hermitian, "Check if the operator is Hermitian.") + .def("n_qubits", &Operator::n_qubits, "Get the number of qubits the operator acts on.") + .def("terms", &Operator::terms, "Get the list of Pauli terms that make up the operator.") + .def("to_string", &Operator::to_string, "Get string representation of the operator.") + .def("add_operator", + nb::overload_cast(&Operator::add_operator), + "Add a Pauli operator to this operator.") + .def( + "add_random_operator", + [](Operator& op, std::uint64_t operator_count, std::optional seed) { + return op.add_random_operator(operator_count, + seed.value_or(std::random_device{}())); + }, + "operator_count"_a, + "seed"_a = std::nullopt, + "Add a specified number of random Pauli operators to this operator. An optional seed " + "can be provided for reproducibility.") + .def("optimize", &Operator::optimize, "Optimize the operator by combining like terms.") + .def("get_dagger", + &Operator::get_dagger, + "Get the adjoint (Hermitian conjugate) of the operator.") + .def("apply_to_state", &Operator::apply_to_state, "Apply the operator to a state vector.") + .def("get_expectation_value", + &Operator::get_expectation_value, + "Get the expectation value of the operator with respect to a state vector.") + .def("get_transition_amplitude", + &Operator::get_transition_amplitude, + "Get the transition amplitude of the operator between two state vectors.") + .def(nb::self *= Complex()) + .def(nb::self * Complex()) + .def(+nb::self) + .def(-nb::self) + .def(nb::self += nb::self) + .def(nb::self + nb::self) + .def(nb::self -= nb::self) + .def(nb::self - nb::self) + .def(nb::self * nb::self) + .def(nb::self *= nb::self) + .def(nb::self += PauliOperator()) + .def(nb::self + PauliOperator()) + .def(nb::self -= PauliOperator()) + .def(nb::self - PauliOperator()) + .def(nb::self *= PauliOperator()) + .def(nb::self * PauliOperator()); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/operator/pauli_operator.cpp b/scaluq/operator/pauli_operator.cpp index 65815487..14406fb6 100644 --- a/scaluq/operator/pauli_operator.cpp +++ b/scaluq/operator/pauli_operator.cpp @@ -6,7 +6,8 @@ namespace scaluq { PauliOperator::Data::Data(Complex coef) : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) {} -PauliOperator::Data::Data(std::string_view pauli_string, Complex coef) : _coef(coef) { +PauliOperator::Data::Data(std::string_view pauli_string, Complex coef) + : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) { auto ss = std::stringstream(std::string(pauli_string)); while (1) { char pauli; @@ -29,7 +30,7 @@ PauliOperator::Data::Data(std::string_view pauli_string, Complex coef) : _coef(c } PauliOperator::Data::Data(const std::vector& pauli_id_par_qubit, Complex coef) - : _coef(coef) { + : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) { for (std::uint64_t i = 0; i < pauli_id_par_qubit.size(); ++i) { add_single_pauli(i, pauli_id_par_qubit[i]); } @@ -38,7 +39,7 @@ PauliOperator::Data::Data(const std::vector& pauli_id_par_qubit, PauliOperator::Data::Data(const std::vector& target_qubit_list, const std::vector& pauli_id_list, Complex coef) - : _coef(coef) { + : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) { if (target_qubit_list.size() != pauli_id_list.size()) { throw std::runtime_error( "PauliOperator::PauliOperator: target_qubit_list must have same size to pauli_id_list"); @@ -48,48 +49,46 @@ PauliOperator::Data::Data(const std::vector& target_qubit_list, } } -PauliOperator::Data::Data(const std::vector& bit_flip_mask, - const std::vector& phase_flip_mask, - Complex coef) - : _coef(coef) { - std::uint64_t num_y = 0; - std::uint64_t max_target = 0; - if (auto msb = internal::BitVector(bit_flip_mask).msb(); - msb != std::numeric_limits::max() && max_target < msb) - max_target = msb; - if (auto msb = internal::BitVector(phase_flip_mask).msb(); - msb != std::numeric_limits::max() && max_target < msb) - max_target = msb; - for (std::uint64_t target_idx = 0; target_idx <= max_target; target_idx++) { - if (!bit_flip_mask[target_idx]) { - if (!phase_flip_mask[target_idx]) +PauliOperator::Data::Data(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex coef) + : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) { + for (std::uint64_t target_idx = 0; target_idx < sizeof(std::uint64_t) * 8; target_idx++) { + bool bit_flip = bit_flip_mask >> target_idx & 1; + bool phase_flip = phase_flip_mask >> target_idx & 1; + if (!bit_flip) { + if (!phase_flip) continue; else add_single_pauli(target_idx, 3); } else { - if (!phase_flip_mask[target_idx]) + if (!phase_flip) add_single_pauli(target_idx, 1); else { add_single_pauli(target_idx, 2); - ++num_y; } } } } void PauliOperator::Data::add_single_pauli(std::uint64_t target_qubit, std::uint64_t pauli_id) { + if (target_qubit >= sizeof(std::uint64_t) * 8) { + throw std::runtime_error( + "PauliOperator::Data::add_single_pauli: target_qubit is too large"); + } + if (pauli_id >= 4) { + throw std::runtime_error("PauliOperator::Data::add_single_pauli: pauli_id is invalid"); + } _target_qubit_list.push_back(target_qubit); _pauli_id_list.push_back(pauli_id); - if ((_bit_flip_mask | _phase_flip_mask)[target_qubit]) { + if ((_bit_flip_mask | _phase_flip_mask) >> target_qubit & 1) { throw std::runtime_error( "PauliOperator::Data::add_single_pauli: You cannot add single pauli twice for same " "qubit."); } if (pauli_id == PauliOperator::X || pauli_id == PauliOperator::Y) { - _bit_flip_mask[target_qubit] = true; + _bit_flip_mask |= 1ULL << target_qubit; } if (pauli_id == PauliOperator::Y || pauli_id == PauliOperator::Z) { - _phase_flip_mask[target_qubit] = true; + _phase_flip_mask |= 1ULL << target_qubit; } } @@ -108,52 +107,14 @@ std::string PauliOperator::get_pauli_string() const { return res; } -void PauliOperator::apply_to_state(StateVector& state_vector) const { - if (state_vector.n_qubits() < get_qubit_count()) { - throw std::runtime_error( - "PauliOperator::apply_to_state: n_qubits of state_vector is too small to apply the " - "operator"); - } - std::uint64_t bit_flip_mask = _ptr->_bit_flip_mask.data_raw()[0]; - std::uint64_t phase_flip_mask = _ptr->_phase_flip_mask.data_raw()[0]; - Complex coef = get_coef(); - if (bit_flip_mask == 0) { - Kokkos::parallel_for( - state_vector.dim(), KOKKOS_LAMBDA(std::uint64_t state_idx) { - if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { - state_vector._raw[state_idx] *= -coef; - } else { - state_vector._raw[state_idx] *= coef; - } - }); - Kokkos::fence(); - return; - } - std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; - std::uint64_t global_phase_90rot_count = std::popcount(bit_flip_mask & phase_flip_mask); - Complex global_phase = PHASE_M90ROT().val[global_phase_90rot_count % 4]; - Kokkos::parallel_for( - state_vector.dim() >> 1, KOKKOS_LAMBDA(std::uint64_t state_idx) { - std::uint64_t basis_0 = internal::insert_zero_to_basis_index(state_idx, pivot); - std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; - Complex tmp1 = state_vector._raw[basis_0] * global_phase; - Complex tmp2 = state_vector._raw[basis_1] * global_phase; - if (Kokkos::popcount(basis_0 & phase_flip_mask) & 1) tmp2 = -tmp2; - if (Kokkos::popcount(basis_1 & phase_flip_mask) & 1) tmp1 = -tmp1; - state_vector._raw[basis_0] = tmp2 * coef; - state_vector._raw[basis_1] = tmp1 * coef; - }); - Kokkos::fence(); -} - Complex PauliOperator::get_expectation_value(const StateVector& state_vector) const { if (state_vector.n_qubits() < get_qubit_count()) { throw std::runtime_error( "PauliOperator::get_expectation_value: n_qubits of state_vector is too small to apply " "the operator"); } - std::uint64_t bit_flip_mask = _ptr->_bit_flip_mask.data_raw()[0]; - std::uint64_t phase_flip_mask = _ptr->_phase_flip_mask.data_raw()[0]; + std::uint64_t bit_flip_mask = _ptr->_bit_flip_mask; + std::uint64_t phase_flip_mask = _ptr->_phase_flip_mask; if (bit_flip_mask == 0) { double res; Kokkos::parallel_reduce( @@ -170,7 +131,7 @@ Complex PauliOperator::get_expectation_value(const StateVector& state_vector) co } std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; std::uint64_t global_phase_90rot_count = std::popcount(bit_flip_mask & phase_flip_mask); - Complex global_phase = PHASE_90ROT().val[global_phase_90rot_count % 4]; + Complex global_phase = internal::PHASE_90ROT()[global_phase_90rot_count % 4]; double res; Kokkos::parallel_reduce( state_vector.dim() >> 1, @@ -196,8 +157,8 @@ Complex PauliOperator::get_transition_amplitude(const StateVector& state_vector_ "PauliOperator::get_expectation_value: n_qubits of state_vector is too small to apply " "the operator"); } - std::uint64_t bit_flip_mask = _ptr->_bit_flip_mask.data_raw()[0]; - std::uint64_t phase_flip_mask = _ptr->_phase_flip_mask.data_raw()[0]; + std::uint64_t bit_flip_mask = _ptr->_bit_flip_mask; + std::uint64_t phase_flip_mask = _ptr->_phase_flip_mask; if (bit_flip_mask == 0) { Complex res; Kokkos::parallel_reduce( @@ -214,7 +175,7 @@ Complex PauliOperator::get_transition_amplitude(const StateVector& state_vector_ } std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; std::uint64_t global_phase_90rot_count = std::popcount(bit_flip_mask & phase_flip_mask); - Complex global_phase = PHASE_90ROT().val[global_phase_90rot_count % 4]; + Complex global_phase = internal::PHASE_90ROT()[global_phase_90rot_count % 4]; Complex res; Kokkos::parallel_reduce( state_vector_bra.dim() >> 1, @@ -234,7 +195,7 @@ Complex PauliOperator::get_transition_amplitude(const StateVector& state_vector_ return _ptr->_coef * res; } -[[nodiscard]] ComplexMatrix PauliOperator::get_matrix_ignoring_coef() const { +[[nodiscard]] internal::ComplexMatrix PauliOperator::get_matrix_ignoring_coef() const { std::uint64_t flip_mask, phase_mask, rot90_count; Kokkos::parallel_reduce( Kokkos::RangePolicy(0, _ptr->_pauli_id_list.size()), @@ -258,36 +219,37 @@ Complex PauliOperator::get_transition_amplitude(const StateVector& state_vector_ rot90_count); std::vector rot = {1, -1.i, -1, 1.i}; std::uint64_t matrix_dim = 1ULL << _ptr->_pauli_id_list.size(); - ComplexMatrix mat = ComplexMatrix::Zero(matrix_dim, matrix_dim); + internal::ComplexMatrix mat = internal::ComplexMatrix::Zero(matrix_dim, matrix_dim); for (std::uint64_t index = 0; index < matrix_dim; index++) { const StdComplex sign = 1. - 2. * (Kokkos::popcount(index & phase_mask) % 2); mat(index, index ^ flip_mask) = rot[rot90_count % 4] * sign; } return mat; } -[[nodiscard]] ComplexMatrix PauliOperator::get_matrix() const { +[[nodiscard]] internal::ComplexMatrix PauliOperator::get_matrix() const { return get_matrix_ignoring_coef() * StdComplex(_ptr->_coef); } PauliOperator PauliOperator::operator*(const PauliOperator& target) const { int extra_90rot_cnt = 0; - auto x_left = _ptr->_bit_flip_mask - _ptr->_phase_flip_mask; + auto x_left = _ptr->_bit_flip_mask & ~_ptr->_phase_flip_mask; auto y_left = _ptr->_bit_flip_mask & _ptr->_phase_flip_mask; - auto z_left = _ptr->_phase_flip_mask - _ptr->_bit_flip_mask; - auto x_right = target._ptr->_bit_flip_mask - target._ptr->_phase_flip_mask; + auto z_left = _ptr->_phase_flip_mask & ~_ptr->_bit_flip_mask; + auto x_right = target._ptr->_bit_flip_mask & ~target._ptr->_phase_flip_mask; auto y_right = target._ptr->_bit_flip_mask & target._ptr->_phase_flip_mask; - auto z_right = target._ptr->_phase_flip_mask - target._ptr->_bit_flip_mask; - extra_90rot_cnt += (x_left & y_right).popcount(); // XY = iZ - extra_90rot_cnt += (y_left & z_right).popcount(); // YZ = iX - extra_90rot_cnt += (z_left & x_right).popcount(); // ZX = iY - extra_90rot_cnt -= (x_left & z_right).popcount(); // XZ = -iY - extra_90rot_cnt -= (y_left & x_right).popcount(); // YX = -iZ - extra_90rot_cnt -= (z_left & y_right).popcount(); // ZY = -iX + auto z_right = target._ptr->_phase_flip_mask & ~target._ptr->_bit_flip_mask; + extra_90rot_cnt += std::popcount(x_left & y_right); // XY = iZ + extra_90rot_cnt += std::popcount(y_left & z_right); // YZ = iX + extra_90rot_cnt += std::popcount(z_left & x_right); // ZX = iY + extra_90rot_cnt -= std::popcount(x_left & z_right); // XZ = -iY + extra_90rot_cnt -= std::popcount(y_left & x_right); // YX = -iZ + extra_90rot_cnt -= std::popcount(z_left & y_right); // ZY = -iX extra_90rot_cnt %= 4; if (extra_90rot_cnt < 0) extra_90rot_cnt += 4; - return PauliOperator(_ptr->_bit_flip_mask ^ target._ptr->_bit_flip_mask, - _ptr->_phase_flip_mask ^ target._ptr->_phase_flip_mask, - _ptr->_coef * target._ptr->_coef * PHASE_90ROT().val[extra_90rot_cnt]); + return PauliOperator( + _ptr->_bit_flip_mask ^ target._ptr->_bit_flip_mask, + _ptr->_phase_flip_mask ^ target._ptr->_phase_flip_mask, + _ptr->_coef * target._ptr->_coef * internal::PHASE_90ROT()[extra_90rot_cnt]); } } // namespace scaluq diff --git a/scaluq/operator/pauli_operator.hpp b/scaluq/operator/pauli_operator.hpp index c5d04891..83bb229e 100644 --- a/scaluq/operator/pauli_operator.hpp +++ b/scaluq/operator/pauli_operator.hpp @@ -5,7 +5,7 @@ #include "../state/state_vector.hpp" #include "../types.hpp" -#include "../util/bit_vector.hpp" +#include "apply_pauli.hpp" namespace scaluq { class PauliOperator { @@ -17,7 +17,7 @@ class PauliOperator { friend class Operator; std::vector _target_qubit_list, _pauli_id_list; Complex _coef; - internal::BitVector _bit_flip_mask, _phase_flip_mask; + std::uint64_t _bit_flip_mask, _phase_flip_mask; public: explicit Data(Complex coef = 1.); @@ -26,17 +26,13 @@ class PauliOperator { const std::vector& pauli_id_list, Complex coef = 1.); Data(const std::vector& pauli_id_par_qubit, Complex coef = 1.); - Data(const std::vector& bit_flip_mask, - const std::vector& phase_flip_mask, - Complex coef); + Data(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex coef); void add_single_pauli(std::uint64_t target_qubit, std::uint64_t pauli_id); - Complex get_coef() const { return _coef; } + Complex coef() const { return _coef; } void set_coef(Complex c) { _coef = c; } - const std::vector& get_target_qubit_list() const { - return _target_qubit_list; - } - const std::vector& get_pauli_id_list() const { return _pauli_id_list; } - std::tuple, std::vector> get_XZ_mask_representation() const { + const std::vector& target_qubit_list() const { return _target_qubit_list; } + const std::vector& pauli_id_list() const { return _pauli_id_list; } + std::tuple get_XZ_mask_representation() const { return {_bit_flip_mask, _phase_flip_mask}; } }; @@ -57,20 +53,18 @@ class PauliOperator { : _ptr(std::make_shared(target_qubit_list, pauli_id_list, coef)) {} PauliOperator(const std::vector& pauli_id_par_qubit, Complex coef = 1.) : _ptr(std::make_shared(pauli_id_par_qubit, coef)) {} - PauliOperator(const std::vector& bit_flip_mask, - const std::vector& phase_flip_mask, - Complex coef) + PauliOperator(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex coef = 1.) : _ptr(std::make_shared(bit_flip_mask, phase_flip_mask, coef)) {} - [[nodiscard]] inline Complex get_coef() const { return _ptr->get_coef(); } - [[nodiscard]] inline const std::vector& get_target_qubit_list() const { - return _ptr->get_target_qubit_list(); + [[nodiscard]] inline Complex coef() const { return _ptr->coef(); } + [[nodiscard]] inline const std::vector& target_qubit_list() const { + return _ptr->target_qubit_list(); } - [[nodiscard]] inline const std::vector& get_pauli_id_list() const { - return _ptr->get_pauli_id_list(); + [[nodiscard]] inline const std::vector& pauli_id_list() const { + return _ptr->pauli_id_list(); } - [[nodiscard]] inline std::tuple, std::vector> - get_XZ_mask_representation() const { + [[nodiscard]] inline std::tuple get_XZ_mask_representation() + const { return _ptr->get_XZ_mask_representation(); } [[nodiscard]] std::string get_pauli_string() const; @@ -83,13 +77,21 @@ class PauliOperator { return std::ranges::max(_ptr->_target_qubit_list) + 1; } - void apply_to_state(StateVector& state_vector) const; + void apply_to_state(StateVector& state_vector) const { + if (state_vector.n_qubits() < get_qubit_count()) { + throw std::runtime_error( + "PauliOperator::apply_to_state: n_qubits of state_vector is too small to apply the " + "operator"); + } + internal::apply_pauli( + 0ULL, _ptr->_bit_flip_mask, _ptr->_phase_flip_mask, _ptr->_coef, state_vector); + } [[nodiscard]] Complex get_expectation_value(const StateVector& state_vector) const; [[nodiscard]] Complex get_transition_amplitude(const StateVector& state_vector_bra, const StateVector& state_vector_ket) const; - [[nodiscard]] ComplexMatrix get_matrix() const; - [[nodiscard]] ComplexMatrix get_matrix_ignoring_coef() const; + [[nodiscard]] internal::ComplexMatrix get_matrix() const; + [[nodiscard]] internal::ComplexMatrix get_matrix_ignoring_coef() const; [[nodiscard]] PauliOperator operator*(const PauliOperator& target) const; [[nodiscard]] inline PauliOperator operator*(Complex target) const { @@ -97,4 +99,124 @@ class PauliOperator { } }; +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_operator_pauli_operator_hpp(nb::module_& m) { + nb::enum_(m, "PauliID") + .value("I", PauliOperator::I) + .value("X", PauliOperator::X) + .value("Y", PauliOperator::Y) + .value("Z", PauliOperator::Z) + .export_values(); + + nb::class_( + m, "PauliOperatorData", "Internal data structure for PauliOperator.") + .def(nb::init(), "coef"_a = 1., "Initialize data with coefficient.") + .def(nb::init(), + "pauli_string"_a, + "coef"_a = 1., + "Initialize data with pauli string.") + .def(nb::init&, + const std::vector&, + Complex>(), + "target_qubit_list"_a, + "pauli_id_list"_a, + "coef"_a = 1., + "Initialize data with target qubits and pauli ids.") + .def(nb::init&, Complex>(), + "pauli_id_par_qubit"_a, + "coef"_a = 1., + "Initialize data with pauli ids per qubit.") + .def(nb::init(), + "bit_flip_mask"_a, + "phase_flip_mask"_a, + "coef"_a = 1., + "Initialize data with bit flip and phase flip masks.") + .def(nb::init(), + "data"_a, + "Initialize pauli operator from Data object.") + .def("add_single_pauli", + &PauliOperator::Data::add_single_pauli, + "target_qubit"_a, + "pauli_id"_a, + "Add a single pauli operation to the data.") + .def("coef", &PauliOperator::Data::coef, "Get the coefficient of the Pauli operator.") + .def("set_coef", + &PauliOperator::Data::set_coef, + "c"_a, + "Set the coefficient of the Pauli operator.") + .def("target_qubit_list", + &PauliOperator::Data::target_qubit_list, + "Get the list of target qubits.") + .def("pauli_id_list", &PauliOperator::Data::pauli_id_list, "Get the list of Pauli IDs.") + .def("get_XZ_mask_representation", + &PauliOperator::Data::get_XZ_mask_representation, + "Get the X and Z mask representation as a tuple of vectors."); + + nb::class_( + m, + "PauliOperator", + "Pauli operator as coef and tensor product of single pauli for each qubit.") + .def(nb::init(), "coef"_a = 1., "Initialize operator which just multiplying coef.") + .def(nb::init&, + const std::vector&, + Complex>(), + "target_qubit_list"_a, + "pauli_id_list"_a, + "coef"_a = 1., + "Initialize pauli operator. For each `i`, single pauli correspond to " + "`pauli_id_list[i]` is applied to `target_qubit_list`-th qubit.") + .def(nb::init(), + "pauli_string"_a, + "coef"_a = 1., + "Initialize pauli operator. If `pauli_string` is `\"X0Y2\"`, Pauli-X is applied to " + "0-th qubit and Pauli-Y is applied to 2-th qubit. In `pauli_string`, spaces are " + "ignored.") + .def(nb::init&, Complex>(), + "pauli_id_par_qubit"_a, + "coef"_a = 1., + "Initialize pauli operator. For each `i`, single pauli correspond to " + "`paul_id_per_qubit` is applied to `i`-th qubit.") + .def(nb::init(), + "bit_flip_mask"_a, + "phase_flip_mask"_a, + "coef"_a = 1., + "Initialize pauli operator. For each `i`, single pauli applied to `i`-th qubit is got " + "from `i-th` bit of `bit_flip_mask` and `phase_flip_mask` as follows.\n\n.. " + "csv-table::\n\n \"bit_flip\",\"phase_flip\",\"pauli\"\n \"0\",\"0\",\"I\"\n " + "\"0\",\"1\",\"Z\"\n \"1\",\"0\",\"X\"\n \"1\",\"1\",\"Y\"") + .def("coef", &PauliOperator::coef, "Get property `coef`.") + .def("target_qubit_list", + &PauliOperator::target_qubit_list, + "Get qubits to be applied pauli.") + .def("pauli_id_list", + &PauliOperator::pauli_id_list, + "Get pauli id to be applied. The order is correspond to the result of " + "`target_qubit_list`") + .def("get_XZ_mask_representation", + &PauliOperator::get_XZ_mask_representation, + "Get single-pauli property as binary integer representation. See description of " + "`__init__(bit_flip_mask_py: int, phase_flip_mask_py: int, coef: float=1.)` for " + "details.") + .def("get_pauli_string", + &PauliOperator::get_pauli_string, + "Get single-pauli property as string representation. See description of " + "`__init__(pauli_string: str, coef: float=1.)` for details.") + .def("get_dagger", &PauliOperator::get_dagger, "Get adjoint operator.") + .def("get_qubit_count", + &PauliOperator::get_qubit_count, + "Get num of qubits to applied with, when count from 0-th qubit. Subset of $[0, " + "\\mathrm{qubit_count})$ is the target.") + .def("apply_to_state", &PauliOperator::apply_to_state, "Apply pauli to state vector.") + .def("get_expectation_value", + &PauliOperator::get_expectation_value, + "Get expectation value of measuring state vector. $\\bra{\\psi}P\\ket{\\psi}$.") + .def("get_transition_amplitude", + &PauliOperator::get_transition_amplitude, + "Get transition amplitude of measuring state vector. $\\bra{\\chi}P\\ket{\\psi}$.") + .def(nb::self * nb::self) + .def(nb::self * Complex()); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/state/state_vector.cpp b/scaluq/state/state_vector.cpp index d82e9c83..2913ba21 100644 --- a/scaluq/state/state_vector.cpp +++ b/scaluq/state/state_vector.cpp @@ -12,13 +12,13 @@ StateVector::StateVector(std::uint64_t n_qubits) set_zero_state(); } -void StateVector::set_amplitude_at_index(std::uint64_t index, const Complex& c) { +void StateVector::set_amplitude_at(std::uint64_t index, const Complex& c) { Kokkos::View host_view("single_value"); host_view() = c; Kokkos::deep_copy(Kokkos::subview(_raw, index), host_view()); } -Complex StateVector::get_amplitude_at_index(std::uint64_t index) const { +Complex StateVector::get_amplitude_at(std::uint64_t index) const { Kokkos::View host_view("single_value"); Kokkos::deep_copy(host_view, Kokkos::subview(_raw, index)); return host_view(); @@ -26,7 +26,7 @@ Complex StateVector::get_amplitude_at_index(std::uint64_t index) const { void StateVector::set_zero_state() { Kokkos::deep_copy(_raw, 0); - set_amplitude_at_index(0, 1); + set_amplitude_at(0, 1); } void StateVector::set_zero_norm_state() { Kokkos::deep_copy(_raw, 0); } @@ -39,7 +39,7 @@ void StateVector::set_computational_basis(std::uint64_t basis) { "computational basis must be smaller than 2^qubit_count"); } Kokkos::deep_copy(_raw, 0); - set_amplitude_at_index(basis, 1); + set_amplitude_at(basis, 1); } StateVector StateVector::Haar_random_state(std::uint64_t n_qubits, std::uint64_t seed) { @@ -60,7 +60,7 @@ std::uint64_t StateVector::n_qubits() const { return this->_n_qubits; } std::uint64_t StateVector::dim() const { return this->_dim; } -std::vector StateVector::amplitudes() const { +std::vector StateVector::get_amplitudes() const { return internal::convert_device_view_to_host_vector(_raw); } @@ -219,7 +219,7 @@ std::vector StateVector::sampling(std::uint64_t sampling_count, std::string StateVector::to_string() const { std::stringstream os; - auto amp = this->amplitudes(); + auto amp = this->get_amplitudes(); os << " *** Quantum State ***\n"; os << " * Qubit Count : " << _n_qubits << '\n'; os << " * Dimension : " << _dim << '\n'; @@ -233,7 +233,7 @@ std::string StateVector::to_string() const { } return tmp; }(i, _n_qubits) - << ": " << amp[i] << std::endl; + << ": " << amp[i] << (i < _dim - 1 ? "\n" : ""); } return os.str(); } diff --git a/scaluq/state/state_vector.hpp b/scaluq/state/state_vector.hpp index 3ef21324..6ede5460 100644 --- a/scaluq/state/state_vector.hpp +++ b/scaluq/state/state_vector.hpp @@ -15,7 +15,7 @@ class StateVector { public: static constexpr std::uint64_t UNMEASURED = 2; - StateVectorView _raw; + Kokkos::View _raw; StateVector() = default; StateVector(std::uint64_t n_qubits); StateVector(const StateVector& other) = default; @@ -25,12 +25,12 @@ class StateVector { /** * @attention Very slow. You should use load() instead if you can. */ - void set_amplitude_at_index(std::uint64_t index, const Complex& c); + void set_amplitude_at(std::uint64_t index, const Complex& c); /** - * @attention Very slow. You should use amplitudes() instead if you can. + * @attention Very slow. You should use get_amplitudes() instead if you can. */ - [[nodiscard]] Complex get_amplitude_at_index(std::uint64_t index) const; + [[nodiscard]] Complex get_amplitude_at(std::uint64_t index) const; [[nodiscard]] static StateVector Haar_random_state(std::uint64_t n_qubits, std::uint64_t seed = std::random_device()()); @@ -46,7 +46,7 @@ class StateVector { [[nodiscard]] std::uint64_t dim() const; - [[nodiscard]] std::vector amplitudes() const; + [[nodiscard]] std::vector get_amplitudes() const; [[nodiscard]] double get_squared_norm() const; @@ -72,4 +72,99 @@ class StateVector { [[nodiscard]] std::string to_string() const; }; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_state_state_vector_hpp(nb::module_& m) { + nb::class_(m, + "StateVector", + "Vector representation of quantum state.\n\n.. note:: Qubit index is " + "start from 0. If the amplitudes of $\\ket{b_{n-1}\\dots b_0}$ is " + "$b_i$, the state is $\\sum_i b_i 2^i$.") + .def(nb::init(), + "Construct state vector with specified qubits, initialized with computational " + "basis $\\ket{0\\dots0}$.") + .def(nb::init(), "Constructing state vector by copying other state.") + .def_static( + "Haar_random_state", + [](std::uint64_t n_qubits, std::optional seed) { + return StateVector::Haar_random_state(n_qubits, + seed.value_or(std::random_device{}())); + }, + "n_qubits"_a, + "seed"_a = std::nullopt, + "Constructing state vector with Haar random state. If seed is not specified, the value " + "from random device is used.") + .def("set_amplitude_at", + &StateVector::set_amplitude_at, + "Manually set amplitude at one index.") + .def("get_amplitude_at", + &StateVector::get_amplitude_at, + "Get amplitude at one index.\n\n.. note:: If you want to get all amplitudes, you " + "should " + "use `StateVector::get_amplitudes()`.") + .def("set_zero_state", + &StateVector::set_zero_state, + "Initialize with computational basis $\\ket{00\\dots0}$.") + .def("set_zero_norm_state", + &StateVector::set_zero_norm_state, + "Initialize with 0 (null vector).") + .def("set_computational_basis", + &StateVector::set_computational_basis, + "Initialize with computational basis \\ket{\\mathrm{basis}}.") + .def("amplitudes", + &StateVector::get_amplitudes, + "Get all amplitudes with as `list[complex]`.") + .def("n_qubits", &StateVector::n_qubits, "Get num of qubits.") + .def("dim", &StateVector::dim, "Get dimension of the vector ($=2^\\mathrm{n\\_qubits}$).") + .def("get_squared_norm", + &StateVector::get_squared_norm, + "Get squared norm of the state. $\\braket{\\psi|\\psi}$.") + .def("normalize", + &StateVector::normalize, + "Normalize state (let $\\braket{\\psi|\\psi} = 1$ by multiplying coef).") + .def("get_zero_probability", + &StateVector::get_zero_probability, + "Get the probability to observe $\\ket{0}$ at specified index.") + .def("get_marginal_probability", + &StateVector::get_marginal_probability, + "Get the marginal probability to observe as specified. Specify the result as n-length " + "list. `0` and `1` represent the qubit is observed and get the value. `2` represents " + "the qubit is not observed.") + .def("get_entropy", &StateVector::get_entropy, "Get the entropy of the vector.") + .def("add_state_vector", + &StateVector::add_state_vector, + "Add other state vector and make superposition. $\\ket{\\mathrm{this}} " + "\\leftarrow " + "\\ket{\\mathrm{this}} + \\ket{\\mathrm{state}}$.") + .def("add_state_vector_with_coef", + &StateVector::add_state_vector_with_coef, + "add other state vector with multiplying the coef and make superposition. " + "$\\ket{\\mathrm{this}}\\leftarrow\\ket{\\mathrm{this}}+\\mathrm{coef}" + "\\ket{\\mathrm{" + "state}}$.") + .def("multiply_coef", + &StateVector::multiply_coef, + "Multiply coef. " + "$\\ket{\\mathrm{this}}\\leftarrow\\mathrm{coef}\\ket{\\mathrm{this}}$.") + .def( + "sampling", + [](const StateVector& state, + std::uint64_t sampling_count, + std::optional seed) { + return state.sampling(sampling_count, seed.value_or(std::random_device{}())); + }, + "sampling_count"_a, + "seed"_a = std::nullopt, + "Sampling specified times. Result is `list[int]` with the `sampling_count` length.") + .def("to_string", &StateVector::to_string, "Information as `str`.") + .def("load", &StateVector::load, "Load amplitudes of `list[int]` with `dim` length.") + .def("__str__", &StateVector::to_string, "Information as `str`.") + .def_ro_static("UNMEASURED", + &StateVector::UNMEASURED, + "Constant used for `StateVector::get_marginal_probability` to express the " + "the qubit is not measured."); +} +} // namespace internal +#endif } // namespace scaluq diff --git a/scaluq/state/state_vector_batched.cpp b/scaluq/state/state_vector_batched.cpp index 9d029f01..03ccc3ee 100644 --- a/scaluq/state/state_vector_batched.cpp +++ b/scaluq/state/state_vector_batched.cpp @@ -7,7 +7,7 @@ StateVectorBatched::StateVectorBatched(std::uint64_t batch_size, std::uint64_t n : _batch_size(batch_size), _n_qubits(n_qubits), _dim(1ULL << _n_qubits), - _raw(StateVectorBatchedView( + _raw(Kokkos::View( Kokkos::ViewAllocateWithoutInitializing("states"), _batch_size, _dim)) { set_zero_state(); } @@ -38,7 +38,7 @@ void StateVectorBatched::set_state_vector(std::uint64_t batch_id, const StateVec Kokkos::fence(); } -StateVector StateVectorBatched::get_state_vector(std::uint64_t batch_id) const { +StateVector StateVectorBatched::get_state_vector_at(std::uint64_t batch_id) const { StateVector ret(_n_qubits); Kokkos::parallel_for( _dim, KOKKOS_CLASS_LAMBDA(std::uint64_t i) { ret._raw(i) = _raw(batch_id, i); }); @@ -133,7 +133,7 @@ std::vector> StateVectorBatched::sampling(std::uint64 Kokkos::DefaultExecutionSpace::array_layout>(result); } -std::vector> StateVectorBatched::amplitudes() const { +std::vector> StateVectorBatched::get_amplitudes() const { return internal::convert_2d_device_view_to_host_vector(_raw); } @@ -367,7 +367,7 @@ std::string StateVectorBatched::to_string() const { for (std::uint64_t b = 0; b < _batch_size; ++b) { StateVector tmp(_n_qubits); os << "--------------------\n"; - os << " * Batch_id : " << b << '\n'; + os << " * Batch id : " << b << '\n'; os << " * State vector : \n"; for (std::uint64_t i = 0; i < _dim; ++i) { os << @@ -378,7 +378,7 @@ std::string StateVectorBatched::to_string() const { } return tmp; }(i, _n_qubits) - << ": " << states_h(b, i) << std::endl; + << ": " << states_h(b, i) << (b < _batch_size - 1 || i < _dim - 1 ? "\n" : ""); } } return os.str(); diff --git a/scaluq/state/state_vector_batched.hpp b/scaluq/state/state_vector_batched.hpp index 6c4c5048..8d8e86eb 100644 --- a/scaluq/state/state_vector_batched.hpp +++ b/scaluq/state/state_vector_batched.hpp @@ -10,7 +10,7 @@ class StateVectorBatched { std::uint64_t _dim; public: - StateVectorBatchedView _raw; + Kokkos::View _raw; StateVectorBatched() = default; StateVectorBatched(std::uint64_t batch_size, std::uint64_t n_qubits); StateVectorBatched(const StateVectorBatched& other) = default; @@ -25,7 +25,7 @@ class StateVectorBatched { void set_state_vector(const StateVector& state); void set_state_vector(std::uint64_t batch_id, const StateVector& state); - [[nodiscard]] StateVector get_state_vector(std::uint64_t batch_id) const; + [[nodiscard]] StateVector get_state_vector_at(std::uint64_t batch_id) const; void set_zero_state(); void set_computational_basis(std::uint64_t basis); @@ -40,7 +40,7 @@ class StateVectorBatched { bool set_same_state, std::uint64_t seed = std::random_device()()); - [[nodiscard]] std::vector> amplitudes() const; + [[nodiscard]] std::vector> get_amplitudes() const; [[nodiscard]] std::vector get_squared_norm() const; @@ -61,4 +61,109 @@ class StateVectorBatched { std::string to_string() const; friend std::ostream& operator<<(std::ostream& os, const StateVectorBatched& states); }; + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_state_state_vector_batched_hpp(nb::module_& m) { + nb::class_( + m, + "StateVectorBatched", + "Batched vector representation of quantum state.\n\n.. note:: Qubit index is start from 0. " + "If the amplitudes of $\\ket{b_{n-1}\\dots b_0}$ is $b_i$, the state is $\\sum_i b_i " + "2^i$.") + .def(nb::init(), + "Construct batched state vector with specified batch size and qubits.") + .def(nb::init(), + "Constructing batched state vector by copying other batched state.") + .def("n_qubits", &StateVectorBatched::n_qubits, "Get num of qubits.") + .def("dim", + &StateVectorBatched::dim, + "Get dimension of the vector ($=2^\\mathrm{n\\_qubits}$).") + .def("batch_size", &StateVectorBatched::batch_size, "Get batch size.") + .def("set_state_vector", + nb::overload_cast(&StateVectorBatched::set_state_vector), + "Set the state vector for all batches.") + .def("set_state_vector", + nb::overload_cast( + &StateVectorBatched::set_state_vector), + "Set the state vector for a specific batch.") + .def("get_state_vector_at", + &StateVectorBatched::get_state_vector_at, + "Get the state vector for a specific batch.") + .def("set_zero_state", + &StateVectorBatched::set_zero_state, + "Initialize all batches with computational basis $\\ket{00\\dots0}$.") + .def("set_zero_norm_state", + &StateVectorBatched::set_zero_norm_state, + "Initialize with 0 (null vector).") + .def("set_computational_basis", + &StateVectorBatched::set_computational_basis, + "Initialize with computational basis \\ket{\\mathrm{basis}}.") + .def( + "sampling", + [](const StateVectorBatched& states, + std::uint64_t sampling_count, + std::optional seed) { + return states.sampling(sampling_count, seed.value_or(std::random_device{}())); + }, + "sampling_count"_a, + "seed"_a = std::nullopt, + "Sampling specified times. Result is `list[list[int]]` with the `sampling_count` " + "length.") + .def_static( + "Haar_random_states", + [](std::uint64_t batch_size, + std::uint64_t n_qubits, + bool set_same_state, + std::optional seed) { + return StateVectorBatched::Haar_random_states( + batch_size, n_qubits, set_same_state, seed.value_or(std::random_device{}())); + }, + "batch_size"_a, + "n_qubits"_a, + "set_same_state"_a, + "seed"_a = std::nullopt, + "Construct batched state vectors with Haar random states. If seed is not " + "specified, the value from random device is used.") + .def("amplitudes", + &StateVectorBatched::get_amplitudes, + "Get all amplitudes with as `list[list[complex]]`.") + .def("get_squared_norm", + &StateVectorBatched::get_squared_norm, + "Get squared norm of each state in the batch. $\\braket{\\psi|\\psi}$.") + .def("normalize", + &StateVectorBatched::normalize, + "Normalize each state in the batch (let $\\braket{\\psi|\\psi} = 1$ by " + "multiplying coef).") + .def("get_zero_probability", + &StateVectorBatched::get_zero_probability, + "Get the probability to observe $\\ket{0}$ at specified index for each state in " + "the batch.") + .def("get_marginal_probability", + &StateVectorBatched::get_marginal_probability, + "Get the marginal probability to observe as specified for each state in the batch. " + "Specify the result as n-length list. `0` and `1` represent the qubit is observed " + "and get the value. `2` represents the qubit is not observed.") + .def("get_entropy", + &StateVectorBatched::get_entropy, + "Get the entropy of each state in the batch.") + .def("add_state_vector", + &StateVectorBatched::add_state_vector, + "Add other batched state vectors and make superposition. $\\ket{\\mathrm{this}} " + "\\leftarrow \\ket{\\mathrm{this}} + \\ket{\\mathrm{states}}$.") + .def("add_state_vector_with_coef", + &StateVectorBatched::add_state_vector_with_coef, + "Add other batched state vectors with multiplying the coef and make superposition. " + "$\\ket{\\mathrm{this}}\\leftarrow\\ket{\\mathrm{this}}+\\mathrm{coef}" + "\\ket{\\mathrm{states}}$.") + .def("load", + &StateVectorBatched::load, + "Load batched amplitudes from `list[list[complex]]`.") + .def("copy", &StateVectorBatched::copy, "Create a copy of the batched state vector.") + .def("to_string", &StateVectorBatched::to_string, "Information as `str`.") + .def("__str__", &StateVectorBatched::to_string, "Information as `str`."); +} +} // namespace internal +#endif + } // namespace scaluq diff --git a/scaluq/types.hpp b/scaluq/types.hpp index a36aa85d..66741688 100644 --- a/scaluq/types.hpp +++ b/scaluq/types.hpp @@ -17,34 +17,38 @@ inline void finalize() { Kokkos::finalize(); } inline bool is_initialized() { return Kokkos::is_initialized(); } inline bool is_finalized() { return Kokkos::is_finalized(); } +using StdComplex = std::complex; using Complex = Kokkos::complex; using namespace std::complex_literals; -using StdComplex = std::complex; +namespace internal { +template +constexpr bool lazy_false_v = false; // Used for lazy evaluation in static_assert. using ComplexMatrix = Eigen::Matrix; using SparseComplexMatrix = Eigen::SparseMatrix; using Matrix = Kokkos::View; -using StateVectorView = Kokkos::View; -using StateVectorBatchedView = Kokkos::View; - -struct array_4 { - Complex val[4]; -}; - -struct matrix_2_2 { - Complex val[2][2]; -}; - -struct matrix_4_4 { - Complex val[4][4]; -}; - -struct diagonal_matrix_2_2 { - Complex val[2]; -}; - +using Matrix2x2 = Kokkos::Array, 2>; +using Matrix4x4 = Kokkos::Array, 4>; +using DiagonalMatrix2x2 = Kokkos::Array; +} // namespace internal + +#ifdef SCALUQ_USE_NANOBIND +namespace internal { +void bind_types_hpp(nb::module_& m) { + m.def("finalize", + &finalize, + "Terminate the Kokkos execution environment. Release the resources.\n\n.. note:: " + "Finalization fails if there exists `StateVector` allocated. You must use " + "`StateVector` only inside inner scopes than the usage of `finalize` or delete all of " + "existing `StateVector`.\n\n.. note:: This is " + "automatically called when the program exits. If you call this manually, you cannot use " + "most of scaluq's functions until the program exits."); + m.def("is_finalized", &is_initialized, "Return true if `finalize()` is already called."); +} +} // namespace internal +#endif struct SparseValue { Complex val; uint32_t r, c; diff --git a/scaluq/util/bit_vector.hpp b/scaluq/util/bit_vector.hpp deleted file mode 100644 index 909f4810..00000000 --- a/scaluq/util/bit_vector.hpp +++ /dev/null @@ -1,178 +0,0 @@ -#pragma once - -#include -#include -#include -#include - -#include "../types.hpp" - -namespace scaluq { -namespace internal { -class BitVector { -public: - constexpr static std::uint64_t BIT_SIZE = sizeof(std::uint64_t) * 8; - - BitVector(std::uint64_t sz = 1) : _data((sz + BIT_SIZE - 1) / BIT_SIZE) {} - BitVector(const std::vector& vec) : _data((vec.size() + BIT_SIZE - 1) / BIT_SIZE) { - for (std::uint64_t i = 0; i < vec.size(); ++i) { - set(i, vec[i]); - } - } - - [[nodiscard]] inline const std::vector& data_raw() const { return _data; } - [[nodiscard]] inline std::vector& data_raw() { return _data; } - - [[nodiscard]] inline bool get(std::uint64_t idx) const { - if (idx >= _data.size() * BIT_SIZE) return false; - return _data[idx / BIT_SIZE] >> (idx % BIT_SIZE) & 1ULL; - } - inline void set(std::uint64_t idx, bool b) { - if (idx >= _data.size() * BIT_SIZE) _data.resize(idx / BIT_SIZE + 1); - if (b) - _data[idx / BIT_SIZE] |= 1ULL << (idx % BIT_SIZE); - else - _data[idx / BIT_SIZE] &= ~(1ULL << (idx % BIT_SIZE)); - } - - template - class _Reference { - friend BitVector; - - public: - _Reference& operator=(bool b) { - static_assert(!Const); - _container.set(_idx, b); - return *this; - } - _Reference& operator&=(bool b) { - static_assert(!Const); - _container.set(_idx, b && _container.get(_idx)); - return *this; - } - _Reference& operator|=(bool b) { - static_assert(!Const); - _container.set(_idx, b || _container.get(_idx)); - return *this; - } - _Reference& operator^=(bool b) { - static_assert(!Const); - _container.set(_idx, b ^ _container.get(_idx)); - return *this; - } - operator bool() const { return _container.get(_idx); } - - private: - using ContainerReference = std::conditional_t; - ContainerReference _container; - const int _idx; - - _Reference(ContainerReference container, int idx) : _container(container), _idx(idx) {} - }; - - using ConstReference = _Reference; - using Reference = _Reference; - - [[nodiscard]] inline ConstReference operator[](int idx) const { - return ConstReference(*this, idx); - } - [[nodiscard]] inline Reference operator[](int idx) { return Reference(*this, idx); } - - inline BitVector& operator&=(const BitVector& rhs) { - if (rhs._data.size() < _data.size()) { - _data.resize(rhs._data.size()); - } - for (std::uint64_t i = 0; i < _data.size(); i++) _data[i] &= rhs._data[i]; - return *this; - } - inline BitVector operator&(const BitVector& rhs) const { return BitVector(*this) &= rhs; } - inline BitVector& operator|=(const BitVector& rhs) { - if (rhs._data.size() > _data.size()) { - _data.resize(rhs._data.size()); - } - for (std::uint64_t i = 0; i < rhs._data.size(); i++) _data[i] |= rhs._data[i]; - return *this; - } - inline BitVector operator|(const BitVector& rhs) const { return BitVector(*this) |= rhs; } - inline BitVector& operator^=(const BitVector& rhs) { - if (rhs._data.size() > _data.size()) { - _data.resize(rhs._data.size()); - } - for (std::uint64_t i = 0; i < rhs._data.size(); i++) _data[i] ^= rhs._data[i]; - return *this; - } - inline BitVector operator^(const BitVector& rhs) const { return BitVector(*this) ^= rhs; } - inline BitVector& operator-=(const BitVector& rhs) { - for (std::uint64_t i = 0; i < std::min(_data.size(), rhs._data.size()); i++) - _data[i] &= ~rhs._data[i]; - return *this; - } - inline BitVector operator-(const BitVector& rhs) const { return BitVector(*this) -= rhs; } - - inline std::weak_ordering operator<=>(const BitVector& other) const { - std::uint64_t sz = std::max(_data.size(), other._data.size()); - for (std::uint64_t i = sz; i-- != 0;) { - std::uint64_t l = i >= _data.size() ? 0ULL : _data[i]; - std::uint64_t r = i >= other._data.size() ? 0ULL : other._data[i]; - if (l != r) return l <=> r; - if (i == 0) break; - } - return std::weak_ordering::equivalent; - } - inline bool operator==(const BitVector& other) const { - std::uint64_t sz = std::max(_data.size(), other._data.size()); - for (std::uint64_t i = sz; i-- != 0;) { - std::uint64_t l = i >= _data.size() ? 0ULL : _data[i]; - std::uint64_t r = i >= other._data.size() ? 0ULL : other._data[i]; - if (l != r) return false; - } - return true; - } - - operator std::vector() const { - std::vector vec(_data.size() * BIT_SIZE); - for (std::uint64_t i = 0; i < vec.size(); ++i) { - vec[i] = get(i); - } - return vec; - } - - inline bool empty() const { - return std::ranges::all_of(_data, [](std::uint64_t x) { return x == 0; }); - } - inline std::uint64_t msb() const { - for (std::uint64_t i = _data.size() - 1; i != std::numeric_limits::max(); - i--) { - if (_data[i] != 0) return (i + 1) * BIT_SIZE - std::countl_zero(_data[i]) - 1; - } - return std::numeric_limits::max(); - } - inline std::uint64_t countr_zero() const { - std::uint64_t res = 0; - for (std::uint64_t i = 0; i < _data.size(); i++) { - std::uint64_t to_add = std::countr_zero(_data[i]); - res += to_add; - if (to_add < BIT_SIZE) break; - } - return res; - } - inline std::uint64_t countr_one() const { - std::uint64_t res = 0; - for (std::uint64_t i = 0; i < _data.size(); i++) { - std::uint64_t to_add = std::countr_one(_data[i]); - res += to_add; - if (to_add < BIT_SIZE) break; - } - return res; - } - inline std::uint64_t popcount() const { - std::uint64_t res = 0; - for (std::uint64_t i = 0; i < _data.size(); i++) res += std::popcount(_data[i]); - return res; - } - -private: - std::vector _data; -}; -} // namespace internal -} // namespace scaluq diff --git a/scaluq/util/utility.hpp b/scaluq/util/utility.hpp index 02268865..7b20f60b 100644 --- a/scaluq/util/utility.hpp +++ b/scaluq/util/utility.hpp @@ -67,16 +67,17 @@ inline std::vector mask_to_vector(std::uint64_t mask) { return indices; } -KOKKOS_INLINE_FUNCTION matrix_2_2 matrix_multiply(const matrix_2_2& matrix1, - const matrix_2_2& matrix2) { - return {matrix1.val[0][0] * matrix2.val[0][0] + matrix1.val[0][1] * matrix2.val[1][0], - matrix1.val[0][0] * matrix2.val[0][1] + matrix1.val[0][1] * matrix2.val[1][1], - matrix1.val[1][0] * matrix2.val[0][0] + matrix1.val[1][1] * matrix2.val[1][0], - matrix1.val[1][0] * matrix2.val[0][1] + matrix1.val[1][1] * matrix2.val[1][1]}; +KOKKOS_INLINE_FUNCTION Matrix2x2 matrix_multiply(const Matrix2x2& matrix1, + const Matrix2x2& matrix2) { + return {matrix1[0][0] * matrix2[0][0] + matrix1[0][1] * matrix2[1][0], + matrix1[0][0] * matrix2[0][1] + matrix1[0][1] * matrix2[1][1], + matrix1[1][0] * matrix2[0][0] + matrix1[1][1] * matrix2[1][0], + matrix1[1][0] * matrix2[0][1] + matrix1[1][1] * matrix2[1][1]}; } -inline ComplexMatrix kronecker_product(const ComplexMatrix& lhs, const ComplexMatrix& rhs) { - ComplexMatrix result(lhs.rows() * rhs.rows(), lhs.cols() * rhs.cols()); +inline internal::ComplexMatrix kronecker_product(const internal::ComplexMatrix& lhs, + const internal::ComplexMatrix& rhs) { + internal::ComplexMatrix result(lhs.rows() * rhs.rows(), lhs.cols() * rhs.cols()); for (int i = 0; i < lhs.rows(); i++) { for (int j = 0; j < lhs.cols(); j++) { result.block(i * rhs.rows(), j * rhs.cols(), rhs.rows(), rhs.cols()) = lhs(i, j) * rhs; @@ -85,9 +86,9 @@ inline ComplexMatrix kronecker_product(const ComplexMatrix& lhs, const ComplexMa return result; } -inline ComplexMatrix get_expanded_matrix(const ComplexMatrix& from_matrix, - const std::vector& from_targets, - std::vector& to_targets) { +inline internal::ComplexMatrix get_expanded_matrix(const internal::ComplexMatrix& from_matrix, + const std::vector& from_targets, + std::vector& to_targets) { std::vector targets_map(from_targets.size()); std::ranges::transform(from_targets, targets_map.begin(), [&](std::uint64_t x) { return std::ranges::lower_bound(to_targets, x) - to_targets.begin(); @@ -105,8 +106,8 @@ inline ComplexMatrix get_expanded_matrix(const ComplexMatrix& from_matrix, for (std::uint64_t i : std::views::iota(0ULL, 1ULL << to_targets.size())) { if ((i & targets_idx_mask) == 0) outer_indices.push_back(i); } - ComplexMatrix to_matrix = - ComplexMatrix::Zero(1ULL << to_targets.size(), 1ULL << to_targets.size()); + internal::ComplexMatrix to_matrix = + internal::ComplexMatrix::Zero(1ULL << to_targets.size(), 1ULL << to_targets.size()); for (std::uint64_t i : std::views::iota(0ULL, 1ULL << from_targets.size())) { for (std::uint64_t j : std::views::iota(0ULL, 1ULL << from_targets.size())) { for (std::uint64_t o : outer_indices) { diff --git a/tests/circuit/circuit_test.cpp b/tests/circuit/circuit_test.cpp index acbac053..9d139037 100644 --- a/tests/circuit/circuit_test.cpp +++ b/tests/circuit/circuit_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include "../util/util.hpp" #include "gate/gate_factory.hpp" @@ -22,7 +23,7 @@ TEST(CircuitTest, CircuitBasic) { StateVector state = StateVector::Haar_random_state(n); Eigen::VectorXcd state_eigen(dim); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (std::uint64_t i = 0; i < dim; ++i) state_eigen[i] = state_cp[i]; Circuit circuit(n); @@ -124,11 +125,11 @@ TEST(CircuitTest, CircuitBasic) { state_eigen = get_eigen_matrix_full_qubit_Swap(target, target_sub, n) * state_eigen; target = random.int32() % n; - circuit.add_gate(gate::U1(target, M_PI)); + circuit.add_gate(gate::U1(target, std::numbers::pi)); state_eigen = get_expanded_eigen_matrix_with_identity(target, make_Z(), n) * state_eigen; target = random.int32() % n; - circuit.add_gate(gate::U2(target, 0, M_PI)); + circuit.add_gate(gate::U2(target, 0, std::numbers::pi)); state_eigen = get_expanded_eigen_matrix_with_identity(target, make_H(), n) * state_eigen; target = random.int32() % n; @@ -146,7 +147,7 @@ TEST(CircuitTest, CircuitBasic) { PauliOperator pauli = PauliOperator("I 0 X 1 Y 2 Z 3"); circuit.add_gate(multi_Pauli(pauli)); - ComplexMatrix mat_x(2, 2); + internal::ComplexMatrix mat_x(2, 2); target = random.int32() % n; mat_x << 0, 1, 1, 0; circuit.add_gate(dense_matrix(target, mat_x)); @@ -158,7 +159,7 @@ TEST(CircuitTest, CircuitBasic) { circuit.update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); for (std::uint64_t i = 0; i < dim; ++i) ASSERT_NEAR(std::abs(state_eigen[i] - (CComplex)state_cp[i]), 0, eps); } @@ -170,7 +171,7 @@ TEST(CircuitTest, CircuitRev) { Random random; StateVector state = StateVector::Haar_random_state(n); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); Eigen::VectorXcd state_eigen(dim); for (std::uint64_t i = 0; i < dim; ++i) state_eigen[i] = state_cp[i]; @@ -243,7 +244,7 @@ TEST(CircuitTest, CircuitRev) { /* Observable observable(n); - angle = 2 * PI * random.uniform(); + angle = 2 * std::numbers::pi * random.uniform(); observable.add_operator(1.0, "Z 0"); observable.add_operator(2.0, "Z 0 Z 1"); @@ -256,7 +257,7 @@ TEST(CircuitTest, CircuitRev) { auto revcircuit = circuit.get_inverse(); revcircuit.update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); for (std::uint64_t i = 0; i < dim; ++i) ASSERT_NEAR(abs(state_eigen[i] - (CComplex)state_cp[i]), 0, eps); } diff --git a/tests/circuit/param_circuit_test.cpp b/tests/circuit/param_circuit_test.cpp index 1fd7f553..8f0fdb60 100644 --- a/tests/circuit/param_circuit_test.cpp +++ b/tests/circuit/param_circuit_test.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -25,7 +26,7 @@ TEST(ParamCircuitTest, ApplyParamCircuit) { pkeys[pidx] = "p" + std::to_string(pidx); std::array params = {}; for (std::uint64_t pidx : std::views::iota(std::uint64_t{0}, nparams)) - params[pidx] = random.uniform() * M_PI * 2; + params[pidx] = random.uniform() * std::numbers::pi * 2; std::map pmap; for (std::uint64_t pidx : std::views::iota(std::uint64_t{0}, nparams)) pmap[pkeys[pidx]] = params[pidx]; @@ -37,15 +38,15 @@ TEST(ParamCircuitTest, ApplyParamCircuit) { if (param_gate_kind == 0) { std::uint64_t target = random.int32() % n_qubits; circuit.add_gate(gate::RX(target, coef * params[pidx])); - pcircuit.add_param_gate(gate::PRX(target, coef), pkeys[pidx]); + pcircuit.add_param_gate(gate::ParamRX(target, coef), pkeys[pidx]); } else if (param_gate_kind == 1) { std::uint64_t target = random.int32() % n_qubits; circuit.add_gate(gate::RY(target, coef * params[pidx])); - pcircuit.add_param_gate(gate::PRY(target, coef), pkeys[pidx]); + pcircuit.add_param_gate(gate::ParamRY(target, coef), pkeys[pidx]); } else if (param_gate_kind == 2) { std::uint64_t target = random.int32() % n_qubits; circuit.add_gate(gate::RZ(target, coef * params[pidx])); - pcircuit.add_param_gate(gate::PRZ(target, coef), pkeys[pidx]); + pcircuit.add_param_gate(gate::ParamRZ(target, coef), pkeys[pidx]); } else { std::vector target_vec, pauli_id_vec; for (std::uint64_t target = 0; target < n_qubits; target++) { @@ -54,7 +55,7 @@ TEST(ParamCircuitTest, ApplyParamCircuit) { } PauliOperator pauli(target_vec, pauli_id_vec, 1.0); circuit.add_gate(gate::PauliRotation(pauli, coef * params[pidx])); - pcircuit.add_param_gate(gate::PPauliRotation(pauli, coef), pkeys[pidx]); + pcircuit.add_param_gate(gate::ParamPauliRotation(pauli, coef), pkeys[pidx]); } } else { std::uint64_t control = random.int32() % n_qubits; @@ -69,11 +70,11 @@ TEST(ParamCircuitTest, ApplyParamCircuit) { } StateVector state = StateVector::Haar_random_state(n_qubits); StateVector state_cp = state.copy(); - auto amplitudes_base = state.amplitudes(); + auto amplitudes_base = state.get_amplitudes(); circuit.update_quantum_state(state); pcircuit.update_quantum_state(state_cp, pmap); - auto amplitudes = state.amplitudes(); - auto amplitudes_cp = state_cp.amplitudes(); + auto amplitudes = state.get_amplitudes(); + auto amplitudes_cp = state_cp.get_amplitudes(); for (std::uint64_t idx : std::views::iota(std::uint64_t{0}, 1ULL << n_qubits)) { ASSERT_NEAR(Kokkos::abs(amplitudes[idx] - amplitudes_cp[idx]), 0, eps); } @@ -81,8 +82,8 @@ TEST(ParamCircuitTest, ApplyParamCircuit) { auto ipcircuit = pcircuit.get_inverse(); icircuit.update_quantum_state(state); ipcircuit.update_quantum_state(state_cp, pmap); - amplitudes = state.amplitudes(); - amplitudes_cp = state_cp.amplitudes(); + amplitudes = state.get_amplitudes(); + amplitudes_cp = state_cp.get_amplitudes(); for (std::uint64_t idx : std::views::iota(std::uint64_t{0}, 1ULL << n_qubits)) { ASSERT_NEAR(Kokkos::abs(amplitudes[idx] - amplitudes_base[idx]), 0, eps); ASSERT_NEAR(Kokkos::abs(amplitudes_cp[idx] - amplitudes_base[idx]), 0, eps); @@ -92,9 +93,9 @@ TEST(ParamCircuitTest, ApplyParamCircuit) { TEST(ParamCircuitTest, InsufficientParameterGiven) { Circuit circuit(1); - circuit.add_param_gate(gate::PRX(0), "0"); - circuit.add_param_gate(gate::PRX(0), "1"); - circuit.add_param_gate(gate::PRX(0), "0"); + circuit.add_param_gate(gate::ParamRX(0), "0"); + circuit.add_param_gate(gate::ParamRX(0), "1"); + circuit.add_param_gate(gate::ParamRX(0), "0"); StateVector state(1); ASSERT_NO_THROW(circuit.update_quantum_state(state, {{"0", 0}, {"1", 0}})); ASSERT_NO_THROW(circuit.update_quantum_state(state, {{"0", 0}, {"1", 0}, {"2", 0}})); diff --git a/tests/gate/gate_test.cpp b/tests/gate/gate_test.cpp index d78ad185..87ca3f92 100644 --- a/tests/gate/gate_test.cpp +++ b/tests/gate/gate_test.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -23,14 +24,14 @@ void run_random_gate_apply(std::uint64_t n_qubits) { Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim); for (int repeat = 0; repeat < 10; repeat++) { auto state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } const Gate gate = QuantumGateConstructor(); gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = test_state; @@ -48,15 +49,15 @@ void run_random_gate_apply(std::uint64_t n_qubits) { Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim); for (int repeat = 0; repeat < 10; repeat++) { auto state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } - const double angle = M_PI * random.uniform(); + const double angle = std::numbers::pi * random.uniform(); const Gate gate = QuantumGateConstructor(angle, {}); gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = std::polar(1., angle) * test_state; @@ -76,7 +77,7 @@ void run_random_gate_apply(std::uint64_t n_qubits, Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim); for (int repeat = 0; repeat < 10; repeat++) { auto state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } @@ -84,7 +85,7 @@ void run_random_gate_apply(std::uint64_t n_qubits, const std::uint64_t target = random.int64() % n_qubits; const Gate gate = QuantumGateConstructor(target, {}); gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state; @@ -103,17 +104,17 @@ void run_random_gate_apply(std::uint64_t n_qubits, Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim); for (int repeat = 0; repeat < 10; repeat++) { auto state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } - const double angle = M_PI * random.uniform(); + const double angle = std::numbers::pi * random.uniform(); const auto matrix = matrix_factory(angle); const std::uint64_t target = random.int64() % n_qubits; const Gate gate = QuantumGateConstructor(target, angle, {}); gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state; @@ -132,20 +133,20 @@ void run_random_gate_apply_IBMQ( Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim); for (int repeat = 0; repeat < 10; repeat++) { auto state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int gate_type = 0; gate_type < 3; gate_type++) { for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } - double theta = M_PI * random.uniform(); - double phi = M_PI * random.uniform(); - double lambda = M_PI * random.uniform(); + double theta = std::numbers::pi * random.uniform(); + double phi = std::numbers::pi * random.uniform(); + double lambda = std::numbers::pi * random.uniform(); if (gate_type == 0) { theta = 0; phi = 0; } else if (gate_type == 1) { - theta = M_PI / 2; + theta = std::numbers::pi / 2; } const auto matrix = matrix_factory(theta, phi, lambda); const std::uint64_t target = random.int64() % n_qubits; @@ -158,7 +159,7 @@ void run_random_gate_apply_IBMQ( gate = gate::U3(target, theta, phi, lambda, {}); } gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state; @@ -180,7 +181,7 @@ void run_random_gate_apply_two_target(std::uint64_t n_qubits) { auto state = StateVector::Haar_random_state(n_qubits); for (int g = 0; g < 2; g++) { Gate gate; - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } @@ -196,7 +197,7 @@ void run_random_gate_apply_two_target(std::uint64_t n_qubits) { func_eig = get_eigen_matrix_full_qubit_CZ; } gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); Eigen::MatrixXcd test_mat = func_eig(control, target, n_qubits); test_state = test_mat * test_state; @@ -209,7 +210,7 @@ void run_random_gate_apply_two_target(std::uint64_t n_qubits) { for (int repeat = 0; repeat < 10; repeat++) { auto state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (int i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } @@ -219,7 +220,7 @@ void run_random_gate_apply_two_target(std::uint64_t n_qubits) { if (target1 == target2) target1 = (target1 + 1) % n_qubits; auto gate = gate::Swap(target1, target2); gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); Eigen::MatrixXcd test_mat = get_eigen_matrix_full_qubit_Swap(target1, target2, n_qubits); test_state = test_mat * test_state; @@ -239,7 +240,7 @@ void run_random_gate_apply_pauli(std::uint64_t n_qubits) { // Test for PauliGate for (int repeat = 0; repeat < 10; repeat++) { StateVector state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); auto state_bef = state.copy(); for (std::uint64_t i = 0; i < dim; i++) { @@ -277,7 +278,7 @@ void run_random_gate_apply_pauli(std::uint64_t n_qubits) { Gate pauli_gate = gate::Pauli(pauli); pauli_gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = matrix * test_state; // check if the state is updated correctly @@ -285,10 +286,10 @@ void run_random_gate_apply_pauli(std::uint64_t n_qubits) { ASSERT_NEAR(std::abs((CComplex)state_cp[i] - test_state[i]), 0, eps); } - auto state_bef_cp = state_bef.amplitudes(); + auto state_bef_cp = state_bef.get_amplitudes(); Gate pauli_gate_inv = pauli_gate->get_inverse(); pauli_gate_inv->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); // check if the state is restored correctly for (std::uint64_t i = 0; i < dim; i++) { @@ -299,13 +300,13 @@ void run_random_gate_apply_pauli(std::uint64_t n_qubits) { // Test for PauliRotationGate for (int repeat = 0; repeat < 10; repeat++) { StateVector state = StateVector::Haar_random_state(n_qubits); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); auto state_bef = state.copy(); assert(test_state.size() == (int)state_cp.size()); for (std::uint64_t i = 0; i < dim; i++) { test_state[i] = state_cp[i]; } - const double angle = M_PI * random.uniform(); + const double angle = std::numbers::pi * random.uniform(); std::vector target_vec, pauli_id_vec; for (std::uint64_t target = 0; target < n_qubits; target++) { target_vec.emplace_back(target); @@ -333,11 +334,11 @@ void run_random_gate_apply_pauli(std::uint64_t n_qubits) { } } matrix = std::cos(angle / 2) * Eigen::MatrixXcd::Identity(dim, dim) - - Complex(0, 1) * std::sin(angle / 2) * matrix; + StdComplex(0, 1) * std::sin(angle / 2) * matrix; PauliOperator pauli(target_vec, pauli_id_vec, 1.0); Gate pauli_gate = gate::PauliRotation(pauli, angle); pauli_gate->update_quantum_state(state); - state_cp = state.amplitudes(); + state_cp = state.get_amplitudes(); test_state = matrix * test_state; assert((int)state_cp.size() == test_state.size()); // check if the state is updated correctly @@ -346,8 +347,8 @@ void run_random_gate_apply_pauli(std::uint64_t n_qubits) { } Gate pauli_gate_inv = pauli_gate->get_inverse(); pauli_gate_inv->update_quantum_state(state); - state_cp = state.amplitudes(); - auto state_bef_cp = state_bef.amplitudes(); + state_cp = state.get_amplitudes(); + auto state_bef_cp = state_bef.get_amplitudes(); // check if the state is restored correctly for (std::uint64_t i = 0; i < dim; i++) { ASSERT_NEAR(std::abs((CComplex)(state_cp[i] - state_bef_cp[i])), 0, eps); @@ -613,3 +614,151 @@ TEST(GateTest, ApplyProbablisticGate) { ASSERT_GT(i_cnt, 0); ASSERT_LT(x_cnt, i_cnt); } + +void test_gate(Gate gate_control, + Gate gate_simple, + std::uint64_t n_qubits, + std::uint64_t control_mask) { + StateVector state = StateVector::Haar_random_state(n_qubits); + auto amplitudes = state.get_amplitudes(); + StateVector state_controlled(n_qubits - std::popcount(control_mask)); + std::vector amplitudes_controlled(state_controlled.dim()); + for (std::uint64_t i : std::views::iota(0ULL, state_controlled.dim())) { + amplitudes_controlled[i] = + amplitudes[internal::insert_zero_at_mask_positions(i, control_mask) | control_mask]; + } + state_controlled.load(amplitudes_controlled); + gate_control->update_quantum_state(state); + gate_simple->update_quantum_state(state_controlled); + amplitudes = state.get_amplitudes(); + amplitudes_controlled = state_controlled.get_amplitudes(); + for (std::uint64_t i : std::views::iota(0ULL, state_controlled.dim())) { + ASSERT_NEAR( + Kokkos::abs(amplitudes_controlled[i] - + amplitudes[internal::insert_zero_at_mask_positions(i, control_mask) | + control_mask]), + 0., + eps); + } +} + +template +void test_standard_gate_control(Factory factory, std::uint64_t n) { + Random random; + std::vector shuffled(n); + std::iota(shuffled.begin(), shuffled.end(), 0ULL); + for (std::uint64_t i : std::views::iota(0ULL, n) | std::views::reverse) { + std::uint64_t j = random.int32() % (i + 1); + if (i != j) std::swap(shuffled[i], shuffled[j]); + } + std::vector targets(num_target); + for (std::uint64_t i : std::views::iota(0ULL, num_target)) { + targets[i] = shuffled[i]; + } + std::uint64_t num_control = random.int32() % (n - num_target + 1); + std::vector controls(num_control); + for (std::uint64_t i : std::views::iota(0ULL, num_control)) { + controls[i] = shuffled[num_target + i]; + } + std::uint64_t control_mask = 0ULL; + for (std::uint64_t c : controls) control_mask |= 1ULL << c; + std::vector angles(num_rotation); + for (double& angle : angles) angle = random.uniform() * std::numbers::pi * 2; + if constexpr (num_target == 0 && num_rotation == 1) { + Gate g1 = factory(angles[0], controls); + Gate g2 = factory(angles[0], {}); + test_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 0) { + Gate g1 = factory(targets[0], controls); + Gate g2 = + factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), {}); + test_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 1) { + Gate g1 = factory(targets[0], angles[0], controls); + Gate g2 = factory( + targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), angles[0], {}); + test_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 2) { + Gate g1 = factory(targets[0], angles[0], angles[1], controls); + Gate g2 = factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), + angles[0], + angles[1], + {}); + test_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 3) { + Gate g1 = factory(targets[0], angles[0], angles[1], angles[2], controls); + Gate g2 = factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), + angles[0], + angles[1], + angles[2], + {}); + test_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 2 && num_rotation == 0) { + Gate g1 = factory(targets[0], targets[1], controls); + Gate g2 = factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), + targets[1] - std::popcount(control_mask & ((1ULL << targets[1]) - 1)), + {}); + test_gate(g1, g2, n, control_mask); + } else { + FAIL(); + } +} + +template +void test_pauli_control(std::uint64_t n) { + PauliOperator::Data data1, data2; + std::vector controls; + std::uint64_t control_mask = 0; + std::uint64_t num_control = 0; + Random random; + for (std::uint64_t i : std::views::iota(0ULL, n)) { + std::uint64_t dat = random.int32() % 12; + if (dat < 4) { + data1.add_single_pauli(i, dat); + data2.add_single_pauli(i - num_control, dat); + } else if (dat < 8) { + controls.push_back(i); + control_mask |= 1ULL << i; + num_control++; + } + } + if constexpr (!rotation) { + Gate g1 = gate::Pauli(PauliOperator(data1), controls); + Gate g2 = gate::Pauli(PauliOperator(data2), {}); + test_gate(g1, g2, n, control_mask); + } else { + double angle = random.uniform() * std::numbers::pi * 2; + Gate g1 = gate::PauliRotation(PauliOperator(data1), angle, controls); + Gate g2 = gate::PauliRotation(PauliOperator(data2), angle, {}); + test_gate(g1, g2, n, control_mask); + } +} + +TEST(GateTest, Control) { + std::uint64_t n = 10; + for ([[maybe_unused]] std::uint64_t _ : std::views::iota(0, 10)) { + test_standard_gate_control<0, 1>(gate::GlobalPhase, n); + test_standard_gate_control<1, 0>(gate::X, n); + test_standard_gate_control<1, 0>(gate::Y, n); + test_standard_gate_control<1, 0>(gate::Z, n); + test_standard_gate_control<1, 0>(gate::S, n); + test_standard_gate_control<1, 0>(gate::Sdag, n); + test_standard_gate_control<1, 0>(gate::T, n); + test_standard_gate_control<1, 0>(gate::Tdag, n); + test_standard_gate_control<1, 0>(gate::SqrtX, n); + test_standard_gate_control<1, 0>(gate::SqrtXdag, n); + test_standard_gate_control<1, 0>(gate::SqrtY, n); + test_standard_gate_control<1, 0>(gate::SqrtYdag, n); + test_standard_gate_control<1, 0>(gate::P0, n); + test_standard_gate_control<1, 0>(gate::P1, n); + test_standard_gate_control<1, 1>(gate::RX, n); + test_standard_gate_control<1, 1>(gate::RY, n); + test_standard_gate_control<1, 1>(gate::RZ, n); + test_standard_gate_control<1, 1>(gate::U1, n); + test_standard_gate_control<1, 2>(gate::U2, n); + test_standard_gate_control<1, 3>(gate::U3, n); + test_standard_gate_control<2, 0>(gate::Swap, n); + test_pauli_control(n); + test_pauli_control(n); + } +} diff --git a/tests/gate/merge_test.cpp b/tests/gate/merge_test.cpp index 590f36de..98234d16 100644 --- a/tests/gate/merge_test.cpp +++ b/tests/gate/merge_test.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -29,15 +30,17 @@ TEST(GateTest, MergeGate) { gates.push_back(gate::SqrtYdag(target)); gates.push_back(gate::P0(target)); gates.push_back(gate::P1(target)); - gates.push_back(gate::RX(target, random.uniform() * PI() * 2)); - gates.push_back(gate::RY(target, random.uniform() * PI() * 2)); - gates.push_back(gate::RZ(target, random.uniform() * PI() * 2)); - gates.push_back(gate::U1(target, random.uniform() * PI() * 2)); - gates.push_back(gate::U2(target, random.uniform() * PI() * 2, random.uniform() * PI() * 2)); + gates.push_back(gate::RX(target, random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::RY(target, random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::RZ(target, random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::U1(target, random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::U2(target, + random.uniform() * std::numbers::pi * 2, + random.uniform() * std::numbers::pi * 2)); gates.push_back(gate::U3(target, - random.uniform() * PI() * 2, - random.uniform() * PI() * 2, - random.uniform() * PI() * 2)); + random.uniform() * std::numbers::pi * 2, + random.uniform() * std::numbers::pi * 2, + random.uniform() * std::numbers::pi * 2)); gates.push_back( gate::OneTargetMatrix(target, {std::array{Complex(random.uniform(), random.uniform()), @@ -49,7 +52,7 @@ TEST(GateTest, MergeGate) { gates.push_back(gate::Swap(target, target ^ 1)); } gates.push_back(gate::I()); - gates.push_back(gate::GlobalPhase(random.uniform() * PI() * 2)); + gates.push_back(gate::GlobalPhase(random.uniform() * std::numbers::pi * 2)); gates.push_back( gate::TwoTargetMatrix(0, 1, @@ -73,11 +76,11 @@ TEST(GateTest, MergeGate) { gates.push_back(gate::Pauli(PauliOperator("Z 0", random.uniform()))); gates.push_back(gate::Pauli(PauliOperator("Z 1", random.uniform()))); gates.push_back(gate::PauliRotation(PauliOperator("X 0 Y 1", random.uniform()), - random.uniform() * PI() * 2)); - gates.push_back( - gate::PauliRotation(PauliOperator("Z 0", random.uniform()), random.uniform() * PI() * 2)); - gates.push_back( - gate::PauliRotation(PauliOperator("Z 1", random.uniform()), random.uniform() * PI() * 2)); + random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::PauliRotation(PauliOperator("Z 0", random.uniform()), + random.uniform() * std::numbers::pi * 2)); + gates.push_back(gate::PauliRotation(PauliOperator("Z 1", random.uniform()), + random.uniform() * std::numbers::pi * 2)); for (auto&& g1 : gates) { for (auto&& g2 : gates) { std::uint64_t n = 2; diff --git a/tests/gate/param_gate_test.cpp b/tests/gate/param_gate_test.cpp index 526a50b7..2d5e5c6f 100644 --- a/tests/gate/param_gate_test.cpp +++ b/tests/gate/param_gate_test.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -25,14 +26,14 @@ void test_apply_parametric_single_pauli_rotation(std::uint64_t n_qubits, auto state_bef = state.copy(); const std::uint64_t target = random.int32() % n_qubits; - const double param = M_PI * random.uniform(); - const double pcoef = random.uniform() * 2 - 1; - const Gate gate = factory_fixed(target, pcoef * param, {}); - const ParamGate pgate = factory_parametric(target, pcoef, {}); + const double param = std::numbers::pi * random.uniform(); + const double param_coef = random.uniform() * 2 - 1; + const Gate gate = factory_fixed(target, param_coef * param, {}); + const ParamGate pgate = factory_parametric(target, param_coef, {}); gate->update_quantum_state(state); pgate->update_quantum_state(state_cp, param); - auto state_amp = state.amplitudes(); - auto state_cp_amp = state_cp.amplitudes(); + auto state_amp = state.get_amplitudes(); + auto state_cp_amp = state_cp.get_amplitudes(); for (std::uint64_t i = 0; i < dim; i++) { ASSERT_NEAR(Kokkos::abs(state_cp_amp[i] - state_amp[i]), 0, eps); @@ -40,8 +41,8 @@ void test_apply_parametric_single_pauli_rotation(std::uint64_t n_qubits, ParamGate pgate_inv = pgate->get_inverse(); pgate_inv->update_quantum_state(state, param); - state_amp = state.amplitudes(); - auto state_bef_amp = state_bef.amplitudes(); + state_amp = state.get_amplitudes(); + auto state_bef_amp = state_bef.get_amplitudes(); for (std::uint64_t i = 0; i < dim; i++) { ASSERT_NEAR(Kokkos::abs(state_bef_amp[i] - state_amp[i]), 0, eps); } @@ -56,8 +57,8 @@ void test_apply_parametric_multi_pauli_rotation(std::uint64_t n_qubits) { StateVector state = StateVector::Haar_random_state(n_qubits); auto state_cp = state.copy(); auto state_bef = state.copy(); - const double param = M_PI * random.uniform(); - const double pcoef = random.uniform() * 2 - 1; + const double param = std::numbers::pi * random.uniform(); + const double param_coef = random.uniform() * 2 - 1; std::vector target_vec, pauli_id_vec; for (std::uint64_t target = 0; target < n_qubits; target++) { target_vec.emplace_back(target); @@ -65,20 +66,20 @@ void test_apply_parametric_multi_pauli_rotation(std::uint64_t n_qubits) { } PauliOperator pauli(target_vec, pauli_id_vec, 1.0); - Gate gate = gate::PauliRotation(pauli, pcoef * param); - ParamGate pgate = gate::PPauliRotation(pauli, pcoef); + Gate gate = gate::PauliRotation(pauli, param_coef * param); + ParamGate pgate = gate::ParamPauliRotation(pauli, param_coef); gate->update_quantum_state(state); pgate->update_quantum_state(state_cp, param); - auto state_amp = state.amplitudes(); - auto state_cp_amp = state_cp.amplitudes(); + auto state_amp = state.get_amplitudes(); + auto state_cp_amp = state_cp.get_amplitudes(); // check if the state is updated correctly for (std::uint64_t i = 0; i < dim; i++) { ASSERT_NEAR(Kokkos::abs(state_cp_amp[i] - state_amp[i]), 0, eps); } ParamGate pgate_inv = pgate->get_inverse(); pgate_inv->update_quantum_state(state, param); - state_amp = state.amplitudes(); - auto state_bef_amp = state_bef.amplitudes(); + state_amp = state.get_amplitudes(); + auto state_bef_amp = state_bef.get_amplitudes(); // check if the state is restored correctly for (std::uint64_t i = 0; i < dim; i++) { ASSERT_NEAR(Kokkos::abs((state_bef_amp[i] - state_amp[i])), 0, eps); @@ -86,24 +87,24 @@ void test_apply_parametric_multi_pauli_rotation(std::uint64_t n_qubits) { } } -TEST(ParamGateTest, ApplyPRXGate) { - test_apply_parametric_single_pauli_rotation(5, &gate::RX, &gate::PRX); +TEST(ParamGateTest, ApplyParamRXGate) { + test_apply_parametric_single_pauli_rotation(5, &gate::RX, &gate::ParamRX); } -TEST(ParamGateTest, ApplyPRYGate) { - test_apply_parametric_single_pauli_rotation(5, &gate::RX, &gate::PRX); +TEST(ParamGateTest, ApplyParamRYGate) { + test_apply_parametric_single_pauli_rotation(5, &gate::RX, &gate::ParamRX); } -TEST(ParamGateTest, ApplyPRZGate) { - test_apply_parametric_single_pauli_rotation(5, &gate::RX, &gate::PRX); +TEST(ParamGateTest, ApplyParamRZGate) { + test_apply_parametric_single_pauli_rotation(5, &gate::RX, &gate::ParamRX); } -TEST(ParamGateTest, ApplyPPauliRotationGate) { test_apply_parametric_multi_pauli_rotation(5); } +TEST(ParamGateTest, ApplyParamPauliRotationGate) { test_apply_parametric_multi_pauli_rotation(5); } -TEST(ParamGateTest, ApplyPProbablisticGate) { - auto probgate = gate::PProbablistic({.1, .9}, {gate::PRX(0), gate::I()}); +TEST(ParamGateTest, ApplyParamProbablisticGate) { + auto probgate = gate::ParamProbablistic({.1, .9}, {gate::ParamRX(0), gate::I()}); std::uint64_t x_cnt = 0, i_cnt = 0; StateVector state(1); for ([[maybe_unused]] auto _ : std::views::iota(0, 100)) { std::uint64_t before = state.sampling(1)[0]; - probgate->update_quantum_state(state, scaluq::PI()); + probgate->update_quantum_state(state, std::numbers::pi); std::uint64_t after = state.sampling(1)[0]; if (before != after) { x_cnt++; @@ -116,3 +117,89 @@ TEST(ParamGateTest, ApplyPProbablisticGate) { ASSERT_GT(i_cnt, 0); ASSERT_LT(x_cnt, i_cnt); } + +void test_gate(ParamGate gate_control, + ParamGate gate_simple, + std::uint64_t n_qubits, + std::uint64_t control_mask, + double param) { + StateVector state = StateVector::Haar_random_state(n_qubits); + auto amplitudes = state.get_amplitudes(); + StateVector state_controlled(n_qubits - std::popcount(control_mask)); + std::vector amplitudes_controlled(state_controlled.dim()); + for (std::uint64_t i : std::views::iota(0ULL, state_controlled.dim())) { + amplitudes_controlled[i] = + amplitudes[internal::insert_zero_at_mask_positions(i, control_mask) | control_mask]; + } + state_controlled.load(amplitudes_controlled); + gate_control->update_quantum_state(state, param); + gate_simple->update_quantum_state(state_controlled, param); + amplitudes = state.get_amplitudes(); + amplitudes_controlled = state_controlled.get_amplitudes(); + for (std::uint64_t i : std::views::iota(0ULL, state_controlled.dim())) { + ASSERT_NEAR( + Kokkos::abs(amplitudes_controlled[i] - + amplitudes[internal::insert_zero_at_mask_positions(i, control_mask) | + control_mask]), + 0., + eps); + } +} + +template +void test_param_rotation_control(Factory factory, std::uint64_t n) { + std::cerr << "prx" << std::endl; + Random random; + std::vector shuffled(n); + std::iota(shuffled.begin(), shuffled.end(), 0ULL); + for (std::uint64_t i : std::views::iota(0ULL, n) | std::views::reverse) { + std::uint64_t j = random.int32() % (i + 1); + if (i != j) std::swap(shuffled[i], shuffled[j]); + } + std::uint64_t target = shuffled[0]; + std::uint64_t num_control = random.int32() % n; + std::vector controls(num_control); + for (std::uint64_t i : std::views::iota(0ULL, num_control)) { + controls[i] = shuffled[1 + i]; + } + std::uint64_t control_mask = 0ULL; + for (std::uint64_t c : controls) control_mask |= 1ULL << c; + double param = random.uniform() * std::numbers::pi * 2; + ParamGate g1 = factory(target, 1., controls); + ParamGate g2 = factory(target - std::popcount(control_mask & ((1ULL << target) - 1)), 1., {}); + test_gate(g1, g2, n, control_mask, param); +} + +void test_ppauli_control(std::uint64_t n) { + std::cerr << "ppauli" << std::endl; + PauliOperator::Data data1, data2; + std::vector controls; + std::uint64_t control_mask = 0; + std::uint64_t num_control = 0; + Random random; + for (std::uint64_t i : std::views::iota(0ULL, n)) { + std::uint64_t dat = random.int32() % 12; + if (dat < 4) { + data1.add_single_pauli(i, dat); + data2.add_single_pauli(i - num_control, dat); + } else if (dat < 8) { + controls.push_back(i); + control_mask |= 1ULL << i; + num_control++; + } + } + double param = random.uniform() * std::numbers::pi * 2; + ParamGate g1 = gate::ParamPauliRotation(PauliOperator(data1), 1., controls); + ParamGate g2 = gate::ParamPauliRotation(PauliOperator(data2), 1., {}); + test_gate(g1, g2, n, control_mask, param); +} + +TEST(ParamGateTest, Control) { + std::uint64_t n = 10; + for ([[maybe_unused]] std::uint64_t _ : std::views::iota(0, 10)) { + test_param_rotation_control(gate::ParamRX, n); + test_param_rotation_control(gate::ParamRY, n); + test_param_rotation_control(gate::ParamRZ, n); + test_ppauli_control(n); + } +} diff --git a/tests/operator/test_operator.cpp b/tests/operator/test_operator.cpp index bf0b3a00..1b2dfeb0 100644 --- a/tests/operator/test_operator.cpp +++ b/tests/operator/test_operator.cpp @@ -56,7 +56,7 @@ TEST(OperatorTest, CheckExpectationValue) { generate_random_observable_with_eigen(n, random); auto state = StateVector::Haar_random_state(n); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim); for (std::uint64_t i = 0; i < dim; ++i) test_state[i] = state_cp[i]; @@ -78,11 +78,11 @@ TEST(OperatorTest, CheckTransitionAmplitude) { generate_random_observable_with_eigen(n, random); auto state_bra = StateVector::Haar_random_state(n); - auto state_bra_cp = state_bra.amplitudes(); + auto state_bra_cp = state_bra.get_amplitudes(); Eigen::VectorXcd test_state_bra = Eigen::VectorXcd::Zero(dim); for (std::uint64_t i = 0; i < dim; ++i) test_state_bra[i] = state_bra_cp[i]; auto state_ket = StateVector::Haar_random_state(n); - auto state_ket_cp = state_ket.amplitudes(); + auto state_ket_cp = state_ket.get_amplitudes(); Eigen::VectorXcd test_state_ket = Eigen::VectorXcd::Zero(dim); for (std::uint64_t i = 0; i < dim; ++i) test_state_ket[i] = state_ket_cp[i]; @@ -150,7 +150,7 @@ TEST(OperatorTest, ApplyToStateTest) { }()); Operator op(n_qubits); - op.add_operator({std::vector{1, 0, 0}, std::vector{0, 1, 0}, Complex(2)}); + op.add_operator({0b001, 0b010, Complex(2)}); op.add_operator({"X 2 Y 1", 1}); op.apply_to_state(state_vector); @@ -164,7 +164,7 @@ TEST(OperatorTest, ApplyToStateTest) { Complex(-14, 0), Complex(-12, 1), }; - ASSERT_EQ(state_vector.amplitudes(), expected); + ASSERT_EQ(state_vector.get_amplitudes(), expected); } TEST(OperatorTest, Optimize) { @@ -180,7 +180,7 @@ TEST(OperatorTest, Optimize) { {"X 0 Y 1", 10.}, {"Y 0 Z 1", 2.}, {"Z 1", 7.}}; std::vector> test; for (const auto& pauli : op.terms()) { - test.emplace_back(pauli.get_pauli_string(), pauli.get_coef()); + test.emplace_back(pauli.get_pauli_string(), pauli.coef()); } std::ranges::sort(expected, [](const auto& l, const auto& r) { return l.first < r.first; }); std::ranges::sort(test, [](const auto& l, const auto& r) { return l.first < r.first; }); diff --git a/tests/operator/test_pauli_operator.cpp b/tests/operator/test_pauli_operator.cpp index 3c1357a6..ead4a05c 100644 --- a/tests/operator/test_pauli_operator.cpp +++ b/tests/operator/test_pauli_operator.cpp @@ -13,15 +13,15 @@ TEST(PauliOperatorTest, ContainsExtraWhitespace) { PauliOperator expected = PauliOperator("X 0", 1.0); PauliOperator pauli_whitespace = PauliOperator("X 0 ", 1.0); - EXPECT_EQ(1, pauli_whitespace.get_target_qubit_list().size()); - EXPECT_EQ(1, pauli_whitespace.get_pauli_id_list().size()); + EXPECT_EQ(1, pauli_whitespace.target_qubit_list().size()); + EXPECT_EQ(1, pauli_whitespace.pauli_id_list().size()); EXPECT_EQ(expected.get_pauli_string(), pauli_whitespace.get_pauli_string()); } TEST(PauliOperatorTest, EmptyStringConstructsIdentity) { const auto identity = PauliOperator("", 1.0); - ASSERT_EQ(0, identity.get_target_qubit_list().size()); - ASSERT_EQ(0, identity.get_pauli_id_list().size()); + ASSERT_EQ(0, identity.target_qubit_list().size()); + ASSERT_EQ(0, identity.pauli_id_list().size()); ASSERT_EQ("", identity.get_pauli_string()); } @@ -62,7 +62,7 @@ TEST(PauliOperatorTest, SpacedPauliString) { double coef = 2.0; std::string Pauli_string = "X 0 Y 1 "; PauliOperator pauli = PauliOperator(Pauli_string, coef); - size_t PauliSize = pauli.get_target_qubit_list().size(); + size_t PauliSize = pauli.target_qubit_list().size(); ASSERT_EQ(PauliSize, 2); } @@ -95,7 +95,7 @@ TEST_P(PauliOperatorMultiplyTest, MultiplyTest) { const auto p = GetParam(); PauliOperator res = p.op1 * p.op2; EXPECT_EQ(p.expected.get_pauli_string(), res.get_pauli_string()); - EXPECT_EQ(p.expected.get_coef(), res.get_coef()); + EXPECT_EQ(p.expected.coef(), res.coef()); } TEST_P(PauliOperatorMultiplyTest, MultiplyAssignmentTest) { @@ -103,7 +103,7 @@ TEST_P(PauliOperatorMultiplyTest, MultiplyAssignmentTest) { PauliOperator res = p.op1; res = res * p.op2; EXPECT_EQ(p.expected.get_pauli_string(), res.get_pauli_string()); - EXPECT_EQ(p.expected.get_coef(), res.get_coef()); + EXPECT_EQ(p.expected.coef(), res.coef()); } INSTANTIATE_TEST_CASE_P( @@ -154,7 +154,7 @@ INSTANTIATE_TEST_CASE_P( PauliOperatorMultiplyTest, []() { double coef = 2.0; - std::uint64_t MAX_TERM = 100; + std::uint64_t MAX_TERM = 64; std::string pauli_string_x = ""; std::string pauli_string_y = ""; std::string pauli_string_z = ""; @@ -187,8 +187,8 @@ TEST(PauliOperatorTest, ApplyToStateTest) { return tmp; }()); - PauliOperator op(std::vector{1, 0, 0}, std::vector{0, 1, 0}, Complex(2)); + PauliOperator op(0b001, 0b010, Complex(2)); op.apply_to_state(state_vector); std::vector expected = {2, 0, -6, -4, 10, 8, -14, -12}; - ASSERT_EQ(state_vector.amplitudes(), expected); + ASSERT_EQ(state_vector.get_amplitudes(), expected); } diff --git a/tests/state/state_vector_batched_test.cpp b/tests/state/state_vector_batched_test.cpp index d156a8f0..354a555b 100644 --- a/tests/state/state_vector_batched_test.cpp +++ b/tests/state/state_vector_batched_test.cpp @@ -33,7 +33,7 @@ TEST(StateVectorBatchedTest, LoadAndAmplitues) { StateVectorBatched states(batch_size, n_qubits); states.load(states_h); - auto amps = states.amplitudes(); + auto amps = states.get_amplitudes(); for (std::uint64_t b = 0; b < batch_size; ++b) { for (std::uint64_t i = 0; i < dim; ++i) { ASSERT_EQ(amps[b][i].real(), b * dim + i); @@ -49,30 +49,30 @@ TEST(StateVectorBatchedTest, OperateState) { auto states_cp = states.copy(); for (std::uint64_t b = 0; b < batch_size; ++b) { - ASSERT_TRUE(same_state(states.get_state_vector(b), states_cp.get_state_vector(b))); + ASSERT_TRUE(same_state(states.get_state_vector_at(b), states_cp.get_state_vector_at(b))); } states.add_state_vector(states_add); for (std::uint64_t b = 0; b < batch_size; ++b) { - auto v = states_cp.get_state_vector(b); - v.add_state_vector(states_add.get_state_vector(b)); - ASSERT_TRUE(same_state(v, states.get_state_vector(b))); + auto v = states_cp.get_state_vector_at(b); + v.add_state_vector(states_add.get_state_vector_at(b)); + ASSERT_TRUE(same_state(v, states.get_state_vector_at(b))); } states_cp = states.copy(); states.add_state_vector_with_coef(coef, states_add); for (std::uint64_t b = 0; b < batch_size; ++b) { - auto v = states_cp.get_state_vector(b); - v.add_state_vector_with_coef(coef, states_add.get_state_vector(b)); - ASSERT_TRUE(same_state(v, states.get_state_vector(b))); + auto v = states_cp.get_state_vector_at(b); + v.add_state_vector_with_coef(coef, states_add.get_state_vector_at(b)); + ASSERT_TRUE(same_state(v, states.get_state_vector_at(b))); } states_cp = states.copy(); states.multiply_coef(coef); for (std::uint64_t b = 0; b < batch_size; ++b) { - auto v = states_cp.get_state_vector(b); + auto v = states_cp.get_state_vector_at(b); v.multiply_coef(coef); - ASSERT_TRUE(same_state(v, states.get_state_vector(b))); + ASSERT_TRUE(same_state(v, states.get_state_vector_at(b))); } } @@ -83,7 +83,7 @@ TEST(StateVectorBatchedTest, ZeroProbs) { for (std::uint64_t i = 0; i < n_qubits; ++i) { auto zero_probs = states.get_zero_probability(i); for (std::uint64_t b = 0; b < batch_size; ++b) { - auto state = states.get_state_vector(b); + auto state = states.get_state_vector_at(b); ASSERT_NEAR(zero_probs[b], state.get_zero_probability(i), eps); } } @@ -101,7 +101,7 @@ TEST(StateVectorBatchedTest, MarginalProbs) { } auto mg_probs = states.get_marginal_probability(targets); for (std::uint64_t b = 0; b < batch_size; ++b) { - auto state = states.get_state_vector(b); + auto state = states.get_state_vector_at(b); ASSERT_NEAR(mg_probs[b], state.get_marginal_probability(targets), eps); } } @@ -113,7 +113,7 @@ TEST(StateVectorBatchedTest, Entropy) { auto entropies = states.get_entropy(); for (std::uint64_t b = 0; b < batch_size; ++b) { - auto state = states.get_state_vector(b); + auto state = states.get_state_vector_at(b); ASSERT_NEAR(entropies[b], state.get_entropy(), eps); } } diff --git a/tests/state/state_vector_test.cpp b/tests/state/state_vector_test.cpp index 30c5e895..e40356d1 100644 --- a/tests/state/state_vector_test.cpp +++ b/tests/state/state_vector_test.cpp @@ -25,9 +25,9 @@ TEST(StateVectorTest, HaarRandomStateNorm) { TEST(StateVectorTest, OperationAtIndex) { auto state = StateVector::Haar_random_state(10); for (std::uint64_t i = 0; i < state.dim(); ++i) { - state.set_amplitude_at_index(i, 1); - ASSERT_NEAR(state.get_amplitude_at_index(i).real(), 1, eps); - ASSERT_NEAR(state.get_amplitude_at_index(i).imag(), 0., eps); + state.set_amplitude_at(i, 1); + ASSERT_NEAR(state.get_amplitude_at(i).real(), 1, eps); + ASSERT_NEAR(state.get_amplitude_at(i).imag(), 0., eps); } } @@ -35,8 +35,8 @@ TEST(StateVectorTest, CopyState) { const int n = 5; const auto state = StateVector::Haar_random_state(n); StateVector state_cp = state.copy(); - auto vec1 = state.amplitudes(); - auto vec2 = state_cp.amplitudes(); + auto vec1 = state.get_amplitudes(); + auto vec2 = state_cp.get_amplitudes(); ASSERT_EQ(vec1, vec2); } @@ -45,7 +45,7 @@ TEST(StateVectorTest, ZeroNormState) { StateVector state(StateVector::Haar_random_state(n)); state.set_zero_norm_state(); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (std::uint64_t i = 0; i < state.dim(); ++i) { ASSERT_EQ((CComplex)state_cp[i], CComplex(0, 0)); @@ -57,7 +57,7 @@ TEST(StateVectorTest, ComputationalBasisState) { StateVector state(StateVector::Haar_random_state(n)); state.set_computational_basis(31); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); for (std::uint64_t i = 0; i < state.dim(); ++i) { if (i == 31) { @@ -90,10 +90,10 @@ TEST(StateVectorTest, AddState) { const std::uint64_t n = 10; StateVector state1(StateVector::Haar_random_state(n)); StateVector state2(StateVector::Haar_random_state(n)); - auto vec1 = state1.amplitudes(); - auto vec2 = state2.amplitudes(); + auto vec1 = state1.get_amplitudes(); + auto vec2 = state2.get_amplitudes(); state1.add_state_vector(state2); - auto new_vec = state1.amplitudes(); + auto new_vec = state1.get_amplitudes(); for (std::uint64_t i = 0; i < state1.dim(); ++i) { CComplex res = new_vec[i], val = (CComplex)vec1[i] + (CComplex)vec2[i]; @@ -107,11 +107,11 @@ TEST(StateVectorTest, AddStateWithCoef) { const std::uint64_t n = 10; StateVector state1(StateVector::Haar_random_state(n)); StateVector state2(StateVector::Haar_random_state(n)); - auto vec1 = state1.amplitudes(); - auto vec2 = state2.amplitudes(); + auto vec1 = state1.get_amplitudes(); + auto vec2 = state2.get_amplitudes(); state1.add_state_vector_with_coef(coef, state2); - auto new_vec = state1.amplitudes(); + auto new_vec = state1.get_amplitudes(); for (std::uint64_t i = 0; i < state1.dim(); ++i) { CComplex res = new_vec[i], val = (CComplex)vec1[i] + coef * (CComplex)vec2[i]; @@ -125,9 +125,9 @@ TEST(StateVectorTest, MultiplyCoef) { const CComplex coef(0.5, 0.2); StateVector state(StateVector::Haar_random_state(n)); - auto vec = state.amplitudes(); + auto vec = state.get_amplitudes(); state.multiply_coef(coef); - auto new_vec = state.amplitudes(); + auto new_vec = state.get_amplitudes(); for (std::uint64_t i = 0; i < state.dim(); ++i) { CComplex res = new_vec[i], val = coef * (CComplex)vec[i]; @@ -160,7 +160,7 @@ TEST(StateVectorTest, EntropyCalculation) { StateVector state(n); for (std::uint64_t rep = 0; rep < max_repeat; ++rep) { state = StateVector::Haar_random_state(n); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); ASSERT_NEAR(state.get_squared_norm(), 1, eps); Eigen::VectorXcd test_state(dim); for (std::uint64_t i = 0; i < dim; ++i) test_state[i] = (CComplex)state_cp[i]; @@ -181,7 +181,7 @@ TEST(StateVectorTest, GetMarginalProbability) { const std::uint64_t n = 2; const std::uint64_t dim = 1 << n; StateVector state(StateVector::Haar_random_state(n)); - auto state_cp = state.amplitudes(); + auto state_cp = state.get_amplitudes(); std::vector probs; for (std::uint64_t i = 0; i < dim; ++i) { probs.push_back(internal::squared_norm(state_cp[i])); diff --git a/tests/util/util.hpp b/tests/util/util.hpp index 1ef0fbe7..ba9e3e03 100644 --- a/tests/util/util.hpp +++ b/tests/util/util.hpp @@ -10,8 +10,8 @@ using namespace std::complex_literals; using namespace scaluq; inline bool same_state(const StateVector& s1, const StateVector& s2, const double eps = 1e-12) { - auto s1_cp = s1.amplitudes(); - auto s2_cp = s2.amplitudes(); + auto s1_cp = s1.get_amplitudes(); + auto s2_cp = s2.get_amplitudes(); assert(s1.n_qubits() == s2.n_qubits()); for (std::uint64_t i = 0; i < s1.dim(); ++i) { if (std::abs((std::complex)s1_cp[i] - (std::complex)s2_cp[i]) > eps) @@ -23,8 +23,8 @@ inline bool same_state(const StateVector& s1, const StateVector& s2, const doubl inline bool same_state_except_global_phase(const StateVector& s1, const StateVector& s2, const double eps = 1e-12) { - auto s1_cp = s1.amplitudes(); - auto s2_cp = s2.amplitudes(); + auto s1_cp = s1.get_amplitudes(); + auto s2_cp = s2.get_amplitudes(); assert(s1.n_qubits() == s2.n_qubits()); std::uint64_t significant = 0; for (std::uint64_t i = 0; i < s1.dim(); ++i) {