added multiControlledMultiRotatePauli

TysonRayJones · TysonRayJones · commit 64c73198631d · 2021-07-11T22:58:46.000+01:00
diff --git a/QuEST/include/QuEST.h b/QuEST/include/QuEST.h
@@ -4027,6 +4027,107 @@ void multiRotatePauli(Qureg qureg, int* targetQubits, enum pauliOpType* targetPa
  */
 void multiControlledMultiRotateZ(Qureg qureg, int* controlQubits, int numControls, int* targetQubits, int numTargets, qreal angle);
 
+/** Apply a multi-controlled multi-target multi-Pauli rotation, also known as a 
+ * controlled Pauli gadget.
+ * This is the unitary 
+ * \f[ 
+ *    |1\rangle\langle 1|^{\otimes\, \text{numControls}} \; \otimes \,
+ *     \exp \left( - i \, \frac{\theta}{2} \; \bigotimes_{j}^{\text{numTargets}} \hat{\sigma}_j\right)
+ *     \;\;+\;\; \sum\limits_{k=0}^{2^{\,\text{numControls}} - 2} |k\rangle\langle k| \otimes \text{I}
+ * \f]
+ * where \f$\hat{\sigma}_j\f$ are the Pauli operators (::pauliOpType) in `targetPaulis`, which operate 
+ * upon the corresponding qubits in `targetQubits`.
+ *
+    \f[
+                \begin{tikzpicture}[scale=.5]
+                \node[draw=none] at (-4, 1) {targets};
+                \node[draw=none] at (-4, 5) {controls};
+                
+                \node[draw=none] at (0, 8) {$\vdots$};
+                \draw (0, 7) -- (0, 6);
+                
+                \draw (-2.5, 6) -- (2.5, 6);
+                \draw[fill=black] (0, 6) circle (.2);
+                \draw (0, 6) -- (0, 4);         
+                
+                \draw (-2.5, 4) -- (2.5, 4);
+                \draw[fill=black] (0, 4) circle (.2);
+                \draw(0, 4) -- (0, 3);
+
+                \draw (-2.5,0) -- (-1.5, 0);
+                \draw (1.5, 0) -- (2.5, 0);
+                \draw (-2.5,2) -- (-1.5, 2);
+                \draw (1.5, 2) -- (2.5, 2);
+                \draw (-1.5,-1)--(-1.5,3)--(1.5,3)--(1.5,-1);
+                \node[draw=none] at (0, 1) {$e^{-i\frac{\theta}{2} \bigotimes\limits_j \hat{\sigma}_j }$};
+                \node[draw=none] at (0, -1) {$\vdots$};
+                
+                \end{tikzpicture}
+    \f]
+ * 
+ * > All qubits not appearing in \p targetQubits and \p controlQubits are assumed to receive the identity operator.
+ *
+ * For example:
+ * ```
+ *     int numCtrls = 1;
+ *     int numTargs = 4;
+ *     int ctrls[] = {4};
+ *     int targs[] = {0,1,2,3};
+ *     
+ *     pauliOpType paulis[] = {PAULI_X, PAULI_Y, PAULI_Z, PAULI_I};
+ *     
+ *     multiControlledMultiRotatePauli(
+ *         qureg, ctrls, numCtrls, targs, paulis, numTargs, 0.1);
+ * ```
+ * effects 
+ * \f[
+ *    |1\rangle\langle 1 | \otimes \exp\left( -i \, (0.1/2) \, X_0 \, Y_1 \, Z_2 \right) \, \text{I}_3
+ *    \;\; + \;\; |0\rangle\langle 0| \otimes \text{I}^{\otimes 4}
+ * \f] 
+ * on \p qureg, where unspecified qubits (along with those targeted by `PAULI_I`) are 
+ * assumed to receive the identity operator (excluded from exponentiation). 
+ *
+ * > This means specifying `PAULI_I` does *not* induce a global phase factor \f$\exp(-i \theta/2)\f$.
+ * > Hence, if all \p targetPaulis are identity, then this function does nothing to \p qureg.
+ * > Specifying `PAULI_I` on a qubit is superfluous but allowed for convenience.
+ *
+ * This function effects the controlled Pauli gadget by first (controlled) 
+ * rotating the qubits which are targeted with either `X` or `Y` into alternate basis, 
+ * performing multiControlledMultiRotateZ() on all target qubits, then restoring 
+ * the original basis. 
+ *
+ * @see
+ * - multiControlledMultiRotateZ()
+ * - multiRotatePauli()
+ * - multiRotateZ()
+ * - rotateX()
+ * - rotateY()
+ * - rotateZ()
+ * - rotateAroundAxis()
+ *
+ * @ingroup unitary
+ * @param[in,out] qureg object representing the set of all qubits
+ * @param[in] controlQubits list of the indices of qubits to control upon
+ * @param[in] numControls length of length `controlQubits`
+ * @param[in] targetQubits a list of the indices of the target qubits 
+ * @param[in] targetPaulis a list of the Pauli operators around which to rotate the target qubits
+ * @param[in] numTargets length of list `targetQubits`
+ * @param[in] angle the angle by which the multi-qubit state is rotated around the Z axis
+ * @throws invalidQuESTInputError()
+ * - if any qubit in \p controlQubits and \p targetQubits is invalid, i.e. outside <b>[0, </b>`qureg.numQubitsRepresented`<b>)</b>
+ * - if \p controlQubits or \p targetQubits contain any repetitions
+ * - if any qubit in \p controlQubits is also in \p targetQubits (and vice versa)
+ * - if \p numTargets <b>< 1</b>
+ * - if \p numControls <b>< 1</b> (use multiRotateZ() for no controls)
+ * - if any element of \p targetPaulis is not one of `PAULI_I`, `PAULI_X`, `PAULI_Y`, `PAULI_Z`
+ * @throws segmentation-fault
+ * - if \p controlQubits contains fewer elements than \p numControls
+ * - if \p targetQubits contains fewer elements than \p numTargets
+ * - if \p targetPaulis contains fewer elements than \p numTargets
+ * @author Tyson Jones
+ */
+void multiControlledMultiRotatePauli(Qureg qureg, int* controlQubits, int numControls, int* targetQubits, enum pauliOpType* targetPaulis, int numTargets, qreal angle);
+
 /** Computes the expected value of a product of Pauli operators.
  * Letting \f$ \sigma = \otimes_j \hat{\sigma}_j \f$ be the operators indicated by \p pauliCodes 
  * and acting on qubits \p targetQubits, this function computes \f$ \langle \psi | \sigma | \psi \rangle \f$ 
diff --git a/QuEST/src/QuEST.c b/QuEST/src/QuEST.c
@@ -703,6 +703,27 @@ void multiRotatePauli(Qureg qureg, int* targetQubits, enum pauliOpType* targetPa
         numTargets, angle);
 }
 
+void multiControlledMultiRotatePauli(Qureg qureg, int* controlQubits, int numControls, int* targetQubits, enum pauliOpType* targetPaulis, int numTargets, qreal angle) {
+    validateMultiControlsMultiTargets(qureg, controlQubits, numControls, targetQubits, numTargets, __func__);
+    validatePauliCodes(targetPaulis, numTargets, __func__);
+    
+    int conj=0;
+    long long int ctrlMask = getQubitBitMask(controlQubits, numControls);
+    statevec_multiControlledMultiRotatePauli(qureg, ctrlMask, targetQubits, targetPaulis, numTargets, angle, conj);
+    if (qureg.isDensityMatrix) {
+        conj = 1;
+        int shift = qureg.numQubitsRepresented;
+        shiftIndices(targetQubits, numTargets, shift);
+        statevec_multiControlledMultiRotatePauli(qureg, ctrlMask<<shift, targetQubits, targetPaulis, numTargets, angle, conj);
+        shiftIndices(targetQubits, numTargets, -shift);
+    }
+    
+    // @TODO: create actual QASM
+    qasm_recordComment(qureg, 
+        "Here a %d-control %d-target multiControlledMultiRotatePauli of angle %g was performed (QASM not yet implemented)",
+        numControls, numTargets, angle);
+}
+
 void applyPhaseFunc(Qureg qureg, int* qubits, int numQubits, enum bitEncoding encoding, qreal* coeffs, qreal* exponents, int numTerms) {
     validateStateVecQureg(qureg, __func__);
     validateMultiQubits(qureg, qubits, numQubits, __func__);
diff --git a/QuEST/src/QuEST_common.c b/QuEST/src/QuEST_common.c
@@ -447,6 +447,47 @@ void statevec_multiRotatePauli(
     }
 }
 
+void statevec_multiControlledMultiRotatePauli(
+    Qureg qureg, long long int ctrlMask, int* targetQubits, enum pauliOpType* targetPaulis, int numTargets, qreal angle,
+    int applyConj
+) {
+    qreal fac = 1/sqrt(2);
+    qreal sgn = (applyConj)? 1 : -1;
+    ComplexMatrix2 uRx = {.real={{fac,0},{0,fac}}, .imag={{0,sgn*fac},{sgn*fac,0}}}; // Rx(pi/2)* rotates Z -> Y
+    ComplexMatrix2 uRy = {.real={{fac,fac},{-fac,fac}}, .imag={{0,0},{0,0}}};        // Ry(-pi/2) rotates Z -> X
+    
+    // this function is controlled on the all-one state, so no ctrl flips
+    long long int ctrlFlipMask = 0;
+    
+    // mask may be modified to remove superfluous Identity ops
+    long long int targMask = getQubitBitMask(targetQubits, numTargets);
+    
+    // rotate basis so that exp(Z) will effect exp(Y) and exp(X)
+    for (int t=0; t < numTargets; t++) {
+        if (targetPaulis[t] == PAULI_I)
+            targMask -= 1LL << targetQubits[t]; // remove target from mask
+        if (targetPaulis[t] == PAULI_X)
+            statevec_multiControlledUnitary(qureg, ctrlMask, ctrlFlipMask, targetQubits[t], uRy);
+        if (targetPaulis[t] == PAULI_Y)
+            statevec_multiControlledUnitary(qureg, ctrlMask, ctrlFlipMask, targetQubits[t], uRx);
+        // (targetPaulis[t] == 3) is Z basis
+    }
+    
+    // does nothing if there are no qubits to 'rotate'
+    if (targMask != 0)
+        statevec_multiControlledMultiRotateZ(qureg, ctrlMask, targMask, (applyConj)? -angle : angle);
+        
+    // undo X and Y basis rotations
+    uRx.imag[0][1] *= -1; uRx.imag[1][0] *= -1;
+    uRy.real[0][1] *= -1; uRy.real[1][0] *= -1;
+    for (int t=0; t < numTargets; t++) {
+        if (targetPaulis[t] == PAULI_X)
+            statevec_multiControlledUnitary(qureg, ctrlMask, ctrlFlipMask, targetQubits[t], uRy);
+        if (targetPaulis[t] == PAULI_Y)
+            statevec_multiControlledUnitary(qureg, ctrlMask, ctrlFlipMask, targetQubits[t], uRx);
+    }
+}
+
 /* produces both pauli|qureg> or pauli * qureg (as a density matrix) */
 void statevec_applyPauliProd(Qureg workspace, int* targetQubits, enum pauliOpType* pauliCodes, int numTargets) {
     
diff --git a/QuEST/src/QuEST_internal.h b/QuEST/src/QuEST_internal.h
@@ -253,6 +253,8 @@ void statevec_multiControlledMultiRotateZ(Qureg qureg, long long int ctrlMask, l
 
 void statevec_multiRotatePauli(Qureg qureg, int* targetQubits, enum pauliOpType* targetPaulis, int numTargets, qreal angle, int applyConj);
 
+void statevec_multiControlledMultiRotatePauli(Qureg qureg, long long int ctrlMask, int* targetQubits, enum pauliOpType* targetPaulis, int numTargets, qreal angle, int applyConj);
+
 void statevec_setWeightedQureg(Complex fac1, Qureg qureg1, Complex fac2, Qureg qureg2, Complex facOut, Qureg out);
 
 void statevec_applyPauliSum(Qureg inQureg, enum pauliOpType* allCodes, qreal* termCoeffs, int numSumTerms, Qureg outQureg);
diff --git a/tests/test_unitaries.cpp b/tests/test_unitaries.cpp
@@ -1115,6 +1115,169 @@ TEST_CASE( "multiControlledMultiQubitUnitary", "[unitaries]" ) {
 
 
 
+/** @sa multiControlledMultiRotatePauli
+ * @ingroup unittest 
+ * @author Tyson Jones 
+ */
+TEST_CASE( "multiControlledMultiRotatePauli", "[unitaries]" ) {
+        
+    PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
+    qreal param = getRandomReal(-4*M_PI, 4*M_PI);
+        
+    SECTION( "correctness" ) {
+        
+        // try all possible numbers of targets and controls
+        int numTargs = GENERATE_COPY( range(1,NUM_QUBITS) ); // leave space for min 1 control qubit
+        int maxNumCtrls = NUM_QUBITS - numTargs;
+        int numCtrls = GENERATE_COPY( range(1,maxNumCtrls+1) );
+        
+        // generate all possible valid qubit arrangements
+        int* targs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numTargs) );
+        int* ctrls = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numCtrls, targs, numTargs) );
+        
+        /* it's too expensive to try ALL Pauli sequences, via 
+         *      pauliOpType* paulis = GENERATE_COPY( pauliseqs(numTargs) );.
+         * Furthermore, take(10, pauliseqs(numTargs)) will try the same pauli codes.
+         * Hence, we instead opt to randomly generate pauliseqs
+         */
+        pauliOpType paulis[numTargs];
+        for (int i=0; i<numTargs; i++)
+            paulis[i] = (pauliOpType) getRandomInt(0,4);
+
+        // exclude identities from reference matrix exp (they apply unwanted global phase)
+        int refTargs[numTargs];
+        int numRefTargs = 0;
+
+        QMatrix xMatr{{0,1},{1,0}};
+        QMatrix yMatr{{0,-qcomp(0,1)},{qcomp(0,1),0}};
+        QMatrix zMatr{{1,0},{0,-1}};
+        
+        // build correct reference matrix by pauli-matrix exponentiation...
+        QMatrix pauliProd{{1}};
+        for (int i=0; i<numTargs; i++) {
+            QMatrix fac;
+            if (paulis[i] == PAULI_I) continue; // exclude I-targets from ref list
+            if (paulis[i] == PAULI_X) fac = xMatr;
+            if (paulis[i] == PAULI_Y) fac = yMatr;
+            if (paulis[i] == PAULI_Z) fac = zMatr;
+            pauliProd = getKroneckerProduct(fac, pauliProd);
+            
+            // include this target in ref list
+            refTargs[numRefTargs++] = targs[i];
+        }
+
+        // produces exp(-i param/2 pauliProd), unless pauliProd = I
+        QMatrix op;
+        if (numRefTargs > 0)
+            op = getExponentialOfPauliMatrix(param, pauliProd);
+        
+        SECTION( "state-vector" ) {
+            
+            multiControlledMultiRotatePauli(quregVec, ctrls, numCtrls, targs, paulis, numTargs, param);
+            if (numRefTargs > 0)
+                applyReferenceOp(refVec, ctrls, numCtrls, refTargs, numRefTargs, op);
+            REQUIRE( areEqual(quregVec, refVec) );
+        }
+        SECTION( "density-matrix" ) {
+            
+            multiControlledMultiRotatePauli(quregMatr, ctrls, numCtrls, targs, paulis, numTargs, param);
+            if (numRefTargs > 0)
+                applyReferenceOp(refMatr, ctrls, numCtrls, refTargs, numRefTargs, op);
+            REQUIRE( areEqual(quregMatr, refMatr, 10*REAL_EPS) );
+        }
+    }
+    SECTION( "input validation" ) {
+        
+        // test all validation on both state-vector and density-matrix.
+        // want GENERATE_COPY( quregVec, quregMatr ), but too lazy to patch
+        // using github.com/catchorg/Catch2/issues/1809
+        Qureg regs[] = {quregVec, quregMatr};
+        Qureg qureg = regs[GENERATE(0,1)];
+        
+        // over-sized array to prevent seg-fault in case of validation fail below
+        pauliOpType paulis[NUM_QUBITS+1];
+        for (int q=0; q<NUM_QUBITS+1; q++)
+            paulis[q] = PAULI_I;
+
+        SECTION( "pauli codes" ) {
+            
+            int numCtrls = 1;
+            int ctrls[] = {3};
+            int numTargs = 3;
+            int targs[3] = {0, 1, 2};
+
+            // make a single Pauli invalid
+            paulis[GENERATE_COPY(range(0,numTargs))] = (pauliOpType) GENERATE( -1, 4 );
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, ctrls, numCtrls, targs, paulis, numTargs, param), Contains("Invalid Pauli code"));
+        }
+        SECTION( "number of targets" ) {
+            
+            // there cannot be more targets than qubits in register
+            // (numTargs=NUM_QUBITS is caught elsewhere, because that implies ctrls are invalid)
+            int numTargs = GENERATE( -1, 0, NUM_QUBITS+1 );
+            int numCtrls = 1;
+            int targs[NUM_QUBITS+1]; // prevents seg-fault if validation doesn't trigger
+            int ctrls[] = {0};
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, ctrls, numCtrls, targs, paulis, numTargs, param), Contains("Invalid number of target") );
+        }
+        SECTION( "repetition in targets" ) {
+            
+            int numCtrls = 1;
+            int numTargs = 3;
+            int ctrls[] = {0};
+            int targs[] = {1,2,2};
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, ctrls, numCtrls, targs, paulis, numTargs, param), Contains("target") && Contains("unique"));
+        }
+        SECTION( "number of controls" ) {
+            
+            int numCtrls = GENERATE( -1, 0, NUM_QUBITS, NUM_QUBITS+1 );
+            int numTargs = 1;
+            int ctrls[NUM_QUBITS+1]; // avoids seg-fault if validation not triggered
+            int targs[1] = {0};
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, ctrls, numCtrls, targs, paulis, numTargs, param), Contains("Invalid number of control"));
+        }
+        SECTION( "repetition in controls" ) {
+            
+            int numCtrls = 3;
+            int numTargs = 1;
+            int ctrls[] = {0,1,1};
+            int targs[] = {3};
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, ctrls, numCtrls, targs, paulis, numTargs, param), Contains("control") && Contains("unique"));
+        }
+        SECTION( "control and target collision" ) {
+            
+            int numCtrls = 3;
+            int numTargs = 3;
+            int ctrls[] = {0,1,2};
+            int targs[] = {3,1,4};
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, ctrls, numCtrls, targs, paulis, numTargs, param), Contains("Control") && Contains("target") && Contains("disjoint"));
+        }
+        SECTION( "qubit indices" ) {
+            
+            // valid inds
+            int numQb = 2;
+            int qb1[2] = {0,1};
+            int qb2[2] = {2,3};
+            
+            // make qb1 invalid
+            int inv = GENERATE( -1, NUM_QUBITS );
+            qb1[GENERATE_COPY(range(0,numQb))] = inv;
+            
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, qb1, numQb, qb2, paulis, numQb, param), Contains("Invalid control") );
+            REQUIRE_THROWS_WITH( multiControlledMultiRotatePauli(qureg, qb2, numQb, qb1, paulis, numQb, param),  Contains("Invalid target") );
+        }
+    }
+    CLEANUP_TEST( quregVec, quregMatr );
+}
+
+
+
 /** @sa multiControlledMultiRotateZ
  * @ingroup unittest 
  * @author Tyson Jones 
