diff --git a/scaluq/constant.hpp b/scaluq/constant.hpp
index cc415f5b..af4fdbcf 100644
--- a/scaluq/constant.hpp
+++ b/scaluq/constant.hpp
@@ -9,11 +9,6 @@ namespace scaluq {
 
 namespace internal {
 
-#define DEF_MATH_CONSTANT(TRAIT, VALUE)                                                        \
-    template <class T>                                                                         \
-    inline constexpr auto TRAIT##_v = std::enable_if_t<std::is_floating_point_v<T>, T>(VALUE); \
-    inline constexpr auto TRAIT = TRAIT##_v<double>
-
 //! inverse square root of 2
 KOKKOS_INLINE_FUNCTION
 double INVERSE_SQRT2() { return Kokkos::numbers::sqrt2 / 2; }
@@ -27,70 +22,106 @@ KOKKOS_INLINE_FUNCTION
 double SINPI8() { return 0.382683432365090; }
 
 //! identity matrix
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 I_GATE() { return {1, 0, 0, 1}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> I_GATE() {
+    return {1, 0, 0, 1};
+}
 //! Pauli matrix X
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 X_GATE() { return {0, 1, 1, 0}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> X_GATE() {
+    return {0, 1, 1, 0};
+}
 //! Pauli matrix Y
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 Y_GATE() { return {0, Complex(0, -1), Complex(0, 1), 0}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> Y_GATE() {
+    return {0, Complex<Fp>(0, -1), Complex<Fp>(0, 1), 0};
+}
 //! Pauli matrix Z
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 Z_GATE() { return {1, 0, 0, -1}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> Z_GATE() {
+    return {1, 0, 0, -1};
+}
 
 //! list of Pauli matrix I,X,Y,Z
 // std::array<Matrix2x2, 4> PAULI_MATRIX = {I_GATE, X_GATE,
 // Y_GATE, Z_GATE};
 
 //! S-gate
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 S_GATE_MATRIX() { return {1, 0, 0, Complex(0, 1)}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> S_GATE_MATRIX() {
+    return {1, 0, 0, Complex<Fp>(0, 1)};
+}
 //! Sdag-gate
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 S_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(0, -1)}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> S_DAG_GATE_MATRIX() {
+    return {1, 0, 0, Complex<Fp>(0, -1)};
+}
 //! T-gate
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 T_GATE_MATRIX() { return {1, 0, 0, Complex(INVERSE_SQRT2(), INVERSE_SQRT2())}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> T_GATE_MATRIX() {
+    return {1, 0, 0, Complex<Fp>(INVERSE_SQRT2(), INVERSE_SQRT2())};
+}
 //! Tdag-gate
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 T_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(INVERSE_SQRT2(), -INVERSE_SQRT2())}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> T_DAG_GATE_MATRIX() {
+    return {1, 0, 0, Complex<Fp>(INVERSE_SQRT2(), -INVERSE_SQRT2())};
+}
 //! Hadamard gate
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 HADAMARD_MATRIX() {
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> HADAMARD_MATRIX() {
     return {INVERSE_SQRT2(), INVERSE_SQRT2(), INVERSE_SQRT2(), -INVERSE_SQRT2()};
 }
 //! square root of X gate
-KOKKOS_INLINE_FUNCTION
-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)};
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> SQRT_X_GATE_MATRIX() {
+    return {Complex<Fp>(0.5, 0.5),
+            Complex<Fp>(0.5, -0.5),
+            Complex<Fp>(0.5, -0.5),
+            Complex<Fp>(0.5, 0.5)};
 }
 //! square root of Y gate
-KOKKOS_INLINE_FUNCTION
-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)};
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> SQRT_Y_GATE_MATRIX() {
+    return {Complex<Fp>(0.5, 0.5),
+            Complex<Fp>(-0.5, -0.5),
+            Complex<Fp>(0.5, 0.5),
+            Complex<Fp>(0.5, 0.5)};
 }
 //! square root dagger of X gate
-KOKKOS_INLINE_FUNCTION
-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)};
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> SQRT_X_DAG_GATE_MATRIX() {
+    return {Complex<Fp>(0.5, -0.5),
+            Complex<Fp>(0.5, 0.5),
+            Complex<Fp>(0.5, 0.5),
+            Complex<Fp>(0.5, -0.5)};
 }
 //! square root dagger of Y gate
-KOKKOS_INLINE_FUNCTION
-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)};
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> SQRT_Y_DAG_GATE_MATRIX() {
+    return {Complex<Fp>(0.5, -0.5),
+            Complex<Fp>(0.5, -0.5),
+            Complex<Fp>(-0.5, 0.5),
+            Complex<Fp>(0.5, -0.5)};
 }
 //! Projection to 0
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 PROJ_0_MATRIX() { return {1, 0, 0, 0}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> PROJ_0_MATRIX() {
+    return {1, 0, 0, 0};
+}
 //! Projection to 1
-KOKKOS_INLINE_FUNCTION
-Matrix2x2 PROJ_1_MATRIX() { return {0, 0, 0, 1}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> PROJ_1_MATRIX() {
+    return {0, 0, 0, 1};
+}
 //! complex values for exp(j * i*pi/4 )
-KOKKOS_INLINE_FUNCTION
-Kokkos::Array<Complex, 4> PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Kokkos::Array<Complex<Fp>, 4> PHASE_90ROT() {
+    return {1., Complex<Fp>(0, 1), -1, Complex<Fp>(0, -1)};
+}
 //! complex values for exp(-j * i*pi/4 )
-KOKKOS_INLINE_FUNCTION
-Kokkos::Array<Complex, 4> PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; }
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Kokkos::Array<Complex<Fp>, 4> PHASE_M90ROT() {
+    return {1., Complex<Fp>(0, -1), -1, Complex<Fp>(0, 1)};
+}
 }  // namespace internal
 }  // namespace scaluq
diff --git a/scaluq/gate/gate.hpp b/scaluq/gate/gate.hpp
index 259b4156..bf652d7f 100644
--- a/scaluq/gate/gate.hpp
+++ b/scaluq/gate/gate.hpp
@@ -229,7 +229,7 @@ class GateBase : public std::enable_shared_from_this<GateBase<_FloatType>> {
     }
 
     [[nodiscard]] virtual std::shared_ptr<const GateBase<Fp>> get_inverse() const = 0;
-    [[nodiscard]] virtual internal::ComplexMatrix get_matrix() const = 0;
+    [[nodiscard]] virtual internal::ComplexMatrix<Fp> get_matrix() const = 0;
 
     virtual void update_quantum_state(StateVector<Fp>& state_vector) const = 0;
 
