Optimize QBdt::TryDecompose()

unitaryfund · Nov 5, 2024 · 1ed2940 · 1ed2940
1 parent 1471c8b
commit 1ed2940
Show file tree

Hide file tree

Showing 3 changed files with 36 additions and 26 deletions.
diff --git a/examples/qbdd_separability.cpp b/examples/qbdd_separability.cpp
@@ -120,8 +120,10 @@ int main(int argc, char* argv[])
 
     for (bitLenInt i = 0U; i < halfWidth; ++i) {
         const bitLenInt subsystemSize = i + 1U;
-        Qrack::QBdtPtr qsim = std::dynamic_pointer_cast<Qrack::QBdt>(Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, width, Qrack::ZERO_BCI));
-        Qrack::QBdtPtr qsim_b = std::dynamic_pointer_cast<Qrack::QBdt>(Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, width - subsystemSize, Qrack::ZERO_BCI));
+        Qrack::QBdtPtr qsim = std::dynamic_pointer_cast<Qrack::QBdt>(
+            Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, width, Qrack::ZERO_BCI));
+        Qrack::QBdtPtr qsim_b = std::dynamic_pointer_cast<Qrack::QBdt>(
+            Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, width - subsystemSize, Qrack::ZERO_BCI));
 
         ghz(qsim);
 
@@ -147,8 +149,10 @@ int main(int argc, char* argv[])
              << " qubit subsystem, nonseparable: " << duration << " seconds, " << (result ? "failure." : "success.")
              << endl;
 
-        qsim = std::dynamic_pointer_cast<Qrack::QBdt>(Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, subsystemSize, Qrack::ZERO_BCI));
-        qsim_b = std::dynamic_pointer_cast<Qrack::QBdt>(Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, width - subsystemSize, Qrack::ZERO_BCI));
+        qsim = std::dynamic_pointer_cast<Qrack::QBdt>(
+            Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, subsystemSize, Qrack::ZERO_BCI));
+        qsim_b = std::dynamic_pointer_cast<Qrack::QBdt>(
+            Qrack::CreateQuantumInterface({ Qrack::QINTERFACE_BDT }, width - subsystemSize, Qrack::ZERO_BCI));
 
         ghz(qsim);
         ghz(qsim_b);

diff --git a/include/qbdt.hpp b/include/qbdt.hpp
@@ -274,6 +274,33 @@ class QBdt : public QParity, public QInterface {
      */
     bool IsSeparable(bitLenInt start);
 
+    bool TryDecompose(bitLenInt start, QInterfacePtr dest, real1_f error_tol = TRYDECOMPOSE_EPSILON)
+    {
+        if (error_tol > TRYDECOMPOSE_EPSILON) {
+            return QInterface::TryDecompose(start, dest, error_tol);
+        }
+
+        const bitLenInt length = dest->GetQubitCount();
+        const bitLenInt nStart = qubitCount - length;
+        const bitLenInt shift = nStart - start;
+        for (bitLenInt i = 0U; i < shift; ++i) {
+            Swap(start + i, qubitCount - (i + 1U));
+        }
+
+        const bool isSeparable = IsSeparable(nStart);
+
+        for (bitLenInt i = shift; i > 0U; --i) {
+            Swap(start + (i - 1U), qubitCount - i);
+        }
+
+        if (isSeparable) {
+            Decompose(start, dest);
+            return true;
+        }
+
+        return false;
+    }
+
     using QInterface::Compose;
     bitLenInt Compose(QBdtPtr toCopy, bitLenInt start);
     bitLenInt Compose(QInterfacePtr toCopy, bitLenInt start)

diff --git a/src/qbdt/tree.cpp b/src/qbdt/tree.cpp
@@ -297,35 +297,14 @@ bool QBdt::IsSeparable(bitLenInt start)
             "QBdt::IsSeparable() start parameter must be at least 1 and less than the QBdt qubit width!");
     }
 
-#if 0
-    // Premature optimization, but works in principle if Swap() is inexpensive:
-    if ((start << 1U) > qubitCount) {
-        const bitLenInt nStart = qubitCount - start;
-        for (size_t i = 0U; i < nStart; ++i) {
-            Swap(i, qubitCount - (i + 1U));
-        }
-        const bool result = IsSeparable(nStart);
-        for (size_t i = 0U; i < nStart; ++i) {
-            Swap(i, qubitCount - (i + 1U));
-        }
-
-        return result;
-    }
-#endif
-
     FlushBuffers();
 
     // If the tree has been fully reduced, this should ALWAYS be the same for ALL branches
     // (that have nonzero amplitude), if-and-only-if the state is separable.
-    // QBdtNodeInterface* subsystemPtr = NULL;
-
-    // Alternatively, using shared_ptr with our custom (non-)equality operator will reduce
-    // the tree in-flight if it is not already fully-reduced, at no significant additional
-    // overhead if the tree is already fully-reduced.
     QBdtNodeInterfacePtr subsystemPtr = NULL;
 
-    bool result = true;
     const bitCapInt qPower = pow2(start);
+    bool result = true;
 
     par_for_qbdt(
         qPower, start,