added controlled Trotter

New API functions:
- applyControlledTrotterizedPauliStrSumGadget()
- applyMultiControlledTrotterizedPauliStrSumGadget()
- applyMultiStateControlledTrotterizedPauliStrSumGadget()
- C++-only std::vector overloads of the latter two.

Additionally:
- renamed the internal constituent functions, like applyFirstOrderTrotter(), to explicit internal_applyFirstOrderTrotterRepetition()
- renamed paulis_getInds() to paulis_getTargetInds()

Note that new validation was required to check that no PauliStrSum non-identity Paulis overlapped the control qubits. This is relatively expensive; we build a PauliStrSum target-mask in O(#terms * #qubits) time whereas the previous most expensive validation (checking PauliStrSum targets do not exceed Qureg) costs O(#terms * log(#qubits)). Such costs are still completely occluded by those of simulating/processing a PauliStrSum in the backend, but might still attract lazy evaluation of the target-mask which is bound to the PauliStrSum instance. We have deferred any such optimisation and the associated struct changes since it necessitates an update to the PauliStrSum design (like new sync functions)

Author: TysonRayJones

quest/include/operations.h

@@ -2387,6 +2387,30 @@ void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace);
 void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps);
 
 
+/// @notyetdoced
+/// @notyettested
+/// @see
+///  - applyTrotterizedPauliStrSumGadget()
+///  - applyControlledCompMatr1()
+void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps);
+
+
+/// @notyetdoced
+/// @notyettested
+/// @see
+///  - applyTrotterizedPauliStrSumGadget()
+///  - applyMultiControlledCompMatr1()
+void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps);
+
+
+/// @notyetdoced
+/// @notyettested
+/// @see
+///  - applyTrotterizedPauliStrSumGadget()
+///  - applyMultiStateControlledCompMatr1()
+void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps);
+
+
 /** @notyettested
  * 
  * A generalisation of applyTrotterizedPauliStrSumGadget() which accepts a complex angle and permits
@@ -2492,6 +2516,26 @@ void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, q
 }
 #endif
 
+#ifdef __cplusplus
+
+
+/// @notyettested
+/// @notyetvalidated
+/// @notyetdoced
+/// @cppvectoroverload
+/// @see applyMultiControlledTrotterizedPauliStrSumGadget()
+void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vector<int> controls, PauliStrSum sum, qreal angle, int order, int reps);
+
+
+/// @notyettested
+/// @notyetvalidated
+/// @notyetdoced
+/// @cppvectoroverload
+/// @see applyMultiStateControlledTrotterizedPauliStrSumGadget()
+void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vector<int> controls, std::vector<int> states, PauliStrSum sum, qreal angle, int order, int reps);
+
+
+#endif // __cplusplus
 
 /** @} */
 

quest/src/api/operations.cpp

@@ -1130,90 +1130,145 @@ void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace) {
     // workspace -> qureg, and qureg -> sum * qureg
 }
 
-void applyFirstOrderTrotter(Qureg qureg, PauliStrSum sum, qcomp angle, bool reverse) {
-
-    // (internal, invoked by applyTrotterizedPauliStrSumGadget)
-
+void internal_applyFirstOrderTrotterRepetition(
+    Qureg qureg, vector<int>& ketCtrls, vector<int>& braCtrls,
+    vector<int>& states, PauliStrSum sum, qcomp angle, bool reverse
+) {
+    // apply each sum term as a gadget, in forward or reverse order
     for (qindex i=0; i<sum.numTerms; i++) {
         int j = reverse? sum.numTerms - i - 1 : i;
+        qcomp coeff = sum.coeffs[j];
+        PauliStr str = sum.strings[j];
 
-        // effect exp(i angle * sum) by undoing gadget pre-factor
-        qcomp arg = angle * sum.coeffs[j] / util_getPhaseFromGateAngle(1);
-        applyNonUnitaryPauliGadget(qureg, sum.strings[j], arg); // caller disabled valiation therein
-    }
-}
+        // effect |psi> -> exp(i angle * sum)|psi>
+        qcomp arg = angle * coeff;
+        localiser_statevec_anyCtrlPauliGadget(qureg, ketCtrls, states, str, arg);
 
-void applyHigherOrderTrotter(Qureg qureg, PauliStrSum sum, qcomp angle, int order) {
+        if (!qureg.isDensityMatrix)
+            continue;
 
-    // (internal, invoked by applyTrotterizedPauliStrSumGadget)
+        // effect rho -> rho dagger(i angle * sum)
+        arg *= paulis_hasOddNumY(str) ? 1 : -1;
+        str = paulis_getShiftedPauliStr(str, qureg.numQubits);
+        localiser_statevec_anyCtrlPauliGadget(qureg, braCtrls, states, str, arg);
+    }
+}
 