diff --git a/scaluq/gate/gate_factory.hpp b/scaluq/gate/gate_factory.hpp
index 68e244c2..e3d39a7a 100644
--- a/scaluq/gate/gate_factory.hpp
+++ b/scaluq/gate/gate_factory.hpp
@@ -143,7 +143,7 @@ inline Gate<Fp> U3(std::uint64_t target,
 }
 template <std::floating_point Fp>
 inline Gate<Fp> OneTargetMatrix(std::uint64_t target,
-                                const std::array<std::array<Complex, 2>, 2>& matrix,
+                                const std::array<std::array<Complex<Fp>, 2>, 2>& matrix,
                                 const std::vector<std::uint64_t>& controls = {}) {
     return internal::GateFactory::create_gate<internal::OneTargetMatrixGateImpl<Fp>>(
         internal::vector_to_mask({target}), internal::vector_to_mask(controls), matrix);
@@ -179,7 +179,7 @@ inline Gate<Fp> Swap(std::uint64_t target1,
 template <std::floating_point Fp>
 inline Gate<Fp> TwoTargetMatrix(std::uint64_t target1,
                                 std::uint64_t target2,
-                                const std::array<std::array<Complex, 4>, 4>& matrix,
+                                const std::array<std::array<Complex<Fp>, 4>, 4>& matrix,
                                 const std::vector<std::uint64_t>& controls = {}) {
     return internal::GateFactory::create_gate<internal::TwoTargetMatrixGateImpl<Fp>>(
         internal::vector_to_mask({target1, target2}), internal::vector_to_mask(controls), matrix);
@@ -200,47 +200,49 @@ inline Gate<Fp> PauliRotation(const PauliOperator<Fp>& pauli,
 }
 template <std::floating_point Fp>
 inline Gate<Fp> DenseMatrix(const std::vector<std::uint64_t>& targets,
-                            const internal::ComplexMatrix& matrix,
+                            const internal::ComplexMatrix<Fp>& matrix,
                             const std::vector<std::uint64_t>& controls = {}) {
     std::uint64_t nqubits = targets.size();
     std::uint64_t dim = 1ULL << nqubits;
     if (static_cast<std::uint64_t>(matrix.rows()) != dim ||
         static_cast<std::uint64_t>(matrix.cols()) != dim) {
         throw std::runtime_error(
-            "gate::DenseMatrix(const std::vector<std::uint64_t>&, const internal::ComplexMatrix&): "
+            "gate::DenseMatrix(const std::vector<std::uint64_t>&, const "
+            "internal::ComplexMatrix<Fp>&): "
             "matrix size must be 2^{n_qubits} x 2^{n_qubits}.");
     }
     if (targets.size() == 0) return I<Fp>();
     if (targets.size() == 1) {
         return OneTargetMatrix<Fp>(
             targets[0],
-            std::array{std::array{Complex(matrix(0, 0)), Complex(matrix(0, 1))},
-                       std::array{Complex(matrix(1, 0)), Complex(matrix(1, 1))}},
+            std::array{std::array{Complex<Fp>(matrix(0, 0)), Complex<Fp>(matrix(0, 1))},
+                       std::array{Complex<Fp>(matrix(1, 0)), Complex<Fp>(matrix(1, 1))}},
             controls);
     }
     if (targets.size() == 2) {
         return TwoTargetMatrix<Fp>(targets[0],
                                    targets[1],
-                                   std::array{std::array{Complex(matrix(0, 0)),
-                                                         Complex(matrix(0, 1)),
-                                                         Complex(matrix(0, 2)),
-                                                         Complex(matrix(0, 3))},
-                                              std::array{Complex(matrix(1, 0)),
-                                                         Complex(matrix(1, 1)),
-                                                         Complex(matrix(1, 2)),
-                                                         Complex(matrix(1, 3))},
-                                              std::array{Complex(matrix(2, 0)),
-                                                         Complex(matrix(2, 1)),
-                                                         Complex(matrix(2, 2)),
-                                                         Complex(matrix(2, 3))},
-                                              std::array{Complex(matrix(3, 0)),
-                                                         Complex(matrix(3, 1)),
-                                                         Complex(matrix(3, 2)),
-                                                         Complex(matrix(3, 3))}},
+                                   std::array{std::array{Complex<Fp>(matrix(0, 0)),
+                                                         Complex<Fp>(matrix(0, 1)),
+                                                         Complex<Fp>(matrix(0, 2)),
+                                                         Complex<Fp>(matrix(0, 3))},
+                                              std::array{Complex<Fp>(matrix(1, 0)),
+                                                         Complex<Fp>(matrix(1, 1)),
+                                                         Complex<Fp>(matrix(1, 2)),
+                                                         Complex<Fp>(matrix(1, 3))},
+                                              std::array{Complex<Fp>(matrix(2, 0)),
+                                                         Complex<Fp>(matrix(2, 1)),
+                                                         Complex<Fp>(matrix(2, 2)),
+                                                         Complex<Fp>(matrix(2, 3))},
+                                              std::array{Complex<Fp>(matrix(3, 0)),
+                                                         Complex<Fp>(matrix(3, 1)),
+                                                         Complex<Fp>(matrix(3, 2)),
+                                                         Complex<Fp>(matrix(3, 3))}},
                                    controls);
     }
     throw std::runtime_error(
-        "gate::DenseMatrix(const std::vector<std::uint64_t>&, const internal::ComplexMatrix&): "
+        "gate::DenseMatrix(const std::vector<std::uint64_t>&, const "
+        "internal::ComplexMatrix<Fp>&): "
         "DenseMatrix "
         "gate "
         "more "
diff --git a/scaluq/gate/gate_matrix.hpp b/scaluq/gate/gate_matrix.hpp
index 8a82c315..a8f82d70 100644
--- a/scaluq/gate/gate_matrix.hpp
+++ b/scaluq/gate/gate_matrix.hpp
@@ -13,12 +13,12 @@ namespace internal {
 
 template <std::floating_point Fp>
 class OneTargetMatrixGateImpl : public GateBase<Fp> {
-    Matrix2x2 _matrix;
+    Matrix2x2<Fp> _matrix;
 
 public:
     OneTargetMatrixGateImpl(std::uint64_t target_mask,
                             std::uint64_t control_mask,
-                            const std::array<std::array<Complex, 2>, 2>& matrix)
+                            const std::array<std::array<Complex<Fp>, 2>, 2>& matrix)
         : GateBase<Fp>(target_mask, control_mask) {
         _matrix[0][0] = matrix[0][0];
         _matrix[0][1] = matrix[0][1];
@@ -26,7 +26,7 @@ class OneTargetMatrixGateImpl : public GateBase<Fp> {
         _matrix[1][1] = matrix[1][1];
     }
 
-    std::array<std::array<Complex, 2>, 2> matrix() const {
+    std::array<std::array<Complex<Fp>, 2>, 2> matrix() const {
         return {_matrix[0][0], _matrix[0][1], _matrix[1][0], _matrix[1][1]};
     }
 
@@ -34,13 +34,13 @@ class OneTargetMatrixGateImpl : public GateBase<Fp> {
         return std::make_shared<const OneTargetMatrixGateImpl>(
             this->_target_mask,
             this->_control_mask,
-            std::array<std::array<Complex, 2>, 2>{Kokkos::conj(_matrix[0][0]),
-                                                  Kokkos::conj(_matrix[1][0]),
-                                                  Kokkos::conj(_matrix[0][1]),
-                                                  Kokkos::conj(_matrix[1][1])});
+            std::array<std::array<Complex<Fp>, 2>, 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);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << this->_matrix[0][0], this->_matrix[0][1], this->_matrix[1][0], this->_matrix[1][1];
         return mat;
     }
@@ -61,12 +61,12 @@ class OneTargetMatrixGateImpl : public GateBase<Fp> {
 
 template <std::floating_point Fp>
 class TwoTargetMatrixGateImpl : public GateBase<Fp> {
-    Matrix4x4 _matrix;
+    Matrix4x4<Fp> _matrix;
 
 public:
     TwoTargetMatrixGateImpl(std::uint64_t target_mask,
                             std::uint64_t control_mask,
-                            const std::array<std::array<Complex, 4>, 4>& matrix)
+                            const std::array<std::array<Complex<Fp>, 4>, 4>& matrix)
         : GateBase<Fp>(target_mask, control_mask) {
         for (std::uint64_t i : std::views::iota(0, 4)) {
             for (std::uint64_t j : std::views::iota(0, 4)) {
@@ -75,8 +75,8 @@ class TwoTargetMatrixGateImpl : public GateBase<Fp> {
         }
     }
 
-    std::array<std::array<Complex, 4>, 4> matrix() const {
-        std::array<std::array<Complex, 4>, 4> matrix;
+    std::array<std::array<Complex<Fp>, 4>, 4> matrix() const {
+        std::array<std::array<Complex<Fp>, 4>, 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[i][j];
@@ -86,7 +86,7 @@ class TwoTargetMatrixGateImpl : public GateBase<Fp> {
     }
 
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
-        std::array<std::array<Complex, 4>, 4> matrix_dag;
+        std::array<std::array<Complex<Fp>, 4>, 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[j][i]);
@@ -95,8 +95,8 @@ class TwoTargetMatrixGateImpl : public GateBase<Fp> {
         return std::make_shared<const TwoTargetMatrixGateImpl>(
             this->_target_mask, this->_control_mask, matrix_dag);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(4, 4);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> 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],
diff --git a/scaluq/gate/gate_pauli.hpp b/scaluq/gate/gate_pauli.hpp
index 73404de4..50fe2423 100644
--- a/scaluq/gate/gate_pauli.hpp
+++ b/scaluq/gate/gate_pauli.hpp
@@ -25,7 +25,7 @@ class PauliGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override { return this->_pauli.get_matrix(); }
+    internal::ComplexMatrix<Fp> get_matrix() const override { return this->_pauli.get_matrix(); }
 
     void update_quantum_state(StateVector<Fp>& state_vector) const override {
         auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation();
@@ -66,13 +66,13 @@ class PauliRotationGateImpl : public GateBase<Fp> {
             this->_control_mask, _pauli, -_angle);
     }
 
-    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) *
-                  internal::ComplexMatrix::Identity(mat.rows(), mat.cols()) +
-              imag_unit * (StdComplex)Kokkos::sin(-true_angle / 2) * mat;
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat = this->_pauli.get_matrix_ignoring_coef();
+        Complex<Fp> true_angle = _angle * _pauli.coef();
+        StdComplex<Fp> imag_unit(0, 1);
+        mat = (StdComplex<Fp>)Kokkos::cos(-true_angle / 2) *
+                  internal::ComplexMatrix<Fp>::Identity(mat.rows(), mat.cols()) +
+              imag_unit * (StdComplex<Fp>)Kokkos::sin(-true_angle / 2) * mat;
         return mat;
     }
     void update_quantum_state(StateVector<Fp>& state_vector) const override {
diff --git a/scaluq/gate/gate_probablistic.hpp b/scaluq/gate/gate_probablistic.hpp
index 05f28078..8cf51f30 100644
--- a/scaluq/gate/gate_probablistic.hpp
+++ b/scaluq/gate/gate_probablistic.hpp
@@ -73,7 +73,7 @@ class ProbablisticGateImpl : public GateBase<Fp> {
             [](const std::shared_ptr<const GateBase<Fp>>& gate) { return gate->get_inverse(); });
         return std::make_shared<const ProbablisticGateImpl>(_distribution, inv_gate_list);
     }
-    internal::ComplexMatrix get_matrix() const override {
+    internal::ComplexMatrix<Fp> get_matrix() const override {
         throw std::runtime_error(
             "ProbablisticGateImpl::get_matrix(): This function must not be used in "
             "ProbablisticGateImpl.");
diff --git a/scaluq/gate/gate_standard.hpp b/scaluq/gate/gate_standard.hpp
index ba6b7c48..a692d8b4 100644
--- a/scaluq/gate/gate_standard.hpp
+++ b/scaluq/gate/gate_standard.hpp
@@ -14,8 +14,8 @@ class IGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override {
-        return internal::ComplexMatrix::Identity(1, 1);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        return internal::ComplexMatrix<Fp>::Identity(1, 1);
     }
 
     void update_quantum_state(StateVector<Fp>& state_vector) const override {
@@ -44,8 +44,8 @@ class GlobalPhaseGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const GlobalPhaseGateImpl<Fp>>(this->_control_mask, -_phase);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        return internal::ComplexMatrix::Identity(1, 1) * std::exp(std::complex<Fp>(0, _phase));
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        return internal::ComplexMatrix<Fp>::Identity(1, 1) * std::exp(std::complex<Fp>(0, _phase));
     }
 
     void update_quantum_state(StateVector<Fp>& state_vector) const override {
@@ -82,8 +82,8 @@ class XGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 0, 1, 1, 0;
         return mat;
     }
@@ -109,8 +109,8 @@ class YGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 0, -1i, 1i, 0;
         return mat;
     }
@@ -136,8 +136,8 @@ class ZGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, -1;
         return mat;
     }
@@ -163,8 +163,8 @@ class HGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 1, 1, -1;
         mat /= std::sqrt(2);
         return mat;
@@ -208,8 +208,8 @@ class SGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const SdagGateImpl<Fp>>(this->_target_mask, this->_control_mask);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, 1i;
         return mat;
     }
@@ -235,8 +235,8 @@ class SdagGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const SGateImpl<Fp>>(this->_target_mask, this->_control_mask);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, -1i;
         return mat;
     }
@@ -262,8 +262,8 @@ class TGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const TdagGateImpl<Fp>>(this->_target_mask, this->_control_mask);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, (1. + 1i) / std::sqrt(2);
         return mat;
     }
@@ -289,8 +289,8 @@ class TdagGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const TGateImpl<Fp>>(this->_target_mask, this->_control_mask);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, (1. - 1.i) / std::sqrt(2);
         return mat;
     }
@@ -318,9 +318,10 @@ class SqrtXGateImpl : public GateBase<Fp> {
                                                             this->_control_mask);
     }
 
-    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;
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
+        mat << StdComplex<Fp>(0.5, 0.5), StdComplex<Fp>(0.5, -0.5), StdComplex<Fp>(0.5, -0.5),
+            StdComplex<Fp>(0.5, 0.5);
         return mat;
     }
 
@@ -345,9 +346,10 @@ class SqrtXdagGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const SqrtXGateImpl<Fp>>(this->_target_mask, this->_control_mask);
     }
-    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;
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
+        mat << StdComplex<Fp>(0.5, -0.5), StdComplex<Fp>(0.5, 0.5), StdComplex<Fp>(0.5, 0.5),
+            StdComplex<Fp>(0.5, -0.5);
         return mat;
     }
 
@@ -374,8 +376,8 @@ class SqrtYGateImpl : public GateBase<Fp> {
                                                             this->_control_mask);
     }
 
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 0.5 + 0.5i, -0.5 - 0.5i, 0.5 + 0.5i, 0.5 + 0.5i;
         return mat;
     }
@@ -401,8 +403,8 @@ class SqrtYdagGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return std::make_shared<const SqrtYGateImpl<Fp>>(this->_target_mask, this->_control_mask);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 0.5 - 0.5i, 0.5 - 0.5i, -0.5 + 0.5i, 0.5 - 0.5i;
         return mat;
     }
@@ -428,8 +430,8 @@ class P0GateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         throw std::runtime_error("P0::get_inverse: Projection gate doesn't have inverse gate");
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, 0;
         return mat;
     }
@@ -455,8 +457,8 @@ class P1GateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         throw std::runtime_error("P1::get_inverse: Projection gate doesn't have inverse gate");
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 0, 0, 0, 1;
         return mat;
     }
