merge main

qulacs/CMakeLists.txt

@@ -1,15 +1,19 @@
 cmake_minimum_required(VERSION 3.21)
 
 target_sources(qulacs PRIVATE
+    gate/gate_npair_qubit.cpp
     gate/gate_one_control_one_target.cpp
     gate/gate_one_qubit.cpp
     gate/gate_quantum_matrix.cpp
-    gate/update_ops_matrix_dense.cpp
+    gate/gate_two_qubit.cpp
+    gate/update_ops_dense_matrix.cpp
+    gate/update_ops_npair_qubit.cpp
+    gate/update_ops_one_control_one_target.cpp
     gate/update_ops_one_qubit.cpp
+    gate/update_ops_quantum_matrix.cpp
+    gate/update_ops_two_qubit.cpp
     operator/pauli_operator.cpp
     operator/operator.cpp
-    gate/update_ops_one_target_one_control.cpp
-    gate/update_ops_quantum_matrix.cpp
     state/state_vector.cpp
     util/random.cpp
 )

qulacs/all.hpp

@@ -0,0 +1,15 @@
+#include "gate/constant.hpp"
+#include "gate/gate.hpp"
+#include "gate/gate_factory.hpp"
+#include "gate/gate_npair_qubit.hpp"
+#include "gate/gate_one_control_one_target.hpp"
+#include "gate/gate_one_qubit.hpp"
+#include "gate/gate_quantum_matrix.hpp"
+#include "gate/gate_two_qubit.hpp"
+#include "gate/update_ops.hpp"
+#include "operator/constant.hpp"
+#include "operator/operator.hpp"
+#include "operator/pauli_operator.hpp"
+#include "state/state_vector.hpp"
+#include "types.hpp"
+#include "util/random.hpp"

qulacs/gate/constant.hpp

@@ -27,47 +27,72 @@ KOKKOS_INLINE_FUNCTION
 double SINPI8() { return 0.382683432365090; }
 
 //! identity matrix
-matrix_2_2 I_GATE = {1, 0, 0, 1};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 I_GATE() { return {1, 0, 0, 1}; }
 //! Pauli matrix X
-matrix_2_2 X_GATE = {0, 1, 1, 0};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 X_GATE() { return {0, 1, 1, 0}; }
 //! Pauli matrix Y
-matrix_2_2 Y_GATE = {0, Complex(0, -1), Complex(0, 1), 0};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 Y_GATE() { return {0, Complex(0, -1), Complex(0, 1), 0}; }
 //! Pauli matrix Z
-matrix_2_2 Z_GATE = {1, 0, 0, -1};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 Z_GATE() { return {1, 0, 0, -1}; }
 
 //! list of Pauli matrix I,X,Y,Z
 // std::array<matrix_2_2, 4> PAULI_MATRIX = {I_GATE, X_GATE, Y_GATE, Z_GATE};
 
 //! S-gate
