v3 Extend Pauli rotations to accept complex angles #594 (unitaryhack2025)

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

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.

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;
+}

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);

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__);

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)) {

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);
 
 
 /*

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"

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

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);
 
 
 /*

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) */
 
 

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());

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);