Update Simulator Release to check for entanglement

src/Qir/Runtime/lib/Simulators/FullstateSimulator.cpp

@@ -23,7 +23,8 @@
 
 
 #include "FloatUtils.hpp"
-#include "QirTypes.hpp" // TODO: Consider removing dependency on this file.
+#include "QirTypes.hpp"         // TODO: Consider removing dependency on this file.
+#include "QirRuntime.hpp"
 #include "QirRuntimeApi_I.hpp"
 #include "QSharpSimApi_I.hpp"
 #include "SimFactory.hpp"
@@ -225,11 +226,15 @@ namespace Quantum
 
         void ReleaseQubit(Qubit q) override
         {
-            typedef void (*TReleaseQubit)(unsigned, unsigned);
+            typedef bool (*TReleaseQubit)(unsigned, unsigned);
             static TReleaseQubit releaseQubit = reinterpret_cast<TReleaseQubit>(this->GetProc("release"));
 
-            releaseQubit(this->simulatorId, GetQubitId(q)); // Release qubit in the simulator.
-            qubitManager->Release(q);                       // Release it in the qubit manager.
+            // Release qubit in the simulator.
+            if (!releaseQubit(this->simulatorId, GetQubitId(q)))
+            {
+                quantum__rt__fail_cstr("Released qubits are entangled.");
+            }
+            qubitManager->Release(q); // Release it in the qubit manager.
         }
 
         Result Measure(long numBases, PauliId bases[], long numTargets, Qubit targets[]) override

src/Qir/Tests/FullstateSimulator/FullstateSimulatorTests.cpp

@@ -145,6 +145,9 @@ TEST_CASE("Fullstate simulator: ZZ measure", "[fullstate_simulator]")
     Result rOne = iqa->Measure(2, paulis, 2, q);
     REQUIRE(Result_One == sim->GetResultValue(rOne));
 
+    // Disentangle for release.
+    iqa->ControlledX(1, &q[0], q[1]);
+
     sim->ReleaseQubit(q[0]);
     sim->ReleaseQubit(q[1]);
 }
@@ -311,6 +314,10 @@ TEST_CASE("Fullstate simulator: exponents", "[fullstate_simulator]")
     // not crashes? consider it passing
     REQUIRE(true);
 