@@ -483,8 +485,8 @@ class RXGateImpl : public RotationGateBase<Fp> {
         return std::make_shared<const RXGateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -this->_angle);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << std::cos(this->_angle / 2), -1i * std::sin(this->_angle / 2),
             -1i * std::sin(this->_angle / 2), std::cos(this->_angle / 2);
         return mat;
@@ -513,8 +515,8 @@ class RYGateImpl : public RotationGateBase<Fp> {
         return std::make_shared<const RYGateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -this->_angle);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << std::cos(this->_angle / 2), -std::sin(this->_angle / 2), std::sin(this->_angle / 2),
             std::cos(this->_angle / 2);
         return mat;
@@ -543,8 +545,8 @@ class RZGateImpl : public RotationGateBase<Fp> {
         return std::make_shared<const RZGateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -this->_angle);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << std::exp(-0.5i * this->_angle), 0, 0, std::exp(0.5i * this->_angle);
         return mat;
     }
@@ -577,8 +579,8 @@ class U1GateImpl : public GateBase<Fp> {
         return std::make_shared<const U1GateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -_lambda);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << 1, 0, 0, std::exp(1i * _lambda);
         return mat;
     }
@@ -612,8 +614,8 @@ class U2GateImpl : public GateBase<Fp> {
                                                       -_lambda - Kokkos::numbers::pi,
                                                       -_phi + Kokkos::numbers::pi);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> 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.),