-matrix_2_2 S_GATE_MATRIX = {1, 0, 0, Complex(0, 1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 S_GATE_MATRIX() { return {1, 0, 0, Complex(0, 1)}; }
 //! Sdag-gate
-matrix_2_2 S_DAG_GATE_MATRIX = {1, 0, 0, Complex(0, -1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 S_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(0, -1)}; }
 //! T-gate
-matrix_2_2 T_GATE_MATRIX = {
-    COSPI8() - Complex(0, SINPI8()), 0., 0., COSPI8() + Complex(0, SINPI8()) * SINPI8()};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 T_GATE_MATRIX() {
+    return {COSPI8() - Complex(0, SINPI8()), 0., 0., COSPI8() + Complex(0, SINPI8()) * SINPI8()};
+}
 //! Tdag-gate
-matrix_2_2 T_DAG_GATE_MATRIX = {
-    COSPI8() + Complex(0, SINPI8()), 0., 0., COSPI8() - Complex(0, SINPI8())};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 T_DAG_GATE_MATRIX() {
+    return {COSPI8() + Complex(0, SINPI8()), 0., 0., COSPI8() - Complex(0, SINPI8())};
+}
 //! Hadamard gate
-matrix_2_2 HADAMARD_MATRIX = {INVERSE_SQRT2(), INVERSE_SQRT2(), INVERSE_SQRT2(), -INVERSE_SQRT2()};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 HADAMARD_MATRIX() {
+    return {INVERSE_SQRT2(), INVERSE_SQRT2(), INVERSE_SQRT2(), -INVERSE_SQRT2()};
+}
 //! square root of X gate
-matrix_2_2 SQRT_X_GATE_MATRIX = {
-    Complex(0.5, 0.5), Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(0.5, 0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 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
-matrix_2_2 SQRT_Y_GATE_MATRIX = {
-    Complex(0.5, 0.5), Complex(-0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 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
-matrix_2_2 SQRT_X_DAG_GATE_MATRIX = {
-    Complex(0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5), Complex(0.5, -0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 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
-matrix_2_2 SQRT_Y_DAG_GATE_MATRIX = {
-    Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(-0.5, 0.5), Complex(0.5, -0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 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
-matrix_2_2 PROJ_0_MATRIX = {1, 0, 0, 0};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PROJ_0_MATRIX() { return {1, 0, 0, 0}; }
 //! Projection to 1
-matrix_2_2 PROJ_1_MATRIX = {0, 0, 0, 1};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PROJ_1_MATRIX() { return {0, 0, 0, 1}; }
 //! complex values for exp(j * i*pi/4 )
-matrix_2_2 PHASE_90ROT = {1., Complex(0, 1), -1, Complex(0, -1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; }
 //! complex values for exp(-j * i*pi/4 )
-matrix_2_2 PHASE_M90ROT = {1., Complex(0, -1), -1, Complex(0, 1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; }
 }  // namespace qulacs

qulacs/gate/gate.hpp

@@ -4,8 +4,61 @@
 #include "../types.hpp"
 
 namespace qulacs {
-class QuantumGate {
+namespace internal {
+// forward declarations
+class GateBase;
+
+template <typename T>
+concept GateImpl = std::derived_from<T, GateBase>;
+
+template <GateImpl T>
+class GatePtr;
+}  // namespace internal
+using Gate = internal::GatePtr<internal::GateBase>;
+
+namespace internal {
+class GateBase {
 public:
+    virtual ~GateBase() = default;
+
+    [[nodiscard]] virtual std::vector<UINT> get_target_qubit_list() const = 0;
+    [[nodiscard]] virtual std::vector<UINT> get_control_qubit_list() const = 0;
+
+    [[nodiscard]] virtual Gate copy() const = 0;
+    [[nodiscard]] virtual Gate get_inverse() const = 0;
+
     virtual void update_quantum_state(StateVector& state_vector) const = 0;
 };
+
+template <GateImpl T>
+class GatePtr {
+    friend class GateFactory;
+    template <GateImpl U>
+    friend class GatePtr;
+
+private:
+    std::shared_ptr<T> _gate_ptr;
+
+public:
+    GatePtr() = default;
+    GatePtr(const GatePtr& gate) = default;
+    template <GateImpl U>
+    GatePtr(const std::shared_ptr<U>& gate_ptr) {
+        if constexpr (std::is_same_v<T, U>) {
+            _gate_ptr = gate_ptr;
+        } else if constexpr (std::is_same_v<U, internal::GateBase>) {
+            _gate_ptr = std::static_pointer_cast<T>(gate_ptr);
+        } else {
+            if (!(_gate_ptr = std::dynamic_pointer_cast<T>(gate_ptr))) {
+                throw std::runtime_error("invalid gate cast");
+            }
+        }
+    }
+    template <GateImpl U>
+    GatePtr(const GatePtr<U>& gate) : GatePtr(gate._gate_ptr) {}
+
+    T* operator->() const { return _gate_ptr.get(); }
+};
+}  // namespace internal
+
 }  // namespace qulacs

qulacs/gate/gate_factory.hpp

@@ -0,0 +1,78 @@
+#pragma once
+
+#include "gate/gate_npair_qubit.hpp"
+#include "gate/gate_one_control_one_target.hpp"
+#include "gate/gate_one_qubit.hpp"
+#include "gate/gate_quantum_matrix.hpp"
+#include "gate/gate_two_qubit.hpp"
+
+namespace qulacs {
+namespace internal {
+class GateFactory {
+public:
+    template <GateImpl T, typename... Args>
+    static Gate create_gate(Args... args) {
+        return {std::make_shared<T>(args...)};
+    }
+};
+}  // namespace internal
+
+Gate I(UINT target) { return internal::GateFactory::create_gate<internal::IGateImpl>(target); }
+Gate X(UINT target) { return internal::GateFactory::create_gate<internal::XGateImpl>(target); }
+Gate Y(UINT target) { return internal::GateFactory::create_gate<internal::YGateImpl>(target); }
+Gate Z(UINT target) { return internal::GateFactory::create_gate<internal::ZGateImpl>(target); }
+Gate H(UINT target) { return internal::GateFactory::create_gate<internal::HGateImpl>(target); }
+Gate S(UINT target) { return internal::GateFactory::create_gate<internal::SGateImpl>(target); }
+Gate Sdag(UINT target) {
+    return internal::GateFactory::create_gate<internal::SdagGateImpl>(target);
+}
+Gate T(UINT target) { return internal::GateFactory::create_gate<internal::TGateImpl>(target); }
+Gate Tdag(UINT target) {
+    return internal::GateFactory::create_gate<internal::TdagGateImpl>(target);
+}
+Gate sqrtX(UINT target) {
+    return internal::GateFactory::create_gate<internal::sqrtXGateImpl>(target);
+}
+Gate sqrtXdag(UINT target) {
+    return internal::GateFactory::create_gate<internal::sqrtXdagGateImpl>(target);
+}
+Gate sqrtY(UINT target) {
+    return internal::GateFactory::create_gate<internal::sqrtYGateImpl>(target);
+}
+Gate sqrtYdag(UINT target) {
+    return internal::GateFactory::create_gate<internal::sqrtYdagGateImpl>(target);
+}
+Gate P0(UINT target) { return internal::GateFactory::create_gate<internal::P0GateImpl>(target); }
+Gate P1(UINT target) { return internal::GateFactory::create_gate<internal::P1GateImpl>(target); }
+Gate RX(UINT target, double angle) {
+    return internal::GateFactory::create_gate<internal::RXGateImpl>(target, angle);
+}
+Gate RY(UINT target, double angle) {
+    return internal::GateFactory::create_gate<internal::RYGateImpl>(target, angle);
+}
+Gate RZ(UINT target, double angle) {
+    return internal::GateFactory::create_gate<internal::RZGateImpl>(target, angle);
+}
+Gate U1(UINT target, double lambda) {
+    return internal::GateFactory::create_gate<internal::U1GateImpl>(target, lambda);
+}
+Gate U2(UINT target, double phi, double lambda) {
+    return internal::GateFactory::create_gate<internal::U2GateImpl>(target, phi, lambda);
+}
+Gate U3(UINT target, double theta, double phi, double lambda) {
+    return internal::GateFactory::create_gate<internal::U3GateImpl>(target, theta, phi, lambda);
+}
+Gate CNOT(UINT control, UINT target) {
+    return internal::GateFactory::create_gate<internal::CNOTGateImpl>(control, target);
+}
+Gate CZ(UINT control, UINT target) {
+    return internal::GateFactory::create_gate<internal::CZGateImpl>(control, target);
+}
+Gate SWAP(UINT target1, UINT target2) {
+    return internal::GateFactory::create_gate<internal::SWAPGateImpl>(target1, target2);
+}
+Gate FusedSWAP(UINT qubit_index1, UINT qubit_index2, UINT block_size) {
+    return internal::GateFactory::create_gate<internal::FusedSWAPGateImpl>(
+        qubit_index1, qubit_index2, block_size);
+}
+}  // namespace qulacs

qulacs/gate/gate_npair_qubit.cpp

@@ -0,0 +1,11 @@
+#include "gate_npair_qubit.hpp"
+
+#include "update_ops.hpp"
+
+namespace qulacs {
+namespace internal {
+void FusedSWAPGateImpl::update_quantum_state(StateVector& state_vector) const {
+    fusedswap_gate(this->_qubit_index1, this->_qubit_index2, this->_block_size, state_vector);
+}
+}  // namespace internal
+}  // namespace qulacs

qulacs/gate/gate_npair_qubit.hpp

@@ -0,0 +1,46 @@
+#pragma once
+
+#include <cassert>
+#include <ranges>
+
+#include "gate.hpp"
+
+namespace qulacs {
+namespace internal {
+class FusedSWAPGateImpl : public GateBase {
+    UINT _qubit_index1, _qubit_index2, _block_size;
+
+public:
+    FusedSWAPGateImpl(UINT qubit_index1, UINT qubit_index2, UINT block_size)
+        : _qubit_index1(qubit_index1), _qubit_index2(qubit_index2), _block_size(block_size) {
+        UINT upper_index = std::max(qubit_index1, qubit_index2);
+        UINT lower_index = std::min(qubit_index1, qubit_index2);
+        if (upper_index <= (lower_index + block_size - 1)) {
+            throw std::runtime_error(
+                "FusedSwap: upper index must be bigger than lower_index + block_size - 1");
+        }
+    };
+
+    UINT qubit_index1() const { return _qubit_index1; }
+    UINT qubit_index2() const { return _qubit_index2; }
+    UINT block_size() const { return _block_size; }
+
+    std::vector<UINT> get_target_qubit_list() const override {
+        std::vector<UINT> res(_block_size * 2);
+        std::ranges::copy(std::views::iota(_qubit_index1, _qubit_index1 + _block_size),
+                          res.begin());
+        std::ranges::copy(std::views::iota(_qubit_index2, _qubit_index2 + _block_size),
+                          res.begin() + _block_size);
+        return res;
+    }
+    std::vector<UINT> get_control_qubit_list() const override { return {}; }
+
+    Gate copy() const override { return std::make_shared<FusedSWAPGateImpl>(*this); }
+    Gate get_inverse() const override { return std::make_shared<FusedSWAPGateImpl>(*this); }
+
+    void update_quantum_state(StateVector& state_vector) const override;
+};
+}  // namespace internal
+
+using FusedSWAPGate = internal::GatePtr<internal::FusedSWAPGateImpl>;
+}  // namespace qulacs

qulacs/gate/gate_one_control_one_target.cpp

@@ -3,11 +3,13 @@
 #include "update_ops.hpp"
 
 namespace qulacs {
-void CNOT::update_quantum_state(StateVector& state_vector) const {
+namespace internal {
+void CNOTGateImpl::update_quantum_state(StateVector& state_vector) const {
     cnot_gate(this->_control, this->_target, state_vector);
 }
 
-void CZ::update_quantum_state(StateVector& state_vector) const {
+void CZGateImpl::update_quantum_state(StateVector& state_vector) const {
     cz_gate(this->_control, this->_target, state_vector);
 }
+}  // namespace internal
 }  // namespace qulacs