+    // Disentangle for release.
+    iqa->ControlledExp(2, qs, 3, paulis, &qs[2], -0.17);
+    iqa->Exp(2, paulis, qs, -0.42);
+
     for (int i = 0; i < n; i++)
     {
         sim->ReleaseQubit(qs[i]);
@@ -377,6 +384,9 @@ TEST_CASE("Fullstate simulator: get qubit state of Bell state", "[fullstate_simu
     REQUIRE(1.0 == Approx(norm).epsilon(0.0001));
     norm = 0.0;
 
+    // Disentangle to prepare for release.
+    iqa->ControlledX(1, &qs[0], qs[1]);
+
     for (int i = 0; i < n; i++)
     {
         sim->ReleaseQubit(qs[i]);

src/Simulation/Common/Exceptions/ReleasedQubitsAreEntangled.cs

@@ -0,0 +1,19 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Microsoft.Quantum.Simulation.Simulators.Exceptions
+{
+    public class ReleasedQubitsAreEntangled : Exception
+    {
+        public ReleasedQubitsAreEntangled()
+            : base("Released qubits are entangled.")
+        {
+        }
+    }
+}

src/Simulation/Common/Exceptions/ReleasedQubitsAreNotInZeroState.cs

@@ -1,4 +1,4 @@
-// Copyright (c) Microsoft Corporation. All rights reserved.
+// Copyright (c) Microsoft Corporation. All rights reserved.
 // Licensed under the MIT License.
 
 using System;
@@ -9,6 +9,7 @@
 
 namespace Microsoft.Quantum.Simulation.Simulators.Exceptions
 {
+    [Obsolete("This class is deprecated and will be removed in a future release.")]
     public class ReleasedQubitsAreNotInZeroState : Exception
     {
         public ReleasedQubitsAreNotInZeroState()

src/Simulation/Native/src/simulator/capi.cpp

@@ -69,9 +69,9 @@ extern "C"
         Microsoft::Quantum::Simulator::get(id)->allocateQubit(q);
     }
 
-    MICROSOFT_QUANTUM_DECL void release(_In_ unsigned id, _In_ unsigned q)
+    MICROSOFT_QUANTUM_DECL bool release(_In_ unsigned id, _In_ unsigned q)
     {
-        Microsoft::Quantum::Simulator::get(id)->release(q);
+        return Microsoft::Quantum::Simulator::get(id)->release(q);
     }
 
     MICROSOFT_QUANTUM_DECL unsigned num_qubits(_In_ unsigned id)

src/Simulation/Native/src/simulator/capi.hpp

@@ -63,14 +63,14 @@ extern "C"
     MICROSOFT_QUANTUM_DECL bool InjectState(
         _In_ unsigned sid,
         _In_ unsigned n,
-        _In_reads_(n) unsigned* q, // The listed qubits must be unentangled and in state |0>
+        _In_reads_(n) unsigned* q, // The listed qubits must be disentangled and in state |0>
         _In_ double* re, // 2^n real parts of the amplitudes of the superposition the listed qubits should be put into
         _In_ double* im  // 2^n imaginary parts of the amplitudes
     );
 
     // allocate and release
     MICROSOFT_QUANTUM_DECL void allocateQubit(_In_ unsigned sid, _In_ unsigned qid); // NOLINT
-    MICROSOFT_QUANTUM_DECL void release(_In_ unsigned sid, _In_ unsigned q); // NOLINT
+    MICROSOFT_QUANTUM_DECL bool release(_In_ unsigned sid, _In_ unsigned q); // NOLINT
     MICROSOFT_QUANTUM_DECL unsigned num_qubits(_In_ unsigned sid); // NOLINT
 
     // single-qubit gates

src/Simulation/Native/src/simulator/kernels.hpp

@@ -144,6 +144,63 @@ bool isclassical(
     return true;
 }
 
+template <class T, class A>
+std::pair<bool, bool> is_classical_or_entangled(
+    std::vector<std::complex<T>, A> const& wfn,
+    std::size_t q,
+    T eps = 100. * std::numeric_limits<T>::epsilon())
+{
+    // Iterate through the wave function using offstet and mask to check for both whether the qubit
+    // is classical (in the |0⟩ or |1⟩ state) or is entangled with any other qubit. By checking states
+    // where the target qubit is the only difference, the "isclassical" check can speed up the iteration
+    // which must check the whole vector. If the target qubit has a nonzero probability of being measured
+    // as |0⟩ (variable "have0") AND nonzero probability of being measured as |1⟩ (variable "have1") then
+    // we know it is not classical with regard to the computational basis. Further, if the target qubit
+    // is not classical with regard to the computational basis AND for any two states where the only difference
+    // in the states is the target qubit (ie: |101⟩ and |100⟩ for qubit 0) those have nonzero probability, then
+    // the qubit cannot be entangled.
+    std::size_t offset = 1ull << q;
+    bool have0 = false;
+    bool have1 = false;
+    bool notentangled = false;
+
+    std::size_t maski = ~(offset - 1);
+    // When iterating below, the entry in the wave function with index `i + j` and the entry with index
+    // `i + j + offset` represent two states where the only difference is the measured value of the target
+    // qubit. For example, if the target qubit has id 2 (offset = 4) in wave function of size 4, then for
+    // i = 0 and j = 0 the two states checked are 0 and 4, or |0000⟩ and |0100⟩. When i = 8 and j = 1, the
+    // two states are 9 and 13, or |1001⟩ and |1101⟩.
+#ifndef _MSC_VER
+#pragma omp parallel for schedule(static) reduction(|| : have0, have1, notentangled)
+    for (std::intptr_t i = 0; i < static_cast<std::intptr_t>(wfn.size()); i += 2 * offset)
+        for (std::intptr_t j = 0; j < static_cast<std::intptr_t>(offset); ++j)
+        {
+            bool has0 = std::norm(wfn[i + j]) >= eps;
+            bool has1 = std::norm(wfn[i + j + offset]) >= eps;
+            have0 = have0 || has0;
+            have1 = have1 || has1;
+            notentangled = notentangled || (has0 && has1);
+        }
+#else
+#pragma omp parallel for schedule(static) reduction(|| : have0, have1, notentangled)
+    for (std::intptr_t l = 0; l < static_cast<std::intptr_t>(wfn.size()) / 2; l++)
+    {
+        std::intptr_t j = l % offset;
+        std::intptr_t i = ((l & maski) << 1);
+        bool has0 = std::norm(wfn[i + j]) >= eps;
+        bool has1 = std::norm(wfn[i + j + offset]) >= eps;
+        have0 = have0 || has0;
+        have1 = have1 || has1;
+        notentangled = notentangled || (has0 && has1);
+    }
+#endif
+
+    // isclassical = true IFF have0 XOR have1
+    // isentangled = true IFF have0 AND have1 AND for all pairs of states where only the target qubit
+    //               differs one of those states has zero probability.
+    return std::make_pair<bool, bool>(have0 ^ have1, have0 && have1 && !notentangled);
+}
+
 template <class T, class A>
 double jointprobability(std::vector<T, A> const& wfn, std::vector<unsigned> const& qs, bool val = true)
 {

src/Simulation/Native/src/simulator/local_test.cpp

@@ -569,7 +569,7 @@ TEST_CASE("Should fail to inject state if qubits aren't all |0>", "[local_test]"
 
     std::vector<ComplexType> amplitudes_sub = {{amp, 0.0}, {amp, 0.0}, {amp, 0.0}, {0.0, 0.0}};
 
-    // unentangled but not |0>
+    // not entangled but not |0>
     sim.H(qs[1]);
     REQUIRE_FALSE(sim.InjectState(qs, amplitudes));
     REQUIRE_FALSE(sim.InjectState({qs[0], qs[1]}, amplitudes_sub));
@@ -828,6 +828,82 @@ TEST_CASE("test_multicontrol", "[local_test]")
     }
 }
 
+TEST_CASE("test_is_classical_or_entangled", "[local_test]")
+{
+    SimulatorType sim;
+    auto qbits = sim.allocate(4);
+
+    CHECK(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.X(qbits[0]);
+    sim.X(qbits[2]);
+    CHECK(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.X(qbits[0]);
+    sim.H(qbits[0]);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.H(qbits[1]);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.H(qbits[1]);
+    sim.CX(qbits[0], qbits[1]);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.CX(qbits[0], qbits[2]);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.CX(qbits[0], qbits[1]);
+    sim.CX(qbits[0], qbits[2]);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.H(qbits[0]);
+    CHECK(sim.is_classical_or_entangled(qbits[0]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[0]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[1]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[1]).second);
+    CHECK(sim.is_classical_or_entangled(qbits[2]).first);
+    CHECK_FALSE(sim.is_classical_or_entangled(qbits[2]).second);
+
+    sim.release(qbits);
+    CHECK(sim.num_qubits() == 0);
+}
+
 void test_extract_qubits_state_simple(int qubits_number)
 {
     SimulatorType sim;