@@ -650,8 +652,8 @@ class U3GateImpl : public GateBase<Fp> {
         return std::make_shared<const U3GateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -_theta, -_lambda, -_phi);
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat(2, 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> 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.);
@@ -679,8 +681,8 @@ class SwapGateImpl : public GateBase<Fp> {
     std::shared_ptr<const GateBase<Fp>> get_inverse() const override {
         return this->shared_from_this();
     }
-    internal::ComplexMatrix get_matrix() const override {
-        internal::ComplexMatrix mat = internal::ComplexMatrix::Identity(1 << 2, 1 << 2);
+    internal::ComplexMatrix<Fp> get_matrix() const override {
+        internal::ComplexMatrix<Fp> mat = internal::ComplexMatrix<Fp>::Identity(1 << 2, 1 << 2);
         mat << 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1;
         return mat;
     }
diff --git a/scaluq/gate/param_gate.hpp b/scaluq/gate/param_gate.hpp
index 894d8c77..fb3bbf14 100644
--- a/scaluq/gate/param_gate.hpp
+++ b/scaluq/gate/param_gate.hpp
@@ -116,7 +116,7 @@ class ParamGateBase {
     }
 
     [[nodiscard]] virtual std::shared_ptr<const ParamGateBase<Fp>> get_inverse() const = 0;
-    [[nodiscard]] virtual internal::ComplexMatrix get_matrix(double param) const = 0;
+    [[nodiscard]] virtual internal::ComplexMatrix<Fp> get_matrix(double param) const = 0;
 
     virtual void update_quantum_state(StateVector<Fp>& state_vector, double param) const = 0;
 
diff --git a/scaluq/gate/param_gate_pauli.hpp b/scaluq/gate/param_gate_pauli.hpp
index 6e9798c4..0742560f 100644
--- a/scaluq/gate/param_gate_pauli.hpp
+++ b/scaluq/gate/param_gate_pauli.hpp
@@ -28,14 +28,14 @@ class ParamPauliRotationGateImpl : public ParamGateBase<Fp> {
         return std::make_shared<const ParamPauliRotationGateImpl<Fp>>(
             this->_control_mask, _pauli, -this->_pcoef);
     }
-    internal::ComplexMatrix get_matrix(Fp param) const override {
+    internal::ComplexMatrix<Fp> get_matrix(Fp param) const override {
         Fp angle = this->_pcoef * param;
-        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) *
-                  internal::ComplexMatrix::Identity(mat.rows(), mat.cols()) +
-              imag_unit * (StdComplex)Kokkos::sin(-true_angle / 2) * mat;
+        Complex<Fp> true_angle = angle * this->_pauli.coef();
+        internal::ComplexMatrix<Fp> mat = this->_pauli.get_matrix_ignoring_coef();
+        StdComplex<Fp> imag_unit(0, 1);
+        mat = (StdComplex<Fp>)Kokkos::cos(-true_angle / 2) *
+                  internal::ComplexMatrix<Fp>::Identity(mat.rows(), mat.cols()) +
+              imag_unit * (StdComplex<Fp>)Kokkos::sin(-true_angle / 2) * mat;
         return mat;
     }
     void update_quantum_state(StateVector<Fp>& state_vector, Fp param) const override {
diff --git a/scaluq/gate/param_gate_probablistic.hpp b/scaluq/gate/param_gate_probablistic.hpp
index 4b5d5486..117e1c55 100644
--- a/scaluq/gate/param_gate_probablistic.hpp
+++ b/scaluq/gate/param_gate_probablistic.hpp
@@ -76,7 +76,7 @@ class ParamProbablisticGateImpl : public ParamGateBase<Fp> {
             });
         return std::make_shared<const ParamProbablisticGateImpl>(_distribution, inv_gate_list);
     }
-    internal::ComplexMatrix get_matrix(Fp) const override {
+    internal::ComplexMatrix<Fp> get_matrix(Fp) const override {
         throw std::runtime_error(
             "ParamProbablisticGateImpl::get_matrix(): This function must not be used in "
             "ParamProbablisticGateImpl.");
diff --git a/scaluq/gate/param_gate_standard.hpp b/scaluq/gate/param_gate_standard.hpp
index efb2e9bd..23991b0d 100644
--- a/scaluq/gate/param_gate_standard.hpp
+++ b/scaluq/gate/param_gate_standard.hpp
@@ -17,9 +17,9 @@ class ParamRXGateImpl : public ParamGateBase<Fp> {
         return std::make_shared<const ParamRXGateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -this->_pcoef);
     }
-    internal::ComplexMatrix get_matrix(Fp param) const override {
+    internal::ComplexMatrix<Fp> get_matrix(Fp param) const override {
         Fp angle = this->_pcoef * param;
-        internal::ComplexMatrix mat(2, 2);
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << std::cos(angle / 2), -1i * std::sin(angle / 2), -1i * std::sin(angle / 2),
             std::cos(angle / 2);
         return mat;
@@ -47,9 +47,9 @@ class ParamRYGateImpl : public ParamGateBase<Fp> {
         return std::make_shared<const ParamRYGateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -this->_pcoef);
     }
-    internal::ComplexMatrix get_matrix(Fp param) const override {
+    internal::ComplexMatrix<Fp> get_matrix(Fp param) const override {
         Fp angle = this->_pcoef * param;
-        internal::ComplexMatrix mat(2, 2);
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << std::cos(angle / 2), -std::sin(angle / 2), std::sin(angle / 2), std::cos(angle / 2);
         return mat;
     }
@@ -76,9 +76,9 @@ class ParamRZGateImpl : public ParamGateBase<Fp> {
         return std::make_shared<const ParamRZGateImpl<Fp>>(
             this->_target_mask, this->_control_mask, -this->_pcoef);
     }
-    internal::ComplexMatrix get_matrix(Fp param) const override {
+    internal::ComplexMatrix<Fp> get_matrix(Fp param) const override {
         Fp angle = param * this->_pcoef;
-        internal::ComplexMatrix mat(2, 2);
+        internal::ComplexMatrix<Fp> mat(2, 2);
         mat << std::exp(-0.5i * angle), 0, 0, std::exp(0.5i * angle);
         return mat;
     }
diff --git a/scaluq/gate/update_ops.hpp b/scaluq/gate/update_ops.hpp
index 41acc1b9..947804df 100644
--- a/scaluq/gate/update_ops.hpp
+++ b/scaluq/gate/update_ops.hpp
@@ -9,28 +9,28 @@ namespace scaluq {
 namespace internal {
 
 template <std::floating_point Fp>
-inline Matrix2x2 get_IBMQ_matrix(Fp _theta, Fp _phi, Fp _lambda) {
-    Complex exp_val1 = Kokkos::exp(Complex(0, _phi));
-    Complex exp_val2 = Kokkos::exp(Complex(0, _lambda));
-    Complex cos_val = Kokkos::cos(_theta / 2.);
-    Complex sin_val = Kokkos::sin(_theta / 2.);
+inline Matrix2x2<Fp> get_IBMQ_matrix(Fp _theta, Fp _phi, Fp _lambda) {
+    Complex<Fp> exp_val1 = Kokkos::exp(Complex<Fp>(0, _phi));
+    Complex<Fp> exp_val2 = Kokkos::exp(Complex<Fp>(0, _lambda));
+    Complex<Fp> cos_val = Kokkos::cos(_theta / 2.);
+    Complex<Fp> sin_val = Kokkos::sin(_theta / 2.);
     return {cos_val, -exp_val2 * sin_val, exp_val1 * sin_val, exp_val1 * exp_val2 * cos_val};
 }
 
 template <std::floating_point Fp>
 inline void one_target_dense_matrix_gate(std::uint64_t target_mask,
                                          std::uint64_t control_mask,
-                                         const Matrix2x2& matrix,
+                                         const Matrix2x2<Fp>& matrix,
                                          StateVector<Fp>& state) {
     Kokkos::parallel_for(
         state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) {
             std::uint64_t basis_0 =
                 insert_zero_at_mask_positions(it, control_mask | target_mask) | control_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[0][0] * val0 + matrix[0][1] * val1;
-            Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1;
+            Complex<Fp> val0 = state._raw[basis_0];
+            Complex<Fp> val1 = state._raw[basis_1];
+            Complex<Fp> res0 = matrix[0][0] * val0 + matrix[0][1] * val1;
+            Complex<Fp> res1 = matrix[1][0] * val0 + matrix[1][1] * val1;
             state._raw[basis_0] = res0;
             state._raw[basis_1] = res1;
         });
@@ -40,7 +40,7 @@ inline void one_target_dense_matrix_gate(std::uint64_t target_mask,
 template <std::floating_point Fp>
 inline void two_target_dense_matrix_gate(std::uint64_t target_mask,
                                          std::uint64_t control_mask,
-                                         const Matrix4x4& matrix,
+                                         const Matrix4x4<Fp>& matrix,
                                          StateVector<Fp>& state) {
     std::uint64_t lower_target_mask = -target_mask & target_mask;
     std::uint64_t upper_target_mask = target_mask ^ lower_target_mask;
@@ -52,18 +52,18 @@ inline void two_target_dense_matrix_gate(std::uint64_t target_mask,
             std::uint64_t basis_1 = basis_0 | lower_target_mask;
             std::uint64_t basis_2 = basis_0 | upper_target_mask;
             std::uint64_t basis_3 = basis_1 | target_mask;
-            Complex val0 = state._raw[basis_0];
-            Complex val1 = state._raw[basis_1];
-            Complex val2 = state._raw[basis_2];
-            Complex val3 = state._raw[basis_3];
-            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;
+            Complex<Fp> val0 = state._raw[basis_0];
+            Complex<Fp> val1 = state._raw[basis_1];
+            Complex<Fp> val2 = state._raw[basis_2];
+            Complex<Fp> val3 = state._raw[basis_3];
+            Complex<Fp> res0 = matrix[0][0] * val0 + matrix[0][1] * val1 + matrix[0][2] * val2 +
+                               matrix[0][3] * val3;
+            Complex<Fp> res1 = matrix[1][0] * val0 + matrix[1][1] * val1 + matrix[1][2] * val2 +
+                               matrix[1][3] * val3;
+            Complex<Fp> res2 = matrix[2][0] * val0 + matrix[2][1] * val1 + matrix[2][2] * val2 +
+                               matrix[2][3] * val3;
+            Complex<Fp> 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;
@@ -75,7 +75,7 @@ inline void two_target_dense_matrix_gate(std::uint64_t target_mask,
 template <std::floating_point Fp>
 inline void one_target_diagonal_matrix_gate(std::uint64_t target_mask,
                                             std::uint64_t control_mask,
-                                            const DiagonalMatrix2x2& diag,
+                                            const DiagonalMatrix2x2<Fp>& diag,
                                             StateVector<Fp>& state) {
     Kokkos::parallel_for(
         state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) {
@@ -95,7 +95,7 @@ inline void global_phase_gate(std::uint64_t target_mask,
                               std::uint64_t control_mask,
                               Fp angle,
                               StateVector<Fp>& state) {
-    Complex coef = Kokkos::polar(1., angle);
+    Complex<Fp> coef = Kokkos::polar(1., angle);
     Kokkos::parallel_for(
         state.dim() >> std::popcount(control_mask), KOKKOS_LAMBDA(std::uint64_t i) {
             state._raw[insert_zero_at_mask_positions(i, control_mask) | control_mask] *= coef;
@@ -120,8 +120,8 @@ inline void y_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateV
         state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) {
             std::uint64_t i =
                 insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask;
-            state._raw[i] *= Complex(0, 1);
-            state._raw[i | target_mask] *= Complex(0, -1);
+            state._raw[i] *= Complex<Fp>(0, 1);
+            state._raw[i | target_mask] *= Complex<Fp>(0, -1);
             Kokkos::Experimental::swap(state._raw[i], state._raw[i | target_mask]);
         });
     Kokkos::fence();
@@ -133,20 +133,20 @@ inline void z_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateV
         state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) {
             std::uint64_t i =
                 insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask;
-            state._raw[i | target_mask] *= Complex(-1, 0);
+            state._raw[i | target_mask] *= Complex<Fp>(-1, 0);
         });
     Kokkos::fence();
 }
 
 template <std::floating_point Fp>
 inline void h_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, HADAMARD_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, HADAMARD_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
 inline void one_target_phase_gate(std::uint64_t target_mask,
                                   std::uint64_t control_mask,
-                                  Complex phase,
+                                  Complex<Fp> phase,
                                   StateVector<Fp>& state) {
     Kokkos::parallel_for(
         state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) {
@@ -159,20 +159,20 @@ inline void one_target_phase_gate(std::uint64_t target_mask,
 
 template <std::floating_point Fp>
 inline void s_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector<Fp>& state) {
-    one_target_phase_gate(target_mask, control_mask, Complex(0, 1), state);
+    one_target_phase_gate(target_mask, control_mask, Complex<Fp>(0, 1), state);
 }
 
 template <std::floating_point Fp>
 inline void sdag_gate(std::uint64_t target_mask,
                       std::uint64_t control_mask,
                       StateVector<Fp>& state) {
-    one_target_phase_gate(target_mask, control_mask, Complex(0, -1), state);
+    one_target_phase_gate(target_mask, control_mask, Complex<Fp>(0, -1), state);
 }
 
 template <std::floating_point Fp>
 inline void t_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector<Fp>& state) {
     one_target_phase_gate(
-        target_mask, control_mask, Complex(INVERSE_SQRT2(), INVERSE_SQRT2()), state);
+        target_mask, control_mask, Complex<Fp>(INVERSE_SQRT2(), INVERSE_SQRT2()), state);
 }
 
 template <std::floating_point Fp>
@@ -180,45 +180,45 @@ inline void tdag_gate(std::uint64_t target_mask,
                       std::uint64_t control_mask,
                       StateVector<Fp>& state) {
     one_target_phase_gate(
-        target_mask, control_mask, Complex(INVERSE_SQRT2(), -INVERSE_SQRT2()), state);
+        target_mask, control_mask, Complex<Fp>(INVERSE_SQRT2(), -INVERSE_SQRT2()), state);
 }
 
 template <std::floating_point Fp>
 inline void sqrtx_gate(std::uint64_t target_mask,
                        std::uint64_t control_mask,
                        StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_GATE_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_GATE_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
 inline void sqrtxdag_gate(std::uint64_t target_mask,
                           std::uint64_t control_mask,
                           StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_DAG_GATE_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_DAG_GATE_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
 inline void sqrty_gate(std::uint64_t target_mask,
                        std::uint64_t control_mask,
                        StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_GATE_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_GATE_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
 inline void sqrtydag_gate(std::uint64_t target_mask,
                           std::uint64_t control_mask,
                           StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_DAG_GATE_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_DAG_GATE_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
 inline void p0_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, PROJ_0_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, PROJ_0_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
 inline void p1_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector<Fp>& state) {
-    one_target_dense_matrix_gate(target_mask, control_mask, PROJ_1_MATRIX(), state);
+    one_target_dense_matrix_gate(target_mask, control_mask, PROJ_1_MATRIX<Fp>(), state);
 }
 
 template <std::floating_point Fp>
@@ -228,7 +228,7 @@ inline void rx_gate(std::uint64_t target_mask,
                     StateVector<Fp>& state) {
     const Fp cosval = std::cos(angle / 2.);
     const Fp sinval = std::sin(angle / 2.);
-    Matrix2x2 matrix = {cosval, Complex(0, -sinval), Complex(0, -sinval), cosval};
+    Matrix2x2<Fp> matrix = {cosval, Complex<Fp>(0, -sinval), Complex<Fp>(0, -sinval), cosval};
     one_target_dense_matrix_gate(target_mask, control_mask, matrix, state);
 }
 
@@ -239,7 +239,7 @@ inline void ry_gate(std::uint64_t target_mask,
                     StateVector<Fp>& state) {
     const Fp cosval = std::cos(angle / 2.);
     const Fp sinval = std::sin(angle / 2.);
-    Matrix2x2 matrix = {cosval, -sinval, sinval, cosval};
+    Matrix2x2<Fp> matrix = {cosval, -sinval, sinval, cosval};
     one_target_dense_matrix_gate(target_mask, control_mask, matrix, state);
 }
 
@@ -250,7 +250,7 @@ inline void rz_gate(std::uint64_t target_mask,
                     StateVector<Fp>& state) {
     const Fp cosval = std::cos(angle / 2.);
     const Fp sinval = std::sin(angle / 2.);
-    DiagonalMatrix2x2 diag = {Complex(cosval, -sinval), Complex(cosval, sinval)};
+    DiagonalMatrix2x2<double> diag = {Complex<Fp>(cosval, -sinval), Complex<Fp>(cosval, sinval)};
     one_target_diagonal_matrix_gate(target_mask, control_mask, diag, state);
 }
 
@@ -259,7 +259,7 @@ inline void u1_gate(std::uint64_t target_mask,
                     std::uint64_t control_mask,
                     Fp lambda,
                     StateVector<Fp>& state) {
-    Complex exp_val = Kokkos::exp(Complex(0, lambda));
+    Complex<Fp> exp_val = Kokkos::exp(Complex<Fp>(0, lambda));
     Kokkos::parallel_for(
         state.dim() >> (std::popcount(target_mask | control_mask)),
         KOKKOS_LAMBDA(std::uint64_t it) {
diff --git a/scaluq/operator/apply_pauli.hpp b/scaluq/operator/apply_pauli.hpp
index 38c097aa..13d727f5 100644
--- a/scaluq/operator/apply_pauli.hpp
+++ b/scaluq/operator/apply_pauli.hpp
@@ -9,7 +9,7 @@ template <std::floating_point Fp>
 void apply_pauli(std::uint64_t control_mask,
                  std::uint64_t bit_flip_mask,
                  std::uint64_t phase_flip_mask,
-                 Complex coef,
+                 Complex<Fp> coef,
                  StateVector<Fp>& state_vector) {
     if (bit_flip_mask == 0) {
         Kokkos::parallel_for(
@@ -27,14 +27,14 @@ void apply_pauli(std::uint64_t control_mask,
     }
     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];
+    Complex<Fp> global_phase = PHASE_M90ROT<Fp>()[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;
+            Complex<Fp> tmp1 = state_vector._raw[basis_0] * global_phase;
+            Complex<Fp> 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;
@@ -47,16 +47,16 @@ template <std::floating_point Fp>
 void apply_pauli_rotation(std::uint64_t control_mask,
                           std::uint64_t bit_flip_mask,
                           std::uint64_t phase_flip_mask,
-                          Complex coef,
+                          Complex<Fp> coef,
                           Fp angle,
                           StateVector<Fp>& 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);
+    Complex<Fp> true_angle = angle * coef;
+    const Complex<Fp> cosval = Kokkos::cos(-true_angle / 2);
+    const Complex<Fp> 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;
+        const Complex<Fp> cval_min = cosval - Complex<Fp>(0, 1) * sinval;
+        const Complex<Fp> cval_pls = cosval + Complex<Fp>(0, 1) * sinval;
         Kokkos::parallel_for(
             state_vector.dim() >> std::popcount(control_mask), KOKKOS_LAMBDA(std::uint64_t i) {
                 std::uint64_t state_idx =
@@ -83,18 +83,18 @@ void apply_pauli_rotation(std::uint64_t control_mask,
                 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];
+                Complex<Fp> cval_0 = state_vector._raw[basis_0];
+                Complex<Fp> 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];
+                    Complex<Fp>(0, 1) * sinval * cval_1 *
+                        PHASE_M90ROT<Fp>()[(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];
+                    Complex<Fp>(0, 1) * sinval * cval_0 *
+                        PHASE_M90ROT<Fp>()[(global_phase_90_rot_count + bit_parity_1 * 2) % 4];
             });
         Kokkos::fence();
     }
diff --git a/scaluq/operator/operator.hpp b/scaluq/operator/operator.hpp
index 02ecb19d..c68401a8 100644
--- a/scaluq/operator/operator.hpp
+++ b/scaluq/operator/operator.hpp
@@ -52,13 +52,13 @@ class Operator {
                 target_qubit_list[qubit_idx] = qubit_idx;
                 pauli_id_list[qubit_idx] = random.int32() & 0b11;
             }
-            Complex coef = random.uniform() * 2. - 1.;
+            Complex<Fp> coef = random.uniform() * 2. - 1.;
             this->add_operator(PauliOperator(target_qubit_list, pauli_id_list, coef));
         }
     }
 
     void optimize() {
-        std::map<std::tuple<std::uint64_t, std::uint64_t>, Complex> pauli_and_coef;
+        std::map<std::tuple<std::uint64_t, std::uint64_t>, Complex<Fp>> pauli_and_coef;
         for (const auto& pauli : _terms) {
             pauli_and_coef[pauli.get_XZ_mask_representation()] += pauli.coef();
         }
@@ -91,7 +91,7 @@ class Operator {
         state_vector = res;
     }
 
-    [[nodiscard]] Complex get_expectation_value(const StateVector<Fp>& state_vector) const {
+    [[nodiscard]] Complex<Fp> get_expectation_value(const StateVector<Fp>& state_vector) const {
         if (_n_qubits > state_vector.n_qubits()) {
             throw std::runtime_error(
                 "Operator::get_expectation_value: n_qubits of state_vector is too small");
@@ -103,7 +103,7 @@ class Operator {
             terms_view(_terms.data(), nterms);
         Kokkos::View<std::uint64_t*, Kokkos::HostSpace> bmasks_host("bmasks_host", nterms);
         Kokkos::View<std::uint64_t*, Kokkos::HostSpace> pmasks_host("pmasks_host", nterms);
-        Kokkos::View<Complex*, Kokkos::HostSpace> coefs_host("coefs_host", nterms);
+        Kokkos::View<Complex<Fp>*, Kokkos::HostSpace> coefs_host("coefs_host", nterms);
         Kokkos::Experimental::transform(
             Kokkos::DefaultHostExecutionSpace(),
             terms_view,
@@ -121,18 +121,18 @@ class Operator {
             [](const PauliOperator<Fp>& pauli) { return pauli._ptr->_coef; });
         Kokkos::View<std::uint64_t*> bmasks("bmasks", nterms);
         Kokkos::View<std::uint64_t*> pmasks("pmasks", nterms);
-        Kokkos::View<Complex*> coefs("coefs", nterms);
+        Kokkos::View<Complex<Fp>*> coefs("coefs", nterms);
         Kokkos::deep_copy(bmasks, bmasks_host);
         Kokkos::deep_copy(pmasks, pmasks_host);
         Kokkos::deep_copy(coefs, coefs_host);
         std::uint64_t dim = state_vector.dim();
-        Complex res;
+        Complex<Fp> res;
         Kokkos::parallel_reduce(
             Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nterms, dim >> 1}),
-            KOKKOS_LAMBDA(std::uint64_t term_id, std::uint64_t state_idx, Complex & res_lcl) {
+            KOKKOS_LAMBDA(std::uint64_t term_id, std::uint64_t state_idx, Complex<Fp> & res_lcl) {
                 std::uint64_t bit_flip_mask = bmasks[term_id];
                 std::uint64_t phase_flip_mask = pmasks[term_id];
-                Complex coef = coefs[term_id];
+                Complex<Fp> coef = coefs[term_id];
                 if (bit_flip_mask == 0) {
                     std::uint64_t state_idx1 = state_idx << 1;
                     Fp tmp1 = (Kokkos::conj(state_vector._raw[state_idx1]) *
@@ -150,7 +150,8 @@ class Operator {
                         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 = internal::PHASE_90ROT()[global_phase_90rot_count % 4];
+                    Complex<Fp> global_phase =
+                        internal::PHASE_90ROT<Fp>()[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;
                     Fp tmp =
@@ -165,8 +166,8 @@ class Operator {
         return res;
     }
 
-    [[nodiscard]] Complex get_transition_amplitude(const StateVector<Fp>& state_vector_bra,
-                                                   const StateVector<Fp>& state_vector_ket) const {
+    [[nodiscard]] Complex<Fp> get_transition_amplitude(
+        const StateVector<Fp>& state_vector_bra, const StateVector<Fp>& state_vector_ket) const {
         if (state_vector_bra.n_qubits() != state_vector_ket.n_qubits()) {
             throw std::runtime_error(
                 "Operator::get_transition_amplitude: n_qubits of state_vector_bra and "
@@ -180,7 +181,7 @@ class Operator {
         std::uint64_t nterms = _terms.size();
         std::vector<std::uint64_t> bmasks_vector(nterms);
         std::vector<std::uint64_t> pmasks_vector(nterms);
-        std::vector<Complex> coefs_vector(nterms);
+        std::vector<Complex<Fp>> coefs_vector(nterms);
         std::transform(_terms.begin(),
                        _terms.end(),
                        bmasks_vector.begin(),
@@ -197,23 +198,24 @@ class Operator {
             internal::convert_host_vector_to_device_view(bmasks_vector);
         Kokkos::View<std::uint64_t*> pmasks =
             internal::convert_host_vector_to_device_view(pmasks_vector);
-        Kokkos::View<Complex*> coefs = internal::convert_host_vector_to_device_view(coefs_vector);
+        Kokkos::View<Complex<Fp>*> coefs =
+            internal::convert_host_vector_to_device_view(coefs_vector);
         std::uint64_t dim = state_vector_bra.dim();
-        Complex res;
+        Complex<Fp> res;
         Kokkos::parallel_reduce(
             Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nterms, dim >> 1}),
-            KOKKOS_LAMBDA(std::uint64_t term_id, std::uint64_t state_idx, Complex & res_lcl) {
+            KOKKOS_LAMBDA(std::uint64_t term_id, std::uint64_t state_idx, Complex<Fp> & res_lcl) {
                 std::uint64_t bit_flip_mask = bmasks[term_id];
                 std::uint64_t phase_flip_mask = pmasks[term_id];
-                Complex coef = coefs[term_id];
+                Complex<Fp> coef = coefs[term_id];
                 if (bit_flip_mask == 0) {
                     std::uint64_t state_idx1 = state_idx << 1;
-                    Complex tmp1 = (Kokkos::conj(state_vector_bra._raw[state_idx1]) *
-                                    state_vector_ket._raw[state_idx1]);
+                    Complex<Fp> tmp1 = (Kokkos::conj(state_vector_bra._raw[state_idx1]) *
+                                        state_vector_ket._raw[state_idx1]);
                     if (Kokkos::popcount(state_idx1 & phase_flip_mask) & 1) tmp1 = -tmp1;
                     std::uint64_t state_idx2 = state_idx1 | 1;
-                    Complex tmp2 = (Kokkos::conj(state_vector_bra._raw[state_idx2]) *
-                                    state_vector_ket._raw[state_idx2]);
+                    Complex<Fp> tmp2 = (Kokkos::conj(state_vector_bra._raw[state_idx2]) *
+                                        state_vector_ket._raw[state_idx2]);
                     if (Kokkos::popcount(state_idx2 & phase_flip_mask) & 1) tmp2 = -tmp2;
                     res_lcl += coef * (tmp1 + tmp2);
                 } else {
@@ -221,14 +223,15 @@ class Operator {
                         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 = internal::PHASE_90ROT()[global_phase_90rot_count % 4];
+                    Complex<Fp> global_phase =
+                        internal::PHASE_90ROT<Fp>()[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]) *
-                                   state_vector_ket._raw[basis_0] * global_phase;
+                    Complex<Fp> tmp1 = Kokkos::conj(state_vector_bra._raw[basis_1]) *
+                                       state_vector_ket._raw[basis_0] * global_phase;
                     if (Kokkos::popcount(basis_0 & phase_flip_mask) & 1) tmp1 = -tmp1;
-                    Complex tmp2 = Kokkos::conj(state_vector_bra._raw[basis_0]) *
-                                   state_vector_ket._raw[basis_1] * global_phase;
+                    Complex<Fp> tmp2 = Kokkos::conj(state_vector_bra._raw[basis_0]) *
+                                       state_vector_ket._raw[basis_1] * global_phase;
                     if (Kokkos::popcount(basis_1 & phase_flip_mask) & 1) tmp2 = -tmp2;
                     res_lcl += coef * (tmp1 + tmp2);
                 }
@@ -239,20 +242,19 @@ class Operator {
     }
 
     // not implemented yet
-    [[nodiscard]] Complex solve_gound_state_eigenvalue_by_arnoldi_method(
-        const StateVector<Fp>& state, std::uint64_t iter_count, Complex mu = 0.) const;
+    [[nodiscard]] Complex<Fp> solve_gound_state_eigenvalue_by_arnoldi_method(
+        const StateVector<Fp>& state, std::uint64_t iter_count, Complex<Fp> mu = 0.) const;
     // not implemented yet
-    [[nodiscard]] Complex solve_gound_state_eigenvalue_by_power_method(const StateVector<Fp>& state,
-                                                                       std::uint64_t iter_count,
-                                                                       Complex mu = 0.) const;
+    [[nodiscard]] Complex<Fp> solve_gound_state_eigenvalue_by_power_method(
+        const StateVector<Fp>& state, std::uint64_t iter_count, Complex<Fp> mu = 0.) const;
 
-    Operator& operator*=(Complex coef) {
+    Operator& operator*=(Complex<Fp> coef) {
         for (auto& pauli : _terms) {
             pauli = pauli * coef;
         }
         return *this;
     }
-    Operator operator*(Complex coef) const { return Operator(*this) *= coef; }
+    Operator operator*(Complex<Fp> coef) const { return Operator(*this) *= coef; }
     inline Operator operator+() const { return *this; }
     Operator operator-() const { return *this * -1; }
     Operator& operator+=(const Operator& target) {
@@ -336,8 +338,8 @@ void bind_operator_operator_hpp(nb::module_& m) {
         .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 *= Complex<Fp>())
+        .def(nb::self * Complex<Fp>())
         .def(+nb::self)
         .def(-nb::self)
         .def(nb::self += nb::self)
diff --git a/scaluq/operator/pauli_operator.hpp b/scaluq/operator/pauli_operator.hpp
index 493643c1..1ab643b0 100644
--- a/scaluq/operator/pauli_operator.hpp
+++ b/scaluq/operator/pauli_operator.hpp
@@ -22,12 +22,13 @@ class PauliOperator {
         friend class PauliOperator<Fp>;
         friend class Operator<Fp>;
         std::vector<std::uint64_t> _target_qubit_list, _pauli_id_list;
-        Complex _coef;
+        Complex<Fp> _coef;
         std::uint64_t _bit_flip_mask, _phase_flip_mask;
 
     public:
-        explicit Data(Complex coef = 1.) : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) {}
-        Data(std::string_view pauli_string, Complex coef = 1.)
+        explicit Data(Complex<Fp> coef = 1.)
+            : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) {}
+        Data(std::string_view pauli_string, Complex<Fp> coef = 1.)
             : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) {
             auto ss = std::stringstream(std::string(pauli_string));
             while (1) {
@@ -54,7 +55,7 @@ class PauliOperator {
 
         Data(const std::vector<std::uint64_t>& target_qubit_list,
              const std::vector<std::uint64_t>& pauli_id_list,
-             Complex coef = 1.)
+             Complex<Fp> coef = 1.)
             : _coef(coef), _bit_flip_mask(0), _phase_flip_mask(0) {
             if (target_qubit_list.size() != pauli_id_list.size()) {
                 throw std::runtime_error(
@@ -67,14 +68,14 @@ class PauliOperator {
             }
         }
 
-        Data(const std::vector<std::uint64_t>& pauli_id_par_qubit, Complex coef = 1.)
+        Data(const std::vector<std::uint64_t>& pauli_id_par_qubit, Complex<Fp> coef = 1.)
             : _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]);
             }
         }
 
-        Data(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex coef)
+        Data(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex<Fp> 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++) {
@@ -120,8 +121,8 @@ class PauliOperator {
             }
         }
 
-        Complex coef() const { return _coef; }
-        void set_coef(Complex c) { _coef = c; }
+        Complex<Fp> coef() const { return _coef; }
+        void set_coef(Complex<Fp> c) { _coef = c; }
         const std::vector<std::uint64_t>& target_qubit_list() const { return _target_qubit_list; }
         const std::vector<std::uint64_t>& pauli_id_list() const { return _pauli_id_list; }
         std::tuple<std::uint64_t, std::uint64_t> get_XZ_mask_representation() const {
@@ -135,20 +136,20 @@ class PauliOperator {
 public:
     enum PauliID : std::uint64_t { I, X, Y, Z };
 
-    explicit PauliOperator(Complex coef = 1.) : _ptr(std::make_shared<const Data>(coef)) {}
+    explicit PauliOperator(Complex<Fp> coef = 1.) : _ptr(std::make_shared<const Data>(coef)) {}
     explicit PauliOperator(Data data) : _ptr(std::make_shared<const Data>(data)) {}
-    PauliOperator(std::string_view pauli_string, Complex coef = 1.)
+    PauliOperator(std::string_view pauli_string, Complex<Fp> coef = 1.)
         : _ptr(std::make_shared<const Data>(pauli_string, coef)) {}
     PauliOperator(const std::vector<std::uint64_t>& target_qubit_list,
                   const std::vector<std::uint64_t>& pauli_id_list,
-                  Complex coef = 1.)
+                  Complex<Fp> coef = 1.)
         : _ptr(std::make_shared<const Data>(target_qubit_list, pauli_id_list, coef)) {}
-    PauliOperator(const std::vector<std::uint64_t>& pauli_id_par_qubit, Complex coef = 1.)
+    PauliOperator(const std::vector<std::uint64_t>& pauli_id_par_qubit, Complex<Fp> coef = 1.)
         : _ptr(std::make_shared<const Data>(pauli_id_par_qubit, coef)) {}
-    PauliOperator(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex coef = 1.)
+    PauliOperator(std::uint64_t bit_flip_mask, std::uint64_t phase_flip_mask, Complex<Fp> coef = 1.)
         : _ptr(std::make_shared<const Data>(bit_flip_mask, phase_flip_mask, coef)) {}
 
-    [[nodiscard]] inline Complex coef() const { return _ptr->coef(); }
+    [[nodiscard]] inline Complex<Fp> coef() const { return _ptr->coef(); }
     [[nodiscard]] inline const std::vector<std::uint64_t>& target_qubit_list() const {
         return _ptr->target_qubit_list();
     }
@@ -192,7 +193,7 @@ class PauliOperator {
             0ULL, _ptr->_bit_flip_mask, _ptr->_phase_flip_mask, _ptr->_coef, state_vector);
     }
 
-    [[nodiscard]] Complex get_expectation_value(const StateVector<Fp>& state_vector) const {
+    [[nodiscard]] Complex<Fp> get_expectation_value(const StateVector<Fp>& 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 "
@@ -217,7 +218,7 @@ class PauliOperator {
         }
         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 = internal::PHASE_90ROT()[global_phase_90rot_count % 4];
+        Complex<Fp> global_phase = internal::PHASE_90ROT<Fp>()[global_phase_90rot_count % 4];
         Fp res;
         Kokkos::parallel_reduce(
             state_vector.dim() >> 1,
@@ -232,8 +233,8 @@ class PauliOperator {
             res);
         return _ptr->_coef * res;
     }
-    [[nodiscard]] Complex get_transition_amplitude(const StateVector<Fp>& state_vector_bra,
-                                                   const StateVector<Fp>& state_vector_ket) const {
+    [[nodiscard]] Complex<Fp> get_transition_amplitude(
+        const StateVector<Fp>& state_vector_bra, const StateVector<Fp>& state_vector_ket) const {
         if (state_vector_bra.n_qubits() != state_vector_ket.n_qubits()) {
             throw std::runtime_error(
                 "state_vector_bra must have same n_qubits to state_vector_ket.");
@@ -247,12 +248,12 @@ class PauliOperator {
         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;
+            Complex<Fp> res;
             Kokkos::parallel_reduce(
                 state_vector_bra.dim(),
-                KOKKOS_LAMBDA(std::uint64_t state_idx, Complex & sum) {
-                    Complex tmp = Kokkos::conj(state_vector_bra._raw[state_idx]) *
-                                  state_vector_ket._raw[state_idx];
+                KOKKOS_LAMBDA(std::uint64_t state_idx, Complex<Fp> & sum) {
+                    Complex<Fp> tmp = Kokkos::conj(state_vector_bra._raw[state_idx]) *
+                                      state_vector_ket._raw[state_idx];
                     if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) tmp = -tmp;
                     sum += tmp;
                 },
@@ -262,18 +263,18 @@ class PauliOperator {
         }
         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 = internal::PHASE_90ROT()[global_phase_90rot_count % 4];
-        Complex res;
+        Complex<Fp> global_phase = internal::PHASE_90ROT<Fp>()[global_phase_90rot_count % 4];
+        Complex<Fp> res;
         Kokkos::parallel_reduce(
             state_vector_bra.dim() >> 1,
-            KOKKOS_LAMBDA(std::uint64_t state_idx, Complex & sum) {
+            KOKKOS_LAMBDA(std::uint64_t state_idx, Complex<Fp> & sum) {
                 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]) *
-                               state_vector_ket._raw[basis_0] * global_phase;
+                Complex<Fp> tmp1 = Kokkos::conj(state_vector_bra._raw[basis_1]) *
+                                   state_vector_ket._raw[basis_0] * global_phase;
                 if (Kokkos::popcount(basis_0 & phase_flip_mask) & 1) tmp1 = -tmp1;
-                Complex tmp2 = Kokkos::conj(state_vector_bra._raw[basis_0]) *
-                               state_vector_ket._raw[basis_1] * global_phase;
+                Complex<Fp> tmp2 = Kokkos::conj(state_vector_bra._raw[basis_0]) *
+                                   state_vector_ket._raw[basis_1] * global_phase;
                 if (Kokkos::popcount(basis_1 & phase_flip_mask) & 1) tmp2 = -tmp2;
                 sum += tmp1 + tmp2;
             },
@@ -282,11 +283,11 @@ class PauliOperator {
         return _ptr->_coef * res;
     }
 
-    [[nodiscard]] internal::ComplexMatrix get_matrix() const {
-        return get_matrix_ignoring_coef() * StdComplex(_ptr->_coef);
+    [[nodiscard]] internal::ComplexMatrix<Fp> get_matrix() const {
+        return get_matrix_ignoring_coef() * StdComplex<Fp>(_ptr->_coef);
     }
 
-    [[nodiscard]] internal::ComplexMatrix get_matrix_ignoring_coef() const {
+    [[nodiscard]] internal::ComplexMatrix<Fp> get_matrix_ignoring_coef() const {
         std::uint64_t flip_mask, phase_mask, rot90_count;
         Kokkos::parallel_reduce(
             Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(0, _ptr->_pauli_id_list.size()),
@@ -308,11 +309,11 @@ class PauliOperator {
             flip_mask,
             phase_mask,
             rot90_count);
-        std::vector<StdComplex> rot = {1, -1.i, -1, 1.i};
+        std::vector<StdComplex<Fp>> rot = {1, -1.i, -1, 1.i};
         std::uint64_t matrix_dim = 1ULL << _ptr->_pauli_id_list.size();
-        internal::ComplexMatrix mat = internal::ComplexMatrix::Zero(matrix_dim, matrix_dim);
+        internal::ComplexMatrix<Fp> mat = internal::ComplexMatrix<Fp>::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);
+            const StdComplex<Fp> sign = 1. - 2. * (Kokkos::popcount(index & phase_mask) % 2);
             mat(index, index ^ flip_mask) = rot[rot90_count % 4] * sign;
         }
         return mat;
@@ -337,9 +338,9 @@ class PauliOperator {
         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]);
+            _ptr->_coef * target._ptr->_coef * internal::PHASE_90ROT<Fp>()[extra_90rot_cnt]);
     }
-    [[nodiscard]] inline PauliOperator operator*(Complex target) const {
+    [[nodiscard]] inline PauliOperator operator*(Complex<Fp> target) const {
         return PauliOperator(_ptr->_target_qubit_list, _ptr->_pauli_id_list, _ptr->_coef * target);
     }
 };
@@ -356,23 +357,23 @@ void bind_operator_pauli_operator_hpp(nb::module_& m) {
 
     nb::class_<PauliOperator::Data>(
         m, "PauliOperatorData", "Internal data structure for PauliOperator.")
-        .def(nb::init<Complex>(), "coef"_a = 1., "Initialize data with coefficient.")
-        .def(nb::init<std::string_view, Complex>(),
+        .def(nb::init<Complex<Fp>>(), "coef"_a = 1., "Initialize data with coefficient.")
+        .def(nb::init<std::string_view, Complex<Fp>>(),
              "pauli_string"_a,
              "coef"_a = 1.,
              "Initialize data with pauli string.")
         .def(nb::init<const std::vector<std::uint64_t>&,
                       const std::vector<std::uint64_t>&,
-                      Complex>(),
+                      Complex<Fp>>(),
              "target_qubit_list"_a,
              "pauli_id_list"_a,
              "coef"_a = 1.,
              "Initialize data with target qubits and pauli ids.")
-        .def(nb::init<const std::vector<std::uint64_t>&, Complex>(),
+        .def(nb::init<const std::vector<std::uint64_t>&, Complex<Fp>>(),
              "pauli_id_par_qubit"_a,
              "coef"_a = 1.,
              "Initialize data with pauli ids per qubit.")
-        .def(nb::init<std::uint64_t, std::uint64_t, Complex>(),
+        .def(nb::init<std::uint64_t, std::uint64_t, Complex<Fp>>(),
              "bit_flip_mask"_a,
              "phase_flip_mask"_a,
              "coef"_a = 1.,
@@ -402,27 +403,29 @@ void bind_operator_pauli_operator_hpp(nb::module_& m) {
         m,
         "PauliOperator",
         "Pauli operator as coef and tensor product of single pauli for each qubit.")
-        .def(nb::init<Complex>(), "coef"_a = 1., "Initialize operator which just multiplying coef.")
+        .def(nb::init<Complex<Fp>>(),
+             "coef"_a = 1.,
+             "Initialize operator which just multiplying coef.")
         .def(nb::init<const std::vector<std::uint64_t>&,
                       const std::vector<std::uint64_t>&,
-                      Complex>(),
+                      Complex<Fp>>(),
              "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<std::string_view, Complex>(),
+        .def(nb::init<std::string_view, Complex<Fp>>(),
              "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<const std::vector<std::uint64_t>&, Complex>(),
+        .def(nb::init<const std::vector<std::uint64_t>&, Complex<Fp>>(),
              "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<std::uint64_t, std::uint64_t, Complex>(),
+        .def(nb::init<std::uint64_t, std::uint64_t, Complex<Fp>>(),
              "bit_flip_mask"_a,
              "phase_flip_mask"_a,
              "coef"_a = 1.,
@@ -460,7 +463,7 @@ void bind_operator_pauli_operator_hpp(nb::module_& m) {
              &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());
+        .def(nb::self * Complex<Fp>());
 }
 }  // namespace internal
 #endif
diff --git a/scaluq/types.hpp b/scaluq/types.hpp
index 753deb02..62b3446f 100644
--- a/scaluq/types.hpp
+++ b/scaluq/types.hpp
@@ -17,20 +17,29 @@ 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<double>;
-using Complex = Kokkos::complex<double>;
+template <std::floating_point Fp>
+using StdComplex = std::complex<Fp>;
+template <std::floating_point Fp>
+using Complex = Kokkos::complex<Fp>;
 using namespace std::complex_literals;
 
 namespace internal {
 template <typename DummyType>
 constexpr bool lazy_false_v = false;  // Used for lazy evaluation in static_assert.
 
-using ComplexMatrix = Eigen::Matrix<StdComplex, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-using SparseComplexMatrix = Eigen::SparseMatrix<StdComplex>;
+template <std::floating_point Fp>
+using ComplexMatrix =
+    Eigen::Matrix<StdComplex<Fp>, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
 
-using Matrix2x2 = Kokkos::Array<Kokkos::Array<Complex, 2>, 2>;
-using Matrix4x4 = Kokkos::Array<Kokkos::Array<Complex, 4>, 4>;
-using DiagonalMatrix2x2 = Kokkos::Array<Complex, 2>;
+template <std::floating_point Fp>
+using SparseComplexMatrix = Eigen::SparseMatrix<StdComplex<Fp>>;
+
+template <std::floating_point Fp>
+using Matrix2x2 = Kokkos::Array<Kokkos::Array<Complex<Fp>, 2>, 2>;
+template <std::floating_point Fp>
+using Matrix4x4 = Kokkos::Array<Kokkos::Array<Complex<Fp>, 4>, 4>;
+template <std::floating_point Fp>
+using DiagonalMatrix2x2 = Kokkos::Array<Complex<Fp>, 2>;
 }  // namespace internal
 
 #ifdef SCALUQ_USE_NANOBIND
diff --git a/scaluq/util/utility.hpp b/scaluq/util/utility.hpp
index c35a3bf1..3c7d3bd0 100644
--- a/scaluq/util/utility.hpp
+++ b/scaluq/util/utility.hpp
@@ -67,17 +67,19 @@ inline std::vector<std::uint64_t> mask_to_vector(std::uint64_t mask) {
     return indices;
 }
 
-KOKKOS_INLINE_FUNCTION Matrix2x2 matrix_multiply(const Matrix2x2& matrix1,
-                                                 const Matrix2x2& matrix2) {
+template <std::floating_point Fp>
+KOKKOS_INLINE_FUNCTION Matrix2x2<Fp> matrix_multiply(const Matrix2x2<Fp>& matrix1,
+                                                     const Matrix2x2<Fp>& 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 internal::ComplexMatrix kronecker_product(const internal::ComplexMatrix& lhs,
-                                                 const internal::ComplexMatrix& rhs) {
-    internal::ComplexMatrix result(lhs.rows() * rhs.rows(), lhs.cols() * rhs.cols());
+template <std::floating_point Fp>
+inline internal::ComplexMatrix<Fp> kronecker_product(const internal::ComplexMatrix<Fp>& lhs,
+                                                     const internal::ComplexMatrix<Fp>& rhs) {
+    internal::ComplexMatrix<Fp> 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;
@@ -86,9 +88,11 @@ inline internal::ComplexMatrix kronecker_product(const internal::ComplexMatrix&
     return result;
 }
 
-inline internal::ComplexMatrix get_expanded_matrix(const internal::ComplexMatrix& from_matrix,
-                                                   const std::vector<std::uint64_t>& from_targets,
-                                                   std::vector<std::uint64_t>& to_targets) {
+template <std::floating_point Fp>
+inline internal::ComplexMatrix<Fp> get_expanded_matrix(
+    const internal::ComplexMatrix<Fp>& from_matrix,
+    const std::vector<std::uint64_t>& from_targets,
+    std::vector<std::uint64_t>& to_targets) {
     std::vector<std::uint64_t> 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();
@@ -106,8 +110,8 @@ inline internal::ComplexMatrix get_expanded_matrix(const internal::ComplexMatrix
     for (std::uint64_t i : std::views::iota(0ULL, 1ULL << to_targets.size())) {
         if ((i & targets_idx_mask) == 0) outer_indices.push_back(i);
     }
-    internal::ComplexMatrix to_matrix =
-        internal::ComplexMatrix::Zero(1ULL << to_targets.size(), 1ULL << to_targets.size());
+    internal::ComplexMatrix<Fp> to_matrix =
+        internal::ComplexMatrix<Fp>::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 90855d31..77ac4af8 100644
--- a/tests/circuit/circuit_test.cpp
+++ b/tests/circuit/circuit_test.cpp
@@ -147,7 +147,7 @@ TEST(CircuitTest, CircuitBasic) {
     PauliOperator pauli = PauliOperator("I 0 X 1 Y 2 Z 3");
     circuit.add_gate(multi_Pauli(pauli));
 
-    internal::ComplexMatrix mat_x(2, 2);
+    internal::ComplexMatrix<double> mat_x(2, 2);
     target = random.int32() % n;
     mat_x << 0, 1, 1, 0;
     circuit.add_gate(dense_matrix(target, mat_x));
diff --git a/tests/gate/gate_test.cpp b/tests/gate/gate_test.cpp
index 811f624a..dd84dea3 100644
--- a/tests/gate/gate_test.cpp
+++ b/tests/gate/gate_test.cpp
@@ -336,7 +336,7 @@ using CComplex = std::complex<double>;
 //             }
 //         }
 //         matrix = std::cos(angle / 2) * Eigen::MatrixXcd::Identity(dim, dim) -
-//                  StdComplex(0, 1) * std::sin(angle / 2) * matrix;
+//                  StdComplex<double>(0, 1) * std::sin(angle / 2) * matrix;
 //         PauliOperator<double> pauli(target_vec, pauli_id_vec, 1.0);
 //         Gate pauli_gate = gate::PauliRotation(pauli, angle);
 //         pauli_gate->update_quantum_state(state);
@@ -410,7 +410,7 @@ void test_gate(Gate<double> gate_control,
     StateVector<double> state = StateVector<double>::Haar_random_state(n_qubits);
     auto amplitudes = state.get_amplitudes();
     StateVector<double> state_controlled(n_qubits - std::popcount(control_mask));
-    std::vector<Complex> amplitudes_controlled(state_controlled.dim());
+    std::vector<Complex<double>> 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];
diff --git a/tests/gate/merge_test.cpp b/tests/gate/merge_test.cpp
index 98234d16..d5c9afa6 100644
--- a/tests/gate/merge_test.cpp
+++ b/tests/gate/merge_test.cpp
@@ -13,78 +13,78 @@
 using namespace scaluq;
 
 TEST(GateTest, MergeGate) {
-    std::vector<Gate> gates;
+    std::vector<Gate<double>> gates;
     Random random;
     for (std::uint64_t target = 0; target < 2; target++) {
-        gates.push_back(gate::X(target));
-        gates.push_back(gate::Y(target));
-        gates.push_back(gate::Z(target));
-        gates.push_back(gate::H(target));
-        gates.push_back(gate::S(target));
-        gates.push_back(gate::Sdag(target));
-        gates.push_back(gate::T(target));
-        gates.push_back(gate::Tdag(target));
-        gates.push_back(gate::SqrtX(target));
-        gates.push_back(gate::SqrtXdag(target));
-        gates.push_back(gate::SqrtY(target));
-        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() * 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() * 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()),
-                                              Complex(random.uniform(), random.uniform())},
-                                   std::array{Complex(random.uniform(), random.uniform()),
-                                              Complex(random.uniform(), random.uniform())}}));
-        gates.push_back(gate::CX(target, target ^ 1));
-        gates.push_back(gate::CZ(target, target ^ 1));
-        gates.push_back(gate::Swap(target, target ^ 1));
+        gates.push_back(gate::X<double>(target));
+        gates.push_back(gate::Y<double>(target));
+        gates.push_back(gate::Z<double>(target));
+        gates.push_back(gate::H<double>(target));
+        gates.push_back(gate::S<double>(target));
+        gates.push_back(gate::Sdag<double>(target));
+        gates.push_back(gate::T<double>(target));
+        gates.push_back(gate::Tdag<double>(target));
+        gates.push_back(gate::SqrtX<double>(target));
+        gates.push_back(gate::SqrtXdag<double>(target));
+        gates.push_back(gate::SqrtY<double>(target));
+        gates.push_back(gate::SqrtYdag<double>(target));
+        gates.push_back(gate::P0<double>(target));
+        gates.push_back(gate::P1<double>(target));
+        gates.push_back(gate::RX<double>(target, random.uniform() * std::numbers::pi * 2));
+        gates.push_back(gate::RY<double>(target, random.uniform() * std::numbers::pi * 2));
+        gates.push_back(gate::RZ<double>(target, random.uniform() * std::numbers::pi * 2));
+        gates.push_back(gate::U1<double>(target, random.uniform() * std::numbers::pi * 2));
+        gates.push_back(gate::U2<double>(target,
+                                         random.uniform() * std::numbers::pi * 2,
+                                         random.uniform() * std::numbers::pi * 2));
+        gates.push_back(gate::U3<double>(target,
+                                         random.uniform() * std::numbers::pi * 2,
+                                         random.uniform() * std::numbers::pi * 2,
+                                         random.uniform() * std::numbers::pi * 2));
+        gates.push_back(gate::OneTargetMatrix<double>(
+            target,
+            {std::array{Complex<double>(random.uniform(), random.uniform()),
+                        Complex<double>(random.uniform(), random.uniform())},
+             std::array{Complex<double>(random.uniform(), random.uniform()),
+                        Complex<double>(random.uniform(), random.uniform())}}));
+        gates.push_back(gate::CX<double>(target, target ^ 1));
+        gates.push_back(gate::CZ<double>(target, target ^ 1));
+        gates.push_back(gate::Swap<double>(target, target ^ 1));
     }
-    gates.push_back(gate::I());
-    gates.push_back(gate::GlobalPhase(random.uniform() * std::numbers::pi * 2));
-    gates.push_back(
-        gate::TwoTargetMatrix(0,
-                              1,
-                              {std::array{Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform())},
-                               std::array{Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform())},
-                               std::array{Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform())},
-                               std::array{Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform()),
-                                          Complex(random.uniform(), random.uniform())}}));
-    gates.push_back(gate::Pauli(PauliOperator("X 0 Y 1", random.uniform())));
-    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() * 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));
+    gates.push_back(gate::I<double>());
+    gates.push_back(gate::GlobalPhase<double>(random.uniform() * std::numbers::pi * 2));
+    gates.push_back(gate::TwoTargetMatrix<double>(
+        0,
+        1,
+        {std::array{Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform())},
+         std::array{Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform())},
+         std::array{Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform())},
+         std::array{Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform()),
+                    Complex<double>(random.uniform(), random.uniform())}}));
+    gates.push_back(gate::Pauli<double>(PauliOperator<double>("X 0 Y 1", random.uniform())));
+    gates.push_back(gate::Pauli<double>(PauliOperator<double>("Z 0", random.uniform())));
+    gates.push_back(gate::Pauli<double>(PauliOperator<double>("Z 1", random.uniform())));
+    gates.push_back(gate::PauliRotation<double>(PauliOperator<double>("X 0 Y 1", random.uniform()),
+                                                random.uniform() * std::numbers::pi * 2));
+    gates.push_back(gate::PauliRotation<double>(PauliOperator<double>("Z 0", random.uniform()),
+                                                random.uniform() * std::numbers::pi * 2));
+    gates.push_back(gate::PauliRotation<double>(PauliOperator<double>("Z 1", random.uniform()),
+                                                random.uniform() * std::numbers::pi * 2));
     for (auto&& g1 : gates) {
         for (auto&& g2 : gates) {
             std::uint64_t n = 2;
-            auto state1 = StateVector::Haar_random_state(n);
+            auto state1 = StateVector<double>::Haar_random_state(n);
             auto state2 = state1.copy();
             auto [mg, phase] = merge_gate(g1, g2);
             g1->update_quantum_state(state1);
diff --git a/tests/gate/param_gate_test.cpp b/tests/gate/param_gate_test.cpp
index 854749e8..b7fdc02a 100644
--- a/tests/gate/param_gate_test.cpp
+++ b/tests/gate/param_gate_test.cpp
@@ -127,7 +127,7 @@ void test_gate(ParamGate<double> gate_control,
     StateVector state = StateVector<double>::Haar_random_state(n_qubits);
     auto amplitudes = state.get_amplitudes();
     StateVector<double> state_controlled(n_qubits - std::popcount(control_mask));
-    std::vector<Complex> amplitudes_controlled(state_controlled.dim());
+    std::vector<Complex<double>> 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];
diff --git a/tests/operator/test_operator.cpp b/tests/operator/test_operator.cpp
index cfe21852..a9485ec8 100644
--- a/tests/operator/test_operator.cpp
+++ b/tests/operator/test_operator.cpp
@@ -60,8 +60,8 @@ TEST(OperatorTest, CheckExpectationValue) {
         Eigen::VectorXcd test_state = Eigen::VectorXcd::Zero(dim);
         for (std::uint64_t i = 0; i < dim; ++i) test_state[i] = state_cp[i];
 
-        Complex res = rand_observable.get_expectation_value(state);
-        Complex test_res = (test_state.adjoint() * test_rand_observable * test_state)(0, 0);
+        Complex<double> res = rand_observable.get_expectation_value(state);
+        Complex<double> test_res = (test_state.adjoint() * test_rand_observable * test_state)(0, 0);
         ASSERT_NEAR(test_res.real(), res.real(), eps);
         ASSERT_NEAR(res.imag(), 0, eps);
         ASSERT_NEAR(test_res.imag(), 0, eps);
@@ -86,8 +86,9 @@ TEST(OperatorTest, CheckTransitionAmplitude) {
         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];
 
-        Complex res = rand_observable.get_transition_amplitude(state_bra, state_ket);
-        Complex test_res = (test_state_bra.adjoint() * test_rand_observable * test_state_ket)(0, 0);
+        Complex<double> res = rand_observable.get_transition_amplitude(state_bra, state_ket);
+        Complex<double> test_res =
+            (test_state_bra.adjoint() * test_rand_observable * test_state_ket)(0, 0);
         ASSERT_NEAR(test_res.real(), res.real(), eps);
         ASSERT_NEAR(test_res.imag(), res.imag(), eps);
     }
@@ -131,7 +132,7 @@ TEST(OperatorTest, MultiCoefTest) {
 
     for (std::uint64_t repeat = 0; repeat < 10; ++repeat) {
         auto op1 = generate_random_observable_with_eigen(n, random).first;
-        auto coef = Complex(random.normal(), random.normal());
+        auto coef = Complex<double>(random.normal(), random.normal());
         auto op = op1 * coef;
         auto state = StateVector<double>::Haar_random_state(n);
         auto exp1 = op1.get_expectation_value(state);
@@ -144,25 +145,25 @@ TEST(OperatorTest, ApplyToStateTest) {
     const std::uint64_t n_qubits = 3;
     StateVector<double> state_vector(n_qubits);
     state_vector.load([n_qubits] {
-        std::vector<Complex> tmp(1 << n_qubits);
-        for (std::uint64_t i = 0; i < tmp.size(); ++i) tmp[i] = Complex(i, 0);
+        std::vector<Complex<double>> tmp(1 << n_qubits);
+        for (std::uint64_t i = 0; i < tmp.size(); ++i) tmp[i] = Complex<double>(i, 0);
         return tmp;
     }());
 
     Operator<double> op(n_qubits);
-    op.add_operator({0b001, 0b010, Complex(2)});
+    op.add_operator({0b001, 0b010, Complex<double>(2)});
     op.add_operator({"X 2 Y 1", 1});
     op.apply_to_state(state_vector);
 
-    std::vector<Complex> expected = {
-        Complex(2, -6),
-        Complex(0, -7),
-        Complex(-6, 4),
-        Complex(-4, 5),
-        Complex(10, -2),
-        Complex(8, -3),
-        Complex(-14, 0),
-        Complex(-12, 1),
+    std::vector<Complex<double>> expected = {
+        Complex<double>(2, -6),
+        Complex<double>(0, -7),
+        Complex<double>(-6, 4),
+        Complex<double>(-4, 5),
+        Complex<double>(10, -2),
+        Complex<double>(8, -3),
+        Complex<double>(-14, 0),
+        Complex<double>(-12, 1),
     };
     ASSERT_EQ(state_vector.get_amplitudes(), expected);
 }
@@ -176,9 +177,9 @@ TEST(OperatorTest, Optimize) {
     op.add_operator(PauliOperator<double>("Z 1", 4.));
     op.add_operator(PauliOperator<double>("X 0 Y 1", 5.));
     op.optimize();
-    std::vector<std::pair<std::string, Complex>> expected = {
+    std::vector<std::pair<std::string, Complex<double>>> expected = {
         {"X 0 Y 1", 10.}, {"Y 0 Z 1", 2.}, {"Z 1", 7.}};
-    std::vector<std::pair<std::string, Complex>> test;
+    std::vector<std::pair<std::string, Complex<double>>> test;
     for (const auto& pauli : op.terms()) {
         test.emplace_back(pauli.get_pauli_string(), pauli.coef());
     }
diff --git a/tests/operator/test_pauli_operator.cpp b/tests/operator/test_pauli_operator.cpp
index f0cdabc4..6d05930f 100644
--- a/tests/operator/test_pauli_operator.cpp
+++ b/tests/operator/test_pauli_operator.cpp
@@ -192,13 +192,13 @@ TEST(PauliOperatorTest, ApplyToStateTest) {
     const std::uint64_t n_qubits = 3;
     StateVector<double> state_vector(n_qubits);
     state_vector.load([n_qubits] {
-        std::vector<Complex> tmp(1 << n_qubits);
-        for (std::uint64_t i = 0; i < tmp.size(); ++i) tmp[i] = Complex(i, 0);
+        std::vector<Complex<double>> tmp(1 << n_qubits);
+        for (std::uint64_t i = 0; i < tmp.size(); ++i) tmp[i] = Complex<double>(i, 0);
         return tmp;
     }());
 
-    PauliOperator<double> op(0b001, 0b010, Complex(2));
+    PauliOperator<double> op(0b001, 0b010, Complex<double>(2));
     op.apply_to_state(state_vector);
-    std::vector<Complex> expected = {2, 0, -6, -4, 10, 8, -14, -12};
+    std::vector<Complex<double>> expected = {2, 0, -6, -4, 10, 8, -14, -12};
     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 a395ecec..2b66e282 100644
--- a/tests/state/state_vector_batched_test.cpp
+++ b/tests/state/state_vector_batched_test.cpp
@@ -24,7 +24,7 @@ TEST(StateVectorBatchedTest, HaarRandomStateNorm) {
 TEST(StateVectorBatchedTest, LoadAndAmplitues) {
     const std::uint64_t batch_size = 4, n_qubits = 3;
     const std::uint64_t dim = 1 << n_qubits;
-    std::vector states_h(batch_size, std::vector<Complex>(dim));
+    std::vector states_h(batch_size, std::vector<Complex<double>>(dim));
     for (std::uint64_t b = 0; b < batch_size; ++b) {
         for (std::uint64_t i = 0; i < dim; ++i) {
             states_h[b][i] = b * dim + i;
@@ -45,7 +45,7 @@ TEST(StateVectorBatchedTest, OperateState) {
     const std::uint64_t batch_size = 4, n_qubits = 3;
     auto states = StateVectorBatched<double>::Haar_random_state(batch_size, n_qubits, false);
     auto states_add = StateVectorBatched<double>::Haar_random_state(batch_size, n_qubits, false);
-    const Complex coef(2.1, 3.5);
+    const Complex<double> coef(2.1, 3.5);
 
     auto states_cp = states.copy();
     for (std::uint64_t b = 0; b < batch_size; ++b) {
@@ -121,8 +121,8 @@ TEST(StateVectorBatchedTest, Entropy) {
 TEST(StateVectorBatchedTest, Sampling) {
     const std::uint64_t batch_size = 2, n_qubits = 3;
     StateVectorBatched<double> states(batch_size, n_qubits);
-    states.load(
-        std::vector<std::vector<Complex>>{{1, 4, 5, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 6, 4, 1}});
+    states.load(std::vector<std::vector<Complex<double>>>{{1, 4, 5, 0, 0, 0, 0, 0},
+                                                          {0, 0, 0, 0, 0, 6, 4, 1}});
     states.normalize();
     auto result = states.sampling(4096);
     std::vector cnt(2, std::vector<std::uint64_t>(states.dim(), 0));
diff --git a/tests/util/util.hpp b/tests/util/util.hpp
index 5282b54f..73df0df5 100644
--- a/tests/util/util.hpp
+++ b/tests/util/util.hpp
@@ -121,8 +121,8 @@ inline Eigen::MatrixXcd get_expanded_eigen_matrix_with_identity(
     const std::uint64_t right_dim = 1ULL << (qubit_count - target_qubit_index - 1);
     auto left_identity = Eigen::MatrixXcd::Identity(left_dim, left_dim);
     auto right_identity = Eigen::MatrixXcd::Identity(right_dim, right_dim);
-    return internal::kronecker_product(
-        internal::kronecker_product(right_identity, one_target_matrix), left_identity);
+    return internal::kronecker_product<double>(
+        internal::kronecker_product<double>(right_identity, one_target_matrix), left_identity);
 }
 
 // get expanded matrix
@@ -130,7 +130,8 @@ inline Eigen::MatrixXcd get_eigen_matrix_full_qubit_pauli(std::vector<std::uint6
     Eigen::MatrixXcd result = Eigen::MatrixXcd::Identity(1, 1);
     for (std::uint64_t i = 0; i < pauli_ids.size(); ++i) {
         result =
-            internal::kronecker_product(get_eigen_matrix_single_Pauli(pauli_ids[i]), result).eval();
+            internal::kronecker_product<double>(get_eigen_matrix_single_Pauli(pauli_ids[i]), result)
+                .eval();
     }
     return result;
 }