add pauligate, paulirotationgate

qulacs/CMakeLists.txt

@@ -5,12 +5,14 @@ target_sources(qulacs PRIVATE
     gate/gate_npair_qubit.cpp
     gate/gate_one_control_one_target.cpp
     gate/gate_one_qubit.cpp
+    gate/gate_pauli.cpp
     gate/gate_quantum_matrix.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_pauli.cpp
     gate/update_ops_quantum_matrix.cpp
     gate/update_ops_two_qubit.cpp
     operator/pauli_operator.cpp

qulacs/all.hpp

@@ -7,6 +7,7 @@
 #include "gate/gate_npair_qubit.hpp"
 #include "gate/gate_one_control_one_target.hpp"
 #include "gate/gate_one_qubit.hpp"
+#include "gate/gate_pauli.hpp"
 #include "gate/gate_quantum_matrix.hpp"
 #include "gate/gate_two_qubit.hpp"
 #include "gate/update_ops.hpp"
@@ -16,3 +17,4 @@
 #include "state/state_vector.hpp"
 #include "types.hpp"
 #include "util/random.hpp"
+#include "util/bit_vector.hpp"

qulacs/gate/constant.hpp → qulacs/constant.hpp

@@ -3,7 +3,7 @@
 #include <array>
 #include <numbers>
 
-#include "../types.hpp"
+#include "types.hpp"
 
 namespace qulacs {
 //! PI value
@@ -91,8 +91,8 @@ KOKKOS_INLINE_FUNCTION
 matrix_2_2 PROJ_1_MATRIX() { return {0, 0, 0, 1}; }
 //! complex values for exp(j * i*pi/4 )
 KOKKOS_INLINE_FUNCTION
-matrix_2_2 PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; }
+array_4 PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; }
 //! complex values for exp(-j * i*pi/4 )
 KOKKOS_INLINE_FUNCTION
-matrix_2_2 PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; }
+array_4 PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; }
 }  // namespace qulacs

qulacs/gate/gate_pauli.cpp

@@ -0,0 +1,18 @@
+#include "gate_pauli.hpp"
+
+#include "../operator/pauli_operator.hpp"
+#include "update_ops.hpp"
+
+namespace qulacs {
+namespace internal {
+
+void PauliGateImpl::update_quantum_state(StateVector& state_vector) const {
+    pauli_gate(this->_pauli, state_vector);
+}
+
+void PauliRotationGateImpl::update_quantum_state(StateVector& state_vector) const {
+    pauli_rotation_gate(this->_pauli, this->_angle, state_vector);
+}
+
+}  // namespace internal
+}  // namespace qulacs

qulacs/gate/gate_pauli.hpp

@@ -0,0 +1,55 @@
+#pragma once
+
+#include <vector>
+
+#include "../operator/pauli_operator.hpp"
+#include "gate.hpp"
+
+namespace qulacs {
+namespace internal {
+class PauliGateImpl : public GateBase {
+    PauliOperator* _pauli;
+    std::vector<UINT> _target_index_list, _pauli_id_list;
+
+public:
+    PauliGateImpl(PauliOperator* pauli) {
+        _pauli = pauli;
+        _target_index_list = _pauli->get_target_qubit_list();
+        _pauli_id_list = _pauli->get_pauli_id_list();
+    };
+
+    std::vector<UINT> get_target_qubit_list() const override { return _target_index_list; }
+    std::vector<UINT> get_control_qubit_list() const override { return {}; }
+
+    Gate copy() const override { return std::make_shared<PauliGateImpl>(*this); }
+    Gate get_inverse() const override { return std::make_shared<PauliGateImpl>(*this); }
+
+    void update_quantum_state(StateVector& state_vector) const override;
+};
+
+class PauliRotationGateImpl : public GateBase {
+    PauliOperator* _pauli;
+    double _angle;
+    std::vector<UINT> _target_index_list, _pauli_id_list;
+
+public:
+    PauliRotationGateImpl(PauliOperator* pauli, double angle) {
+        _pauli = pauli;
+        _angle = angle;
+        _target_index_list = _pauli->get_target_qubit_list();
+        _pauli_id_list = _pauli->get_pauli_id_list();
+    };
+
+    std::vector<UINT> get_target_qubit_list() const override { return _target_index_list; }
+    std::vector<UINT> get_control_qubit_list() const override { return {}; }
+
+    Gate copy() const override { return std::make_shared<PauliRotationGateImpl>(*this); }
+    Gate get_inverse() const override { return std::make_shared<PauliRotationGateImpl>(*this); }
+
+    void update_quantum_state(StateVector& state_vector) const override;
+};
+}  // namespace internal
+
+using PauliGate = internal::GatePtr<internal::PauliGateImpl>;
+using PauliRotationGate = internal::GatePtr<internal::PauliRotationGateImpl>;
+}  // namespace qulacs

