From ec8d6e185b15f7a92010f15500d1c8462d4c8e00 Mon Sep 17 00:00:00 2001
From: Dan Strano <stranoj@gmail.com>
Date: Sun, 30 Jul 2023 18:59:50 -0400
Subject: [PATCH] Debug rounding in RdmCloneFlush()

---
 include/qstabilizerhybrid.hpp | 28 +++-------------------------
 src/qstabilizerhybrid.cpp     | 14 ++++++++++----
 2 files changed, 13 insertions(+), 29 deletions(-)

diff --git a/include/qstabilizerhybrid.hpp b/include/qstabilizerhybrid.hpp
index fd0b3188f..514f76207 100644
--- a/include/qstabilizerhybrid.hpp
+++ b/include/qstabilizerhybrid.hpp
@@ -264,29 +264,7 @@ class QStabilizerHybrid : public QParity, public QInterface {
             shard->gate[3U] = complex(angleCos, angleSin);
         }
 
-        const complex h[4U] = { SQRT1_2_R1, SQRT1_2_R1, SQRT1_2_R1, -SQRT1_2_R1 };
-        for (size_t i = qubitCount; i < shards.size(); ++i) {
-            // Flush all buffers as close as possible to Clifforrd.
-            const MpsShardPtr& shard = shards[i];
-            shard->Compose(h);
-
-            const real1 angle = (real1)(FractionalRzAngleWithFlush(i, std::arg(shard->gate[3U] / shard->gate[0U])) / 2);
-            if ((2 * abs(angle) / PI_R1) <= FP_NORM_EPSILON) {
-                stabilizer->H(i);
-                stabilizer->ForceM(i, false);
-                stabilizer->Dispose(i, 1U);
-                shards.erase(shards.begin() + i);
-                --ancillaCount;
-
-                continue;
-            }
-            const real1 angleCos = cos(angle);
-            const real1 angleSin = sin(angle);
-            shard->gate[0U] = complex(angleCos, -angleSin);
-            shard->gate[3U] = complex(angleCos, angleSin);
-
-            shard->Compose(h);
-        }
+        RdmCloneFlush();
     }
 
     void CombineAncillae();
@@ -295,11 +273,11 @@ class QStabilizerHybrid : public QParity, public QInterface {
     {
         CombineAncillae();
         QStabilizerHybridPtr clone = std::dynamic_pointer_cast<QStabilizerHybrid>(Clone());
-        clone->RdmCloneFlush();
+        clone->RdmCloneFlush(ONE_R1 / 4);
 
         return clone;
     }
-    void RdmCloneFlush(real1_f threshold = (ONE_R1 / 4));
+    void RdmCloneFlush(real1_f threshold = FP_NORM_EPSILON);
 
     real1_f ApproxCompareHelper(
         QStabilizerHybridPtr toCompare, bool isDiscreteBool, real1_f error_tol = TRYDECOMPOSE_EPSILON);
diff --git a/src/qstabilizerhybrid.cpp b/src/qstabilizerhybrid.cpp
index 69cd589f4..7741ed77e 100644
--- a/src/qstabilizerhybrid.cpp
+++ b/src/qstabilizerhybrid.cpp
@@ -1762,14 +1762,12 @@ void QStabilizerHybrid::CombineAncillae()
         return;
     }
 
-    RdmCloneFlush(FP_NORM_EPSILON);
+    FlushCliffordFromBuffers();
 
     if (!ancillaCount) {
         return;
     }
 
-    FlushCliffordFromBuffers();
-
     // The ancillae sometimes end up in a configuration where measuring an earlier ancilla collapses a later ancilla.
     // If so, we can combine (or cancel) the phase effect on the earlier ancilla and completely separate the later.
     // We must preserve the earlier ancilla's entanglement, besides partial collapse with the later ancilla.
@@ -1865,7 +1863,6 @@ void QStabilizerHybrid::CombineAncillae()
 
 void QStabilizerHybrid::RdmCloneFlush(real1_f threshold)
 {
-    FlushCliffordFromBuffers();
     const complex h[4U] = { SQRT1_2_R1, SQRT1_2_R1, SQRT1_2_R1, -SQRT1_2_R1 };
     for (size_t i = shards.size() - 1U; i >= qubitCount; --i) {
         MpsShardPtr& shard = shards[i];
@@ -1875,6 +1872,15 @@ void QStabilizerHybrid::RdmCloneFlush(real1_f threshold)
         for (int p = 0; p < 2; ++p) {
             shard->Compose(h);
             QStabilizerHybridPtr clone = std::dynamic_pointer_cast<QStabilizerHybrid>(Clone());
+
+            if (clone->stabilizer->IsSeparable(i)) {
+                stabilizer->Dispose(i, 1U);
+                shards.erase(shards.begin() + i);
+                --ancillaCount;
+
+                break;
+            }
+
             clone->stabilizer->H(i);
             clone->stabilizer->ForceM(i, p == 1);