+void internal_applyHigherOrderTrotterRepetition(
+    Qureg qureg, vector<int>& ketCtrls, vector<int>& braCtrls,
+    vector<int>& states, PauliStrSum sum, qcomp angle, int order
+) {
     if (order == 1) {
-        applyFirstOrderTrotter(qureg, sum, angle, false);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, false);
 
     } else if (order == 2) {
-        applyFirstOrderTrotter(qureg, sum, angle/2, false);
-        applyFirstOrderTrotter(qureg, sum, angle/2, true);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, false);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, true);
 
     } else {
         qreal p = 1. / (4 - std::pow(4, 1./(order-1)));
         qcomp a = p * angle;
         qcomp b = (1-4*p) * angle;
 
         int lower = order - 2;
-        applyFirstOrderTrotter(qureg, sum, a, lower);
-        applyFirstOrderTrotter(qureg, sum, a, lower);
-        applyFirstOrderTrotter(qureg, sum, b, lower);
-        applyFirstOrderTrotter(qureg, sum, a, lower);
-        applyFirstOrderTrotter(qureg, sum, a, lower);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower);
+        internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower);
     }
 }
 
-void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps) {
-    validate_quregFields(qureg, __func__);
-    validate_pauliStrSumFields(sum, __func__);
-    validate_pauliStrSumTargets(sum, qureg, __func__);
-    validate_trotterParams(qureg, order, reps, __func__);
-    // sum is permitted to be non-Hermitian
-
+void internal_applyAllTrotterRepetitions(
+    Qureg qureg, int* controls, int* states, int numControls, 
+    PauliStrSum sum, qcomp angle, int order, int reps
+) {
     // exp(i angle sum) = identity when angle=0
     if (angle == qcomp(0,0))
         return;
 
-    // record validation state then disable to avoid repeated
-    // re-validations in each invoked applyPauliGadget() below
-    bool wasValidationEnabled = validateconfig_isEnabled();
-    validateconfig_disable();
+    // prepare control-qubit lists once for all invoked gadgets below
+    auto ketCtrlsVec = util_getVector(controls, numControls);
+    auto braCtrlsVec = (qureg.isDensityMatrix)? util_getBraQubits(ketCtrlsVec, qureg) : vector<int>{};
+    auto statesVec = util_getVector(states, numControls);
 
-    // perform sequence of applyPauliGadget()
-    for (int r=0; r<reps; r++)
-        applyHigherOrderTrotter(qureg, sum, angle/reps, order);
+    qcomp arg = angle / reps;
 
-    // potentially restore validation
-    if (wasValidationEnabled)
-        validateconfig_enable();
+    // perform carefully-ordered sequence of gadgets
+    for (int r=0; r<reps; r++)
+        internal_applyHigherOrderTrotterRepetition(
+            qureg, ketCtrlsVec, braCtrlsVec, statesVec, sum, arg, order);
 
     /// @todo
     /// the accuracy of Trotterisation is greatly improved by randomisation
     /// or (even sub-optimal) grouping into commuting terms. Should we 
     /// implement these above or into another function?
 }
 
-void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps) {
+    validate_quregFields(qureg, __func__);
+    validate_pauliStrSumFields(sum, __func__);
+    validate_pauliStrSumTargets(sum, qureg, __func__);
+    validate_trotterParams(qureg, order, reps, __func__);
+    // sum is permitted to be non-Hermitian
+
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps);
+}
 
-    // validate inputs here despite re-validation below so that func name is correct in error message
+void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumTargets(sum, qureg, __func__);
     validate_trotterParams(qureg, order, reps, __func__);
+    validate_pauliStrSumIsHermitian(sum, __func__);
+
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps);
+}
+
+void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps) {
+    validate_quregFields(qureg, __func__);
+    validate_pauliStrSumFields(sum, __func__);
+    validate_controlAndPauliStrSumTargets(qureg, control, sum, __func__);
+    validate_trotterParams(qureg, order, reps, __func__);
+    validate_pauliStrSumIsHermitian(sum, __func__);
 
-    // crucially, require that sum coefficients are real
+    internal_applyAllTrotterRepetitions(qureg, &control, nullptr, 1, sum, angle, order, reps);
+}
+
+void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps) {
+    validate_quregFields(qureg, __func__);
+    validate_pauliStrSumFields(sum, __func__);
+    validate_controlsAndPauliStrSumTargets(qureg, controls, numControls, sum, __func__);
+    validate_trotterParams(qureg, order, reps, __func__);
     validate_pauliStrSumIsHermitian(sum, __func__);
 
-    applyNonUnitaryTrotterizedPauliStrSumGadget(qureg, sum, angle, order, reps);
+    internal_applyAllTrotterRepetitions(qureg, controls, nullptr, numControls, sum, angle, order, reps);
+}
+
+void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps) {
+    validate_quregFields(qureg, __func__);
+    validate_pauliStrSumFields(sum, __func__);
+    validate_controlsAndPauliStrSumTargets(qureg, controls, numControls, sum, __func__);
+    validate_controlStates(states, numControls, __func__); // permits states==nullptr
+    validate_trotterParams(qureg, order, reps, __func__);
+    validate_pauliStrSumIsHermitian(sum, __func__);
+
+    internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps);
 }
 
 } // end de-mangler
 