qulacs/gate/update_ops.hpp

@@ -1,6 +1,7 @@
 
 #pragma once
 
+#include "../operator/pauli_operator.hpp"
 #include "../state/state_vector.hpp"
 #include "../types.hpp"
 
@@ -60,12 +61,7 @@ void fusedswap_gate(UINT target_qubit_index_0,
                     UINT block_size,
                     StateVector& state);
 
-void pauli_gate(std::vector<UINT> target_qubit_index_list,
-                std::vector<UINT> pauli_id_index,
-                StateVector& state);
+void pauli_gate(PauliOperator* pauli, StateVector& state);
 
-void pauli_rotation_gate(std::vector<UINT> target_qubit_index_list,
-                         std::vector<UINT> pauli_id_index,
-                         double angle,
-                         StateVector& state);
+void pauli_rotation_gate(PauliOperator* pauli, double angle, StateVector& state);
 }  // namespace qulacs

qulacs/gate/update_ops_pauli.cpp

@@ -0,0 +1,63 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_StdAlgorithms.hpp>
+
+#include "../constant.hpp"
+#include "../operator/pauli_operator.hpp"
+#include "../types.hpp"
+#include "update_ops.hpp"
+
+namespace qulacs {
+void pauli_gate(PauliOperator* pauli, StateVector& state) { pauli->apply_to_state(state); }
+
+void pauli_rotation_gate(PauliOperator* pauli, double angle, StateVector& state) {
+    auto [bit_flip_mask_vector, phase_flip_mask_vector] = pauli->get_XZ_mask_representation();
+    UINT bit_flip_mask = bit_flip_mask_vector.data_raw()[0];
+    UINT phase_flip_mask = phase_flip_mask_vector.data_raw()[0];
+    UINT global_phase_90_rot_count = std::popcount(bit_flip_mask & phase_flip_mask);
+    const double cosval = cos(angle / 2);
+    const double sinval = sin(angle / 2);
+    const Complex coef = pauli->get_coef();
+    const auto& amplitudes = state.amplitudes_raw();
+    if (bit_flip_mask == 0) {
+        Kokkos::parallel_for(
+            state.dim(), KOKKOS_LAMBDA(const UINT& state_idx) {
+                if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) {
+                    amplitudes[state_idx] *= cosval - Complex(0, 1) * sinval;
+                } else {
+                    amplitudes[state_idx] *= cosval + Complex(0, 1) * sinval;
+                }
+                amplitudes[state_idx] *= coef;
+            });
+        return;
+    } else {
+        const UINT mask = 1 << bit_flip_mask_vector.msb();
+        const UINT mask_low = mask - 1;
+        const UINT mask_high = ~mask_low;
+        Kokkos::parallel_for(
+            state.dim(), KOKKOS_LAMBDA(const UINT& state_idx) {
+                UINT basis_0 = (state_idx & mask_low) + ((state_idx & mask_high) << 1);
+                UINT basis_1 = basis_0 ^ bit_flip_mask;
+
+                int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) % 2;
+                int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) % 2;
+
+                // fetch values
+                Complex cval_0 = amplitudes[basis_0];
+                Complex cval_1 = amplitudes[basis_1];
+
+                // set values
+                amplitudes[basis_0] =
+                    cosval * cval_0 +
+                    Complex(0, 1) * sinval * cval_1 *
+                        (PHASE_M90ROT()).val[(global_phase_90_rot_count + bit_parity_0 * 2) % 4];
+                amplitudes[basis_1] =
+                    cosval * cval_1 +
+                    Complex(0, 1) * sinval * cval_0 *
+                        (PHASE_M90ROT()).val[(global_phase_90_rot_count + bit_parity_1 * 2) % 4];
+                amplitudes[basis_0] *= coef;
+                amplitudes[basis_1] *= coef;
+            });
+    }
+}
+
+}  // namespace qulacs

