merge_gateのcontrol, DenseMatrix対応

.devcontainer/cpu/devcontainer.json

@@ -31,8 +31,8 @@
 				"C_Cpp.default.includePath": [
 					"scaluq",
 					"/usr/local/include/kokkos",
-					"build/_deps/eigen_fetch-src",
-					"build/_deps/googletest_fetch-src/googletest/include",
+					"build/_deps/eigen-src",
+					"build/_deps/googletest-src/googletest/include",
 					"~/.local/lib/python3.10/site-packages/nanobind/include",
 					"/usr/include/python3.10"
 				],

.devcontainer/gpu/devcontainer.json

@@ -33,8 +33,8 @@
 				"C_Cpp.default.includePath": [
 					"scaluq",
 					"/usr/local/include/kokkos",
-					"build/_deps/eigen_fetch-src",
-					"build/_deps/googletest_fetch-src/googletest/include",
+					"build/_deps/eigen-src",
+					"build/_deps/googletest-src/googletest/include",
 					"~/.local/lib/python3.10/site-packages/nanobind/include",
 					"/usr/include/python3.10"
 				],

CMakeLists.txt

@@ -126,7 +126,7 @@ if ((${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") OR (${CMAKE_CXX_COMPILER_ID} STREQ
     endif()
 
     # Debug options
-    if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+    if ((${CMAKE_SYSTEM_NAME} MATCHES "Darwin") OR CMAKE_CUDA_COMPILER)
         target_compile_options(scaluq PUBLIC $<IF:$<CONFIG:Debug>,-O0 -g,-O3>)
     else()
         target_compile_options(scaluq PUBLIC $<IF:$<CONFIG:Debug>,-O0 -g -fsanitize=address$<COMMA>undefined,-O3>)

exe/main.cpp

@@ -9,41 +9,11 @@ using namespace scaluq;
 using namespace std;
 
 void run() {
-    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;
-
-    Eigen::Matrix<StdComplex, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat(4, 4);
-    mat << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15;
-
-    auto dense_gate = gate::DenseMatrix({1, 3}, mat);
-    std::cout << dense_gate << std::endl;
-
-    auto sparse_gate = gate::SparseMatrix({2, 0}, mat.sparseView());
-    std::cout << sparse_gate << std::endl;
+    internal::SparseComplexMatrix sparse(2, 2);
+    sparse.insert(0, 0) = 2;
+    auto g = gate::SparseMatrix({0}, sparse);
+    auto mat = g->get_matrix();
+    std::cout << mat << endl;
 }
 
 int main() {

pyproject.toml

@@ -62,7 +62,6 @@ ci = [
     "mypy == 1.11.2",
     "scikit-build == 0.17.6",
     "typing_extensions == 4.12.0",
-    "numpy == 1.26.0",
     "nanobind == 2.0.0"
 ]
 

scaluq/CMakeLists.txt

@@ -5,7 +5,7 @@ target_sources(scaluq PRIVATE
     gate/update_ops_dense_matrix.cpp
     gate/update_ops_sparse_matrix.cpp
     gate/update_ops_standard.cpp
-    # gate/merge_gate.cpp
+    gate/merge_gate.cpp
     operator/apply_pauli.cpp
     operator/pauli_operator.cpp
     operator/operator.cpp

scaluq/gate/merge_gate.cpp

@@ -5,19 +5,54 @@
 #include "gate/gate_factory.hpp"
 
 namespace scaluq {
-/**
- * @details ignore global phase (because PauliRoation -> matrix ignore global phase)
- */
+std::pair<Gate, double> merge_gate_dense_matrix(const Gate& gate1, const Gate& gate2) {
+    auto common_control_mask = gate1->control_qubit_mask() & gate2->control_qubit_mask();
+    auto merged_operand_mask =
+        (gate1->operand_qubit_mask() | gate2->operand_qubit_mask()) & ~common_control_mask;
+    auto merged_operand_vector = internal::mask_to_vector(merged_operand_mask);
+    auto matrix1 = internal::get_expanded_matrix(gate1->get_matrix(),
+                                                 gate1->target_qubit_list(),
+                                                 gate1->control_qubit_mask() & ~common_control_mask,
+                                                 merged_operand_vector);
+    auto matrix2 = internal::get_expanded_matrix(gate2->get_matrix(),
+                                                 gate2->target_qubit_list(),
+                                                 gate2->control_qubit_mask() & ~common_control_mask,
+                                                 merged_operand_vector);
+    auto matrix = matrix2 * matrix1;
+    return {gate::DenseMatrix(
+                merged_operand_vector, matrix, internal::mask_to_vector(common_control_mask)),
+            0.};
+}
+
 std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
     GateType gate_type1 = gate1.gate_type();
     GateType gate_type2 = gate2.gate_type();
-    // TODO: Deal with ProbablisticGate
 
-    // Special case: Zero qubit
-    if (gate_type1 == GateType::I) return {gate2, 0.};  // copy can be removed by #125
+    if (gate_type1 == GateType::Probablistic || gate_type2 == GateType::Probablistic) {
+        throw std::runtime_error(
+            "merge_gate(const Gate&, const Gate&): ProbablisticGate is not supported.");
+    }
+
+    if (gate_type1 == GateType::I) return {gate2, 0.};
     if (gate_type2 == GateType::I) return {gate1, 0.};
-    if (gate_type1 == GateType::GlobalPhase) return {gate2, GlobalPhaseGate(gate1)->phase()};
-    if (gate_type2 == GateType::GlobalPhase) return {gate1, GlobalPhaseGate(gate2)->phase()};
+
+    auto gate1_control_mask = gate1->control_qubit_mask();
+    auto gate2_control_mask = gate2->control_qubit_mask();
+    if (gate1_control_mask != gate2_control_mask) return merge_gate_dense_matrix(gate1, gate2);
+    auto control_list = internal::mask_to_vector(gate1_control_mask);
+
+    // Special case: Zero qubit
+    if (gate_type1 == GateType::GlobalPhase || gate_type2 == GateType::GlobalPhase) {
+        const auto& [phase, gate] = [&] {
+            if (gate_type1 == GateType::GlobalPhase)
+                return std::make_pair(GlobalPhaseGate(gate1)->phase(), gate2);
+            return std::make_pair(GlobalPhaseGate(gate2)->phase(), gate1);
+        }();
+        if (gate1_control_mask == 0) return {gate, phase};
+        auto gate_type = gate.gate_type();
+        if (gate_type == GateType::GlobalPhase)
+            return {gate::GlobalPhase(phase + GlobalPhaseGate(gate)->phase(), control_list), 0.};
+    }
 
     // Special case: Pauli
     auto get_pauli_id = [&](GateType gate_type) -> std::optional<std::uint64_t> {
@@ -37,16 +72,40 @@ std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
         if (target1 == target2) {
             if (pauli_id1 == pauli_id2) return {gate::I(), 0.};
             if (pauli_id1 == 1) {
-                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_id2 == 2) {
+                    if (gate1_control_mask == 0) {
+                        return {gate::Z(target1, control_list), -Kokkos::numbers::pi / 2};
+                    }
+                }
+                if (pauli_id2 == 3) {
+                    if (gate1_control_mask == 0) {
+                        return {gate::Y(target1, control_list), Kokkos::numbers::pi / 2};
+                    }
+                }
             }
             if (pauli_id1 == 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_id2 == 3) {
+                    if (gate1_control_mask == 0) {
+                        return {gate::X(target1, control_list), -Kokkos::numbers::pi / 2};
+                    }
+                }
+                if (pauli_id2 == 1) {
+                    if (gate1_control_mask == 0) {
+                        return {gate::Z(target1, control_list), Kokkos::numbers::pi / 2};
+                    }
+                }
             }
             if (pauli_id1 == 3) {
-                if (pauli_id2 == 1) return {gate::Y(target1), -Kokkos::numbers::pi / 2};
-                if (pauli_id2 == 2) return {gate::X(target1), Kokkos::numbers::pi / 2};
+                if (pauli_id2 == 1) {
+                    if (gate1_control_mask == 0) {
+                        return {gate::Y(target1, control_list), -Kokkos::numbers::pi / 2};
+                    }
+                }
+                if (pauli_id2 == 2) {
+                    if (gate1_control_mask == 0) {
+                        return {gate::X(target1, control_list), Kokkos::numbers::pi / 2};
+                    }
+                }
             }
         }
     }
@@ -60,9 +119,11 @@ std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
                           ? PauliGate(gate2)->pauli()
                           : PauliOperator(std::vector{gate2->target_qubit_list()[0]},
                                           std::vector{pauli_id2.value()});
-        return {gate::Pauli(pauli2 * pauli1), 0.};
+        return {gate::Pauli(pauli2 * pauli1, control_list), 0.};
     }
 
+    constexpr double eps = 1e-12;
+
     // Special case: Phase
     auto get_oct_phase = [&](GateType gate_type) -> std::optional<std::uint64_t> {
         if (gate_type == GateType::I) return 0;
@@ -77,11 +138,11 @@ std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
                               std::uint64_t target) -> std::optional<Gate> {
         oct_phase &= 7;
         if (oct_phase == 0) return gate::I();
-        if (oct_phase == 4) return gate::Z(target);
-        if (oct_phase == 2) return gate::S(target);
-        if (oct_phase == 6) return gate::Sdag(target);
-        if (oct_phase == 1) return gate::T(target);
-        if (oct_phase == 7) return gate::Tdag(target);
+        if (oct_phase == 4) return gate::Z(target, control_list);
+        if (oct_phase == 2) return gate::S(target, control_list);
+        if (oct_phase == 6) return gate::Sdag(target, control_list);
+        if (oct_phase == 1) return gate::T(target, control_list);
+        if (oct_phase == 7) return gate::Tdag(target, control_list);
         return std::nullopt;
     };
     auto oct_phase1 = get_oct_phase(gate_type1);
@@ -107,7 +168,11 @@ std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
                             : gate_type2 == GateType::RZ ? RZGate(gate2)->angle()
                                                          : U1Gate(gate2)->lambda();
             double global_phase2 = gate_type2 == GateType::RZ ? -RZGate(gate2)->angle() / 2 : 0.;
-            return {gate::U1(target1, phase1 + phase2), global_phase1 + global_phase2};
+            double global_phase = global_phase1 + global_phase2;
+            if (std::abs(global_phase) < eps) {
+                return {gate::U1(target1, phase1 + phase2, control_list),
+                        global_phase1 + global_phase2};
+            }
         }
     }
 
@@ -127,9 +192,12 @@ std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
         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;
+        double global_phase = global_phase1 + global_phase2;
         if (target1 == target2) {
-            return {gate::RX(target1, rx_param1.value() + rx_param2.value()),
-                    global_phase1 + global_phase2};
+            if (std::abs(global_phase) < eps) {
+                return {gate::RX(target1, rx_param1.value() + rx_param2.value(), control_list),
+                        global_phase1 + global_phase2};
+            }
         }
     }
 
