add expect_statevec

src/qfvm/qfvm.cpp

@@ -242,6 +242,24 @@ simulate_circuit_custate(py::object const& pycircuit,
 }
 #endif
 
+py::object expect_statevec(py::array_t<complex<double>> const&np_inputstate, py::list const paulis)
+{
+    py::buffer_info buf = np_inputstate.request();
+    auto* data_ptr = reinterpret_cast<std::complex<double>*>(buf.ptr);
+    size_t data_size = buf.size;
+    StateVector<double> state(data_ptr, buf.size);
+    py::list pyres;
+    for (auto pauli_h : paulis){
+         py::object pypauli = py::reinterpret_borrow<py::object>(pauli_h);
+         std::vector<pos_t> posv = pypauli.attr("pos").cast<std::vector<pos_t>>();
+         string paulistr = pypauli.attr("paulistr").cast<string>();
+        double expec = state.expect_pauli(paulistr, posv);
+        pyres.attr("append")(expec);
+    }
+    state.move_data_to_python();
+    return pyres;
+}
+
 PYBIND11_MODULE(qfvm, m) {
   m.doc() = "Qfvm simulator";
   m.def("simulate_circuit", &simulate_circuit, "Simulate with circuit",
@@ -252,6 +270,8 @@ PYBIND11_MODULE(qfvm, m) {
         "Simulate with circuit using clifford", py::arg("circuit"),
         py::arg("shots"));
 
+  m.def("expect_statevec", &expect_statevec, "Calculate paulis expectation", py::arg("inputstate"), py::arg("paulis"));
+
 #ifdef _USE_GPU
   m.def("simulate_circuit_gpu", &simulate_circuit_gpu, "Simulate with circuit",
         py::arg("circuit"),

src/qfvm/statevector.hpp

@@ -85,6 +85,9 @@ template <class real_t = double> class StateVector {
                                       uint control_num,
                                       RowMatrixXcd const& mat);
 
+  // Expectation and measurement
+  double expect_pauli(string paulistr, vector<pos_t> const& posv);
+
   // Measure and Reset
   std::pair<uint, double> sample_measure_probs(vector<pos_t> const& qbits);
   vector<double> probabilities() const;
@@ -1264,6 +1267,64 @@ template <typename T> void printVector(const std::vector<T>& vec) {
   std::cout << std::endl;
 }
 
+template <class real_t>
+double StateVector<real_t>::expect_pauli(string paulistr,
+                                         vector<pos_t> const& posv) {
+  size_t flip_mask = 0;
+  size_t z_mask = 0;
+  size_t y_phase_num = 0;
+  uint flip_q = 0;
+
+  for (uint i = 0; i < posv.size(); i++) {
+    uint q = posv[i];
+    switch (paulistr[i]) {
+    case 'I':
+      break;
+    case 'X': {
+      flip_mask += 1ll << q;
+      flip_q = std::max(flip_q, q);
+      break;
+    }
+    case 'Y': {
+      flip_mask += 1ll << q;
+      flip_q = std::max(flip_q, q);
+      y_phase_num++;
+      z_mask += 1ll << q;
+      break;
+    }
+    case 'Z': {
+      z_mask += 1ll << q;
+      break;
+    }
+    }
+  }
+
+  if (!flip_mask) {
+    size_t rsize = size_;
+    double val = 0.;
+#pragma omp paraleel for reduction(+ : val)
+    for (size_t j = 0; j < rsize; ++j) {
+      uint z_phase_num = Qfutil::popcount(j & z_mask) % 2;
+      int sign = 1 - 2 * z_phase_num;
+      val += (data_[j] * std::conj(data_[j])).real() * sign;
+    }
+    return val;
+  } else {
+    double val = 0.;
+    size_t rsize = size_ >> 1;
+#pragma omp parallel for reduction(+ : val)
+    for (size_t j = 0; j < rsize; ++j) {
+      size_t i0 = (j & ((1ll << flip_q) - 1)) | (j >> flip_q << flip_q << 1);
+      size_t i1 = i0 ^ flip_mask;
+      uint z_phase_num = Qfutil::popcount(i0 & z_mask) % 2;
+      uint total_phase_num = z_phase_num * 2 + y_phase_num;
+      std::complex<double> phase = Qfutil::PHASE_YZ[total_phase_num % 4];
+      val += 2. * (data_[i0] * std::conj(data_[i1]) * phase).real();
+    }
+    return val;
+  }
+}
+
 template <class real_t>
 std::pair<uint, double>
 StateVector<real_t>::sample_measure_probs(vector<pos_t> const& qbits) {

src/qfvm/util.h

@@ -5,7 +5,15 @@
 #include <regex>
 #include <time.h>
 #include <type_traits>
+#include <bitset>
+#include <random>
+#include <chrono>
+#include <set>
 #include <vector>
+#include "types.hpp"
+
+#define TICK(x) auto bench_##x = std::chrono::steady_clock::now();
+#define TOCK(x) std::cout << #x ": " << std::chrono::duration_cast<std::chrono::duration<double>>(std::chrono::steady_clock::now() - bench_##x).count() << "s" << std::endl;
 
 namespace Qfutil {
 
@@ -96,4 +104,18 @@ std::vector<real_t> find_numbers(const std::string& str) {
   return res;
 }
 
-} // namespace Qfutil
+
+/*----------------bit function------------------*/
+const std::complex<double> PHASE_YZ[4] = {1, imag_I, -1, -imag_I};
+
+inline static uint popcount(uint x) {
+    x = ((x & 0xaaaaaaaaaaaaaaaaUL) >> 1) + (x & 0x5555555555555555UL);
+    x = ((x & 0xccccccccccccccccUL) >> 2) + (x & 0x3333333333333333UL);
+    x = ((x & 0xf0f0f0f0f0f0f0f0UL) >> 4) + (x & 0x0f0f0f0f0f0f0f0fUL);
+    x = ((x & 0xff00ff00ff00ff00UL) >> 8) + (x & 0x00ff00ff00ff00ffUL);
+    x = ((x & 0xffff0000ffff0000UL) >> 16) + (x & 0x0000ffff0000ffffUL);
+    x = ((x & 0xffffffff00000000UL) >> 32) + (x & 0x00000000ffffffffUL);
+    return (uint)x;
+}
+
+}//namespace Qfutil