Fix deterministic measure of stabilizer

releasenotes/notes/fix_stabilizer_measure-a06d761eba2546d2.yaml

@@ -0,0 +1,5 @@
+---
+fixes:
+  - |
+    Fix of deterministic measure of stabilizer and change of OpenMP
+    parallelization of loop

src/simulators/stabilizer/clifford.hpp

@@ -456,9 +456,8 @@ bool Clifford::measure_and_update(const uint64_t qubit,
 
           t0 = rX & destabilizer_table_[q].Z(i);
           t1 = rZ ^ destabilizer_table_[q].X(i);
-          t1 ^= t0;
 
-          d1 ^= (t0 & d0);
+          d1 ^= (t0 & (~d0));
           d0 ^= t0;
           d1 ^= (t0 & t1);
 
@@ -475,9 +474,8 @@ bool Clifford::measure_and_update(const uint64_t qubit,
 
           t0 = rX & stabilizer_table_[q].Z(i);
           t1 = rZ ^ stabilizer_table_[q].X(i);
-          t1 ^= t0;
 
-          s1 ^= (t0 & s0);
+          s1 ^= (t0 & (~s0));
           s0 ^= t0;
           s1 ^= (t0 & t1);
 
@@ -515,127 +513,79 @@ bool Clifford::measure_and_update(const uint64_t qubit,
   } else {
     // Deterministic outcome
     uint_t outcome = 0;
-    Pauli::Pauli<BV::BinaryVector> accum(num_qubits_);
     uint_t blocks = destabilizer_phases_.blockLength();
 
-    if (blocks < 2) {
+    auto measure_determinisitic_func = [this, blocks, qubit](AER::int_t q) {
+      uint_t accumX = 0;
+      uint_t accumZ = 0;
+      uint_t exponent_l = 0ull;
+      uint_t exponent_h = 0ull;
+      bool accumX_prev = false;
+      bool accumZ_prev = false;
+
       for (uint_t ib = 0; ib < blocks; ib++) {
+        uint_t tl, th, add;
         uint_t destabilizer_mask = destabilizer_table_[qubit].X(ib);
-        uint_t exponent_l = 0ull;
-        uint_t exponent_h = 0ull;
-
-        for (uint_t q = 0; q < num_qubits_; q++) {
-          uint_t tl, th, add;
-          uint_t accumX = 0ull - (uint_t)accum.X[q];
-          uint_t accumZ = 0ull - (uint_t)accum.Z[q];
-
-          tl = accumX & stabilizer_table_[q].Z(ib);
-          th = accumZ ^ stabilizer_table_[q].X(ib);
-
-          add = tl & exponent_l;
-          exponent_l ^= tl;
-          exponent_h ^= add;
-          exponent_h ^= (tl & th);
-
-          tl = stabilizer_table_[q].X(ib) & accumZ;
-          th = stabilizer_table_[q].Z(ib) ^ accumX;
-          th ^= tl;
-
-          add = tl & exponent_l;
-          exponent_l ^= tl;
-          exponent_h ^= add;
-          exponent_h ^= (tl & th);
-
-          add = stabilizer_table_[q].X(ib) & destabilizer_mask;
-          accumX &= AER::Utils::popcount(add) & 1;
-          add = stabilizer_table_[q].Z(ib) & destabilizer_mask;
-          accumZ &= AER::Utils::popcount(add) & 1;
-
-          accum.X.setValue((bool)accumX, q);
-          accum.Z.setValue((bool)accumZ, q);
-        }
-        exponent_h ^= stabilizer_phases_(ib);
-        outcome ^= (exponent_h & destabilizer_mask);
 
-        if ((exponent_l & destabilizer_mask) != 0) {
-          throw std::runtime_error("Clifford: rowsum error");
-        }
-      }
-    } else {
-      uint_t blockSize = destabilizer_phases_.blockSize();
-
-      // loop for cache blocking
-      for (uint_t ii = 0; ii < blocks; ii++) {
-        uint_t destabilizer_mask = destabilizer_table_[qubit].X(ii);
-        if (destabilizer_mask == 0)
-          continue;
-
-        uint_t exponent_l = 0;
-        uint_t exponent_lc = 0;
-        uint_t exponent_h = 0;
-
-        auto measure_determinisitic_func =
-            [this, &accum, &exponent_l, &exponent_lc, &exponent_h, blocks,
-             blockSize, destabilizer_mask, ii](AER::int_t qq) {
-              uint_t qs = qq * blockSize;
-              uint_t qe = qs + blockSize;
-              if (qe > num_qubits_)
-                qe = num_qubits_;
-
-              uint_t local_exponent_l = 0;
-              uint_t local_exponent_h = 0;
-
-              for (uint_t q = qs; q < qe; q++) {
-                uint_t sX = stabilizer_table_[q].X(ii);
-                uint_t sZ = stabilizer_table_[q].Z(ii);
-
-                uint_t accumX = (0ull - (uint_t)accum.X[q]);
-                uint_t accumZ = (0ull - (uint_t)accum.Z[q]);
-
-                // exponents for this block
-                uint_t t0, t1;
-
-                t0 = accumX & sZ;
-                t1 = accumZ ^ sX;
-
-                local_exponent_h ^= (t0 & local_exponent_l);
-                local_exponent_l ^= t0;
-                local_exponent_h ^= (t0 & t1);
-
-                t0 = sX & accumZ;
-                t1 = sZ ^ accumX;
-                t1 ^= t0;
-
-                local_exponent_h ^= (t0 & local_exponent_l);
-                local_exponent_l ^= t0;
-                local_exponent_h ^= (t0 & t1);
-
-                // update accum
-                accumX &= AER::Utils::popcount(sX & destabilizer_mask) & 1;
-                accum.X.setValue((accumX != 0), q);
-                accumZ &= AER::Utils::popcount(sZ & destabilizer_mask) & 1;
-                accum.Z.setValue((accumZ != 0), q);
-              }
-
-#pragma omp atomic
-              exponent_lc |= local_exponent_l;
-#pragma omp atomic
-              exponent_l ^= local_exponent_l;
-#pragma omp atomic
-              exponent_h ^= local_exponent_h;
-            };
-        AER::Utils::apply_omp_parallel_for(
-            (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0,
-            blocks, measure_determinisitic_func, omp_threads_);
-
-        // if exponent_l is 0 and any of local_exponent_l is
-        // 1, then flip exponent_h
-        exponent_h ^= (exponent_lc ^ exponent_l);
-        exponent_h ^= stabilizer_phases_(ii);
-        outcome ^= (exponent_h & destabilizer_mask);
+        uint_t sX = stabilizer_table_[q].X(ib) & destabilizer_mask;
+        uint_t sZ = stabilizer_table_[q].Z(ib) & destabilizer_mask;
+
+        // accumulate for 64 bits block
+        accumX = sX ^ (uint_t)accumX_prev;
+        accumZ = sZ ^ (uint_t)accumZ_prev;
+        accumX ^= (accumX << 1);
+        accumZ ^= (accumZ << 1);
+        accumX ^= (accumX << 2);
+        accumZ ^= (accumZ << 2);
+        accumX ^= (accumX << 4);
+        accumZ ^= (accumZ << 4);
+        accumX ^= (accumX << 8);
+        accumZ ^= (accumZ << 8);
+        accumX ^= (accumX << 16);
+        accumZ ^= (accumZ << 16);
+        accumX ^= (accumX << 32);
+        accumZ ^= (accumZ << 32);
+        // store for next iteration
+        accumX_prev = ((accumX >> 63) & 1) != 0;
+        accumZ_prev = ((accumZ >> 63) & 1) != 0;
+        // correct for this iteration
+        accumX ^= sX;
+        accumZ ^= sZ;
+        accumX &= destabilizer_mask;
+        accumZ &= destabilizer_mask;
+
+        tl = accumX & sZ;
+        th = accumZ ^ sX;
+
+        add = tl & exponent_l;
+        exponent_l ^= tl;
+        exponent_h ^= add;
+        exponent_h ^= (tl & th);
+
+        tl = sX & accumZ;
+        th = sZ ^ accumX;
+
+        add = tl & (~exponent_l);
+        exponent_l ^= tl;
+        exponent_h ^= add;
+        exponent_h ^= (tl & th);
       }
+      // convert 2-bits x 64 integer into bit count here
+      return AER::Utils::popcount(exponent_h) * 2 +
+             AER::Utils::popcount(exponent_l);
+    };
+
+    outcome = AER::Utils::apply_omp_parallel_for_reduction_int(
+        (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0,
+        num_qubits_, measure_determinisitic_func, omp_threads_);
+
+    uint_t stab_h = 0ull;
+    for (uint_t ib = 0; ib < blocks; ib++) {
+      stab_h ^= (destabilizer_table_[qubit].X(ib) & stabilizer_phases_(ib));
     }
-    return ((AER::Utils::popcount(outcome) & 1) != 0);
+    outcome += AER::Utils::popcount(stab_h) * 2;
+
+    return ((outcome & 3) == 2);
   }
 }
 

test/terra/backends/aer_simulator/test_measure.py

@@ -284,6 +284,38 @@ def test_measure_sampling_large_ghz_stabilizer(self, method, device, num_qubits)
         targets["1" * num_qubits] = shots / 2
         self.assertDictAlmostEqual(counts, targets, delta=delta * shots)
 
+    @supported_methods(["stabilizer"])
+    def test_measure_stablizer_deterministic(self, method, device):
+        """Test stabilizer measure for deterministic case"""
+        backend = self.backend(method=method, device=device)
+        shots = 10000
+        qc = QuantumCircuit(5, 1)
+        qc.h([2, 4])
+        qc.cx(2, 0)
+        qc.s(0)
+        qc.cx(4, 2)
+        qc.h(0)
+        qc.cx(2, 3)
+        qc.s(4)
+        qc.cx(1, 0)
+        qc.h([3, 4])
+        qc.cx(3, 2)
+        qc.h(3)
+        qc.cx(0, 3)
+        qc.cx(3, 1)
+        qc.s(0)
+        qc.s(1)
+        qc.h(0)
+        qc.s(0)
+        qc.cx(4, 0)
+        qc.cx(0, 1)
+        qc.measure(1, 0)
+        result = backend.run(qc, shots=shots).result()
+        counts = result.get_counts()
+        self.assertSuccess(result)
+
+        self.assertDictEqual(counts, {"1": shots})
+
     # ---------------------------------------------------------------------
     # Test MPS algorithms for measure
     # ---------------------------------------------------------------------