@@ -149,64 +217,21 @@ std::pair<Gate, double> merge_gate(const Gate& gate1, const Gate& gate2) {
         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;
+        double global_phase = global_phase1 + global_phase2;
         if (target1 == target2) {
-            return {gate::RY(target1, ry_param1.value() + ry_param2.value()),
-                    global_phase1 + global_phase2};
+            if (std::abs(global_phase) < eps) {
+                return {gate::RY(target1, ry_param1.value() + ry_param2.value(), control_list),
+                        global_phase1 + global_phase2};
+            }
         }
     }
 
-    // Special case: CX,CZ,Swap,FusedSwap duplication
-    if (gate_type1 == gate_type2 && gate_type1 == GateType::CX) {
-        CXGate cx1(gate1), cx2(gate2);
-        if (cx1->target() == cx2->target() && cx1->control() == cx2->control())
-            return {gate::I(), 0.};
-    }
-    if (gate_type1 == gate_type2 && gate_type1 == GateType::CZ) {
-        CZGate cz1(gate1), cz2(gate2);
-        if (cz1->target() == cz2->target() && cz1->control() == cz2->control())
-            return {gate::I(), 0.};
-        if (cz1->target() == cz2->control() && cz1->control() == cz2->target())
-            return {gate::I(), 0.};
-    }
+    // Special case: Swap duplication
     if (gate_type1 == gate_type2 && gate_type1 == GateType::Swap) {
-        SwapGate swap1(gate1), swap2(gate2);
-        if (swap1->target1() == swap2->target1() && swap1->target2() == swap2->target2())
-            return {gate::I(), 0.};
-        if (swap1->target1() == swap2->target2() && swap1->target2() == swap2->target1())
-            return {gate::I(), 0.};
-    }
-    if (gate_type1 == gate_type2 && gate_type1 == GateType::FusedSwap) {
-        FusedSwapGate swap1(gate1), swap2(gate2);
-        if (swap1->block_size() == swap2->block_size()) {
-            if (swap1->qubit_index1() == swap2->qubit_index1() &&
-                swap1->qubit_index2() == swap2->qubit_index2())
-                return {gate::I(), 0.};
-            if (swap1->qubit_index1() == swap2->qubit_index2() &&
-                swap1->qubit_index2() == swap2->qubit_index1())
-                return {gate::I(), 0.};
-        }
+        if (gate1->target_qubit_mask() == gate2->target_qubit_mask()) return {gate::I(), 0.};
     }
 
     // General case
-    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<std::uint64_t> 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());
-    std::ranges::sort(merged_targets);
-    merged_targets.erase(std::ranges::unique(merged_targets).begin(), merged_targets.end());
-    if (merged_targets.size() >= 3) {
-        throw std::runtime_error(
-            "gate::merge_gate: Result gate's target size is equal or more than three. This is "
-            "currently not implemented.");
-    }
-    auto matrix1 =
-        internal::get_expanded_matrix(gate1->get_matrix().value(), gate1_targets, merged_targets);
-    auto matrix2 =
-        internal::get_expanded_matrix(gate2->get_matrix().value(), gate2_targets, merged_targets);
-    auto matrix = matrix2 * matrix1;
-    return {gate::DenseMatrix(merged_targets, matrix), 0.};
+    return merge_gate_dense_matrix(gate1, gate2);
 }
 }  // namespace scaluq

scaluq/util/random.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <random>
+#include <ranges>
 
 #include "../types.hpp"
 
@@ -20,6 +21,16 @@ class Random {
 
     [[nodiscard]] std::uint64_t int64() { return this->mt(); }
 
-    [[nodiscard]] std::uint32_t int32() { return this->mt() % UINT32_MAX; }
+    [[nodiscard]] std::uint32_t int32() { return this->mt() & UINT32_MAX; }
+
+    [[nodiscard]] std::vector<std::uint64_t> permutation(std::uint64_t n) {
+        std::vector<std::uint64_t> 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 = int32() % (i + 1);
+            if (i != j) std::swap(shuffled[i], shuffled[j]);
+        }
+        return shuffled;
+    }
 };
 }  // namespace scaluq