added non-unitary Pauli gadgets (#637)

diogomaia00 · TysonRayJones · web-flow · commit 5fa4f609847a · 2025-06-04T16:29:19.000+02:00
as part of unitaryHACK 2025, challenge issue #594 --------- Co-authored-by: Tyson Jones <tyson.jones.input@gmail.com>
diff --git a/AUTHORS.txt b/AUTHORS.txt
@@ -44,6 +44,8 @@ Dr Ian Bush [consultant]
     HPC
 
 External contributors:
+Diogo Pratas Maia
+    added non-unitary Pauli gadget (for unitaryHACK issue #594)
 Mai Đức Khang
     implemented RAM probe (for unitaryHACK issue #600)
 James Richings
diff --git a/README.md b/README.md
@@ -257,6 +257,7 @@ See the [docs](docs/README.md) for enabling acceleration and running the unit te
 
 In addition to QuEST's [authors](AUTHORS.txt), we sincerely thank the following external contributors to QuEST.
 
+- [Diogo Pratas Maia](https://github.com/diogomaia00) for implementing non-unitary Pauli gadgets (unitaryHACK 2025 [#594](https://github.com/QuEST-Kit/QuEST/issues/594)).
 - [Mai Đức Khang](https://github.com/Roll249) for implementing a RAM probe (unitaryHACK 2025 [#600](https://github.com/QuEST-Kit/QuEST/issues/600)). 
 - [James Richings](https://github.com/JPRichings) for patching a v4 overflow bug.
 - [Luc Jaulmes](https://github.com/lucjaulmes) for patching v4's CMake installation.
diff --git a/examples/automated/apply_no_unitary_pauli_gadget.c b/examples/automated/apply_no_unitary_pauli_gadget.c
@@ -0,0 +1,18 @@
+#include "quest.h"
+
+int main() {
+    initQuESTEnv();
+
+    Qureg qureg = createQureg(3);
+    PauliStr str = getInlinePauliStr("XYZ", {0,1,2});
+    qcomp angle = getQcomp(.4, .8);
+
+    initPlusState(qureg);
+    applyNonUnitaryPauliGadget(qureg, str, angle);
+
+    qreal norm = calcTotalProb(qureg);
+    reportScalar("norm", norm);
+
+    finalizeQuESTEnv();
+    return 0;
+}
diff --git a/quest/include/operations.h b/quest/include/operations.h
@@ -5,6 +5,7 @@
  * instead exposed in decoherence.h
  * 
  * @author Tyson Jones
+ * @author Diogo Pratas Maia (non-unitary Pauli gadget)
  * 
  * @defgroup operations Operations
  * @ingroup api
@@ -1853,6 +1854,8 @@ void multiplyPauliGadget(Qureg qureg, PauliStr str, qreal angle);
  */
 void applyPauliGadget(Qureg qureg, PauliStr str, qreal angle);
 
+/// @notyetdoced
+void applyNonUnitaryPauliGadget(Qureg qureg, PauliStr str, qcomp angle);
 
 /// @notyetdoced
 void applyControlledPauliGadget(Qureg qureg, int control, PauliStr str, qreal angle);
diff --git a/quest/src/api/operations.cpp b/quest/src/api/operations.cpp
@@ -1432,6 +1432,22 @@ void applyPauliGadget(Qureg qureg, PauliStr str, qreal angle) {
     applyMultiStateControlledPauliGadget(qureg, nullptr, nullptr, 0, str, angle);
 }
 
+void applyNonUnitaryPauliGadget(Qureg qureg, PauliStr str, qcomp angle) {
+    validate_quregFields(qureg, __func__);
+    validate_pauliStrTargets(qureg, str, __func__);
+
+    qcomp phase = util_getPhaseFromGateAngle(angle);
+    localiser_statevec_anyCtrlPauliGadget(qureg, {}, {}, str, phase);
+
+    if (!qureg.isDensityMatrix)
+        return;
+
+    // conj(e^i(a)XZ) = e^(-i conj(a)XZ) but conj(Y)=-Y, so odd-Y undoes phase negation
+    phase = std::conj(phase) * (paulis_hasOddNumY(str) ? 1 : -1);
+    str = paulis_getShiftedPauliStr(str, qureg.numQubits);
+    localiser_statevec_anyCtrlPauliGadget(qureg, {}, {}, str, phase);
+}
+
 void applyControlledPauliGadget(Qureg qureg, int control, PauliStr str, qreal angle) {
     validate_quregFields(qureg, __func__);
     validate_controlAndPauliStrTargets(qureg, control, str, __func__);
diff --git a/quest/src/core/localiser.cpp b/quest/src/core/localiser.cpp
@@ -1244,7 +1244,7 @@ extern int paulis_getPrefixZSign(Qureg qureg, vector<int> prefixZ) ;
 extern qcomp paulis_getPrefixPaulisElem(Qureg qureg, vector<int> prefixY, vector<int> prefixZ);
 
 
-void anyCtrlZTensorOrGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, bool isGadget, qreal phase) {
+void anyCtrlZTensorOrGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, bool isGadget, qcomp phase) {     
     assertValidCtrlStates(ctrls, ctrlStates);
     setDefaultCtrlStates(ctrls, ctrlStates);
 
@@ -1339,14 +1339,14 @@ void localiser_statevec_anyCtrlPauliTensor(Qureg qureg, vector<int> ctrls, vecto
 }
 
 
-void localiser_statevec_anyCtrlPhaseGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qreal phase) {
+void localiser_statevec_anyCtrlPhaseGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qcomp phase) {
 
     bool isGadget = true;
-    anyCtrlZTensorOrGadget(qureg, ctrls, ctrlStates, targs, isGadget, phase); 
+    anyCtrlZTensorOrGadget(qureg, ctrls, ctrlStates, targs, isGadget, phase);
 }
 
 
-void localiser_statevec_anyCtrlPauliGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, PauliStr str, qreal phase) {
+void localiser_statevec_anyCtrlPauliGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, PauliStr str, qcomp phase) {
 
     // when str=IZ, we must use the above bespoke algorithm
     if (!paulis_containsXOrY(str)) {
diff --git a/quest/src/core/localiser.hpp b/quest/src/core/localiser.hpp
@@ -118,9 +118,9 @@ void localiser_statevec_anyCtrlAnyTargAnyMatr(Qureg qureg, vector<int> ctrls, ve
 
 void localiser_statevec_anyCtrlPauliTensor(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, PauliStr str, qcomp globalFactor=1);
 
-void localiser_statevec_anyCtrlPauliGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, PauliStr str, qreal phase);
+void localiser_statevec_anyCtrlPauliGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, PauliStr str, qcomp phase);
 
-void localiser_statevec_anyCtrlPhaseGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qreal phase);
+void localiser_statevec_anyCtrlPhaseGadget(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qcomp phase);
 
 
 /*
diff --git a/quest/src/core/memory.cpp b/quest/src/core/memory.cpp
@@ -10,6 +10,7 @@
  * allocators in gpu_config.cpp, and use NUMA strategies.
  * 
  * @author Tyson Jones
+ * @author Mai Đức Khang (CPU memory query)
  */
 
 #include "quest/include/types.h"
diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp
@@ -913,7 +913,10 @@ qreal util_getPhaseFromGateAngle(qreal angle) {
     return - angle / 2;
 }
 
+qcomp util_getPhaseFromGateAngle(qcomp angle) {
 
+    return -angle / qcomp(2.0, 0.0);
+}
 
 /*
  * DECOHERENCE FACTORS
diff --git a/quest/src/core/utilities.hpp b/quest/src/core/utilities.hpp
@@ -350,6 +350,7 @@ util_VectorIndexRange util_getLocalIndRangeOfVectorElemsWithinNode(int rank, qin
 
 qreal util_getPhaseFromGateAngle(qreal angle);
 
+qcomp util_getPhaseFromGateAngle(qcomp angle);
 
 
 /*
diff --git a/tests/unit/operations.cpp b/tests/unit/operations.cpp
@@ -1920,6 +1920,35 @@ TEST_CASE( "multiplyPauliStrSum", TEST_CATEGORY LABEL_MIXED_DEPLOY_TAG ) {
 }
 
 
+TEST_CASE( "applyNonUnitaryPauliGadget", TEST_CATEGORY ) {
+
+    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
+
+    SECTION( LABEL_CORRECTNESS ) {
+
+        // prepare a random Pauli string and angle
+        int numTargs = GENERATE_COPY( range(1, numQubits+1) );
+        auto targs = GENERATE_TARGS( numQubits, numTargs );
+        PauliStr str = getRandomPauliStr(targs);
+        qcomp angle = getRandomComplex();
+
+        // prepare the corresponding reference matrix exp(-i angle pauli)
+        auto matrRef = getExponentialOfPauliMatrix(angle, getMatrix(str, numQubits));
+
+        auto testFunc = [&](Qureg qureg, auto& stateRef) {
+            applyNonUnitaryPauliGadget(qureg, str, angle);
+            applyReferenceOperator(stateRef, matrRef);
+        };
+
+        CAPTURE( targs, angle );
+        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
+        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
+    }
+
+    /// @todo input validation
+}
+
+
 /** @} (end defgroup) */
 
 
diff --git a/tests/utils/linalg.cpp b/tests/utils/linalg.cpp
@@ -345,7 +345,7 @@ qmatrix getExponentialOfDiagonalMatrix(qmatrix m) {
 }
 
 
-qmatrix getExponentialOfPauliMatrix(qreal arg, qmatrix m) {
+qmatrix getExponentialOfPauliMatrix(qcomp arg, qmatrix m) {
     
     // exp(-i arg/2 m) where m = prod(paulis)
     qmatrix id = getIdentityMatrix(m.size());
diff --git a/tests/utils/linalg.hpp b/tests/utils/linalg.hpp
@@ -49,7 +49,7 @@ qmatrix getConjugate(qmatrix);
 qmatrix getConjugateTranspose(qmatrix);
 qmatrix getPowerOfDiagonalMatrix(qmatrix diag, qcomp power);
 qmatrix getExponentialOfDiagonalMatrix(qmatrix);
-qmatrix getExponentialOfPauliMatrix(qreal arg, qmatrix pauli);
+qmatrix getExponentialOfPauliMatrix(qcomp arg, qmatrix pauli);
 qmatrix getExponentialOfNormalisedPauliVector(qreal arg, qreal x, qreal y, qreal z);
 qmatrix getOrthonormalisedRows(qmatrix);
 qmatrix getOrthonormalisedRows(qmatrix);

Original file line number	Diff line number	Diff line change
`@@ -5,6 +5,7 @@`
`5`	`5`	`* instead exposed in decoherence.h`
`6`	`6`	`*`
`7`	`7`	`* @author Tyson Jones`
	`8`	`+ * @author Diogo Pratas Maia (non-unitary Pauli gadget)`
`8`	`9`	`*`
`9`	`10`	`* @defgroup operations Operations`
`10`	`11`	`* @ingroup api`
`@@ -1853,6 +1854,8 @@ void multiplyPauliGadget(Qureg qureg, PauliStr str, qreal angle);`
`1853`	`1854`	`*/`
`1854`	`1855`	`void applyPauliGadget(Qureg qureg, PauliStr str, qreal angle);`
`1855`	`1856`
	`1857`	`+/// @notyetdoced`
	`1858`	`+void applyNonUnitaryPauliGadget(Qureg qureg, PauliStr str, qcomp angle);`
`1856`	`1859`
`1857`	`1860`	`/// @notyetdoced`
`1858`	`1861`	`void applyControlledPauliGadget(Qureg qureg, int control, PauliStr str, qreal angle);`
Original file line number	Diff line number	Diff line change
`@@ -10,6 +10,7 @@`
`10`	`10`	`* allocators in gpu_config.cpp, and use NUMA strategies.`
`11`	`11`	`*`
`12`	`12`	`* @author Tyson Jones`
	`13`	`+ * @author Mai Đức Khang (CPU memory query)`
`13`	`14`	`*/`
`14`	`15`
`15`	`16`	`#include "quest/include/types.h"`
Original file line number	Diff line number	Diff line change
`@@ -913,7 +913,10 @@ qreal util_getPhaseFromGateAngle(qreal angle) {`
`913`	`913`	`return - angle / 2;`
`914`	`914`	`}`
`915`	`915`
	`916`	`+qcomp util_getPhaseFromGateAngle(qcomp angle) {`
`916`	`917`
	`918`	`+ return -angle / qcomp(2.0, 0.0);`
	`919`	`+}`
`917`	`920`
`918`	`921`	`/*`
`919`	`922`	`* DECOHERENCE FACTORS`
Original file line number	Diff line number	Diff line change
`@@ -350,6 +350,7 @@ util_VectorIndexRange util_getLocalIndRangeOfVectorElemsWithinNode(int rank, qin`
`350`	`350`
`351`	`351`	`qreal util_getPhaseFromGateAngle(qreal angle);`
`352`	`352`
	`353`	`+qcomp util_getPhaseFromGateAngle(qcomp angle);`
`353`	`354`
`354`	`355`
`355`	`356`	`/*`