+void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, vector<int> controls, PauliStrSum sum, qreal angle, int order, int reps) {
+
+    applyMultiControlledTrotterizedPauliStrSumGadget(qureg, controls.data(), controls.size(), sum, angle, order, reps);
+}
+
+void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, vector<int> controls, vector<int> states, PauliStrSum sum, qreal angle, int order, int reps) {
+    validate_controlsMatchStates(controls.size(), states.size(), __func__);
+
+    applyMultiStateControlledTrotterizedPauliStrSumGadget(qureg, controls.data(), states.data(), controls.size(), sum, angle, order, reps);
+}
+
 
 
 /*

quest/src/api/paulis.cpp

@@ -207,24 +207,41 @@ qcomp paulis_getPrefixPaulisElem(Qureg qureg, vector<int> prefixY, vector<int> p
 }
 
 
-vector<int> paulis_getInds(PauliStr str) {
+vector<int> paulis_getTargetInds(PauliStr str) {
 
     int maxInd = paulis_getIndOfLefmostNonIdentityPauli(str);
 
     vector<int> inds(0);
     inds.reserve(maxInd+1);
 
     for (int i=0; i<=maxInd; i++)
-        if (paulis_getPauliAt(str, i) != 0)
+        if (paulis_getPauliAt(str, i) != 0) // Id
             inds.push_back(i);
 
     return inds;
 }
 
 
+qindex paulis_getTargetBitMask(PauliStr str) {
+    
+    /// @todo 
+    /// would compile-time MAX_NUM_PAULIS_PER_STR bound be faster here,
+    /// since this function is invoked upon every PauliStrSum element?
+    int maxInd = paulis_getIndOfLefmostNonIdentityPauli(str);
+
+    qindex mask = 0;
+
+    for (int i=0; i<=maxInd; i++)
+        if (paulis_getPauliAt(str, i) != 0) // Id
+            mask = flipBit(mask, i);
+
+    return mask;
+}
+
+
 array<vector<int>,3> paulis_getSeparateInds(PauliStr str, Qureg qureg) {
 
-    vector<int> iXYZ = paulis_getInds(str);
+    vector<int> iXYZ = paulis_getTargetInds(str);
     vector<int> iX, iY, iZ;
 
     vector<int>* ptrs[] = {&iX, &iY, &iZ};
@@ -295,6 +312,18 @@ PAULI_MASK_TYPE paulis_getKeyOfSameMixedAmpsGroup(PauliStr str) {
 }
 
 
+qindex paulis_getTargetBitMask(PauliStrSum sum) {
+
+    qindex mask = 0;
+
+    // mask has 1 where any str has a != Id
+    for (int t=0; t<sum.numTerms; t++)
+        mask |= paulis_getTargetBitMask(sum.strings[t]);
+
+    return mask;
+}
+
+
 
 /*
  * PAULI STRING INITIALISATION

quest/src/core/localiser.cpp

@@ -1238,7 +1238,7 @@ template void localiser_statevec_anyCtrlAnyTargAnyMatr(Qureg, vector<int>, vecto
 
 
 extern bool paulis_containsXOrY(PauliStr str);
-extern vector<int> paulis_getInds(PauliStr str);
+extern vector<int> paulis_getTargetInds(PauliStr str);
 extern std::array<vector<int>,3> paulis_getSeparateInds(PauliStr str, Qureg qureg);
 extern int paulis_getPrefixZSign(Qureg qureg, vector<int> prefixZ) ;
 extern qcomp paulis_getPrefixPaulisElem(Qureg qureg, vector<int> prefixY, vector<int> prefixZ);
@@ -1334,7 +1334,7 @@ void localiser_statevec_anyCtrlPauliTensor(Qureg qureg, vector<int> ctrls, vecto
 
         bool isGadget = false;
         qreal phase = 0; // ignored
-        anyCtrlZTensorOrGadget(qureg, ctrls, ctrlStates, paulis_getInds(str), isGadget, phase);
+        anyCtrlZTensorOrGadget(qureg, ctrls, ctrlStates, paulis_getTargetInds(str), isGadget, phase);
     }
 }
 
@@ -1350,7 +1350,7 @@ void localiser_statevec_anyCtrlPauliGadget(Qureg qureg, vector<int> ctrls, vecto
 
     // when str=IZ, we must use the above bespoke algorithm
     if (!paulis_containsXOrY(str)) {
-        localiser_statevec_anyCtrlPhaseGadget(qureg, ctrls, ctrlStates, paulis_getInds(str), phase);
+        localiser_statevec_anyCtrlPhaseGadget(qureg, ctrls, ctrlStates, paulis_getTargetInds(str), phase);
         return;
     }