qulacs/gate/update_ops_two_qubit.cpp

@@ -2,7 +2,7 @@
 #include <Kokkos_StdAlgorithms.hpp>
 
 #include "../types.hpp"
-#include "constant.hpp"
+#include "../constant.hpp"
 #include "update_ops.hpp"
 
 namespace qulacs {

qulacs/operator/pauli_operator.cpp

@@ -1,6 +1,6 @@
 #include "pauli_operator.hpp"
 
-#include "constant.hpp"
+#include "../constant.hpp"
 
 namespace qulacs {
 PauliOperator::PauliOperator(Complex coef) : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) {}

qulacs/util/utility.hpp

@@ -46,4 +46,14 @@ static std::string _check_gt(
                          << "), actual: " << val1 << " vs " << val2 << "\n";
     return error_message_stream.str();
 }
-};  // namespace qulacs
+
+inline static UINT count_population(UINT x) {
+    x = ((x & 0xaaaaaaaaaaaaaaaaUL) >> 1) + (x & 0x5555555555555555UL);
+    x = ((x & 0xccccccccccccccccUL) >> 2) + (x & 0x3333333333333333UL);
+    x = ((x & 0xf0f0f0f0f0f0f0f0UL) >> 4) + (x & 0x0f0f0f0f0f0f0f0fUL);
+    x = ((x & 0xff00ff00ff00ff00UL) >> 8) + (x & 0x00ff00ff00ff00ffUL);
+    x = ((x & 0xffff0000ffff0000UL) >> 16) + (x & 0x0000ffff0000ffffUL);
+    x = ((x & 0xffffffff00000000UL) >> 32) + (x & 0x00000000ffffffffUL);
+    return (UINT)x;
+}
+}  // namespace qulacs

tests/gate/gate_test.cpp

@@ -199,6 +199,51 @@ void run_random_gate_apply_fused(UINT n_qubits, UINT target0, UINT target1, UINT
     }
 }
 
+void run_random_gate_apply_Pauli(
+    UINT n_qubits, std::function<Eigen::MatrixXcd(double, double, double)> matrix_factory) {
+    const int dim = 1ULL << n_qubits;
+    Random random;
+
+    for (int repeat = 0; repeat < 10; repeat++) {
+        auto state = StateVector::Haar_random_state(n_qubits);
+        auto state_cp = state.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();
+            if (gate_type == 0) {
+                theta = 0;
+                phi = 0;
+            } else if (gate_type == 1) {
+                theta = M_PI / 2;
+            }
+            const auto matrix = matrix_factory(theta, phi, lambda);
+            const UINT target = random.int64() % n_qubits;
+            Gate gate;
+            if (gate_type == 0) {
+                gate = U1(target, lambda);
+            } else if (gate_type == 1) {
+                gate = U2(target, phi, lambda);
+            } else {
+                gate = U3(target, theta, phi, lambda);
+            }
+            gate->update_quantum_state(state);
+            state_cp = state.amplitudes();
+
+            test_state =
+                get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state;
+
+            for (int i = 0; i < dim; i++) {
+                ASSERT_NEAR(std::abs((CComplex)state_cp[i] - test_state[i]), 0, eps);
+            }
+        }
+    }
+}
+
 TEST(GateTest, ApplyI) { run_random_gate_apply<I>(5, make_I); }
 TEST(GateTest, ApplyX) { run_random_gate_apply<X>(5, make_X); }
 TEST(GateTest, ApplyY) { run_random_gate_apply<Y>(5, make_Y); }