Add inverse to QFT implementation

README.md

@@ -94,7 +94,7 @@ initRandomPureState(qureg);
 
 applyHadamard(qureg, 0);
 applyControlledRotateX(qureg, 0, 1, angle);
-applyFullQuantumFourierTransform(qureg);
+applyFullQuantumFourierTransform(qureg, inverse);
 
 reportQureg(qureg);
 

docs/tutorial.md

@@ -626,7 +626,8 @@ applyCompMatr1(qureg, 0, m);
 QuEST includes a few convenience functions for effecting [QFT](https://quest-kit.github.io/QuEST/group__op__qft.html) and [Trotter](https://quest-kit.github.io/QuEST/group__op__paulistrsum.html) circuits.
 
 ```cpp
-applyQuantumFourierTransform(qureg, targets, 3);
+bool inverse = false;
+applyQuantumFourierTransform(qureg, targets, 3, inverse);
 
 qreal time = .3;
 int order = 4;
@@ -865,4 +866,4 @@ This is important because it ensures:
 > [!CAUTION]
 > After calling `finalizeQuESTEnv()`, MPI will close and if being accessed directly by the user, will enter an undefined state. Subsequent calls to MPI routines may return gibberish, and distributed machines will lose their ability to communicate. It is recommended to call `finalizeQuESTEnv()` immediately before exiting.
 
-You are now a QuEST expert 🎉 though there are _many_ more functions in the [API](https://quest-kit.github.io/QuEST/group__api.html) not covered here. Go forth and simulate!
+You are now a QuEST expert 🎉 though there are _many_ more functions in the [API](https://quest-kit.github.io/QuEST/group__api.html) not covered here. Go forth and simulate!

quest/include/operations.h

@@ -25,6 +25,8 @@
 #include "quest/include/matrices.h"
 #include "quest/include/channels.h"
 
+#include <stdbool.h>
+
 #ifdef __cplusplus
     #include <vector>
 #endif
@@ -2406,14 +2408,70 @@ extern "C" {
 #endif
 
 
-/// @notyetdoced
-/// @notyetvalidated
-void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets);
+/** @notyetdoced
+ * @notyetvalidated
+ *
+ * Applies the Quantum Fourier Transform to the specified @p targets of @p qureg.
+ * Alternatively, applies the Inverse Quantum Fourier Transform according to @p inverse.
+ * 
+ * @formulae
+ * 
+ * Let @f$ N @f$ = @p numTargets, the @f$ N @f$ qubit Quantum Fourier Transform maps each
+ * computational basis state belonging to the targeted qubits, @f$ \ket{j} @f$, according to
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ * Similarly the Inverse Quantum Fourier Transform maps each basis state like
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{-2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ *
+ * @param[in,out] qureg      the state to modify.
+ * @param[in]     targets    the indices of the target qubits.
+ * @param[in]     numTargets the length of list @p targets
+ * @param[in]     inverse    whether to apply the inverse QFT or forward QFT
+ * @throws @validationerror
+ * - if @p qureg is uninitialised.
+*  - if @p targets are invalid qubit indices.
+ * - if @p numTargets < 1.
+ * @see
+ * - applyFullQuantumFourierTransform()
+ * @author Vasco Ferreira
+ */
+void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets, bool inverse);
 
 
-/// @notyetdoced
-/// @notyetvalidated
-void applyFullQuantumFourierTransform(Qureg qureg);
+/** @notyetdoced
+ * @notyetvalidated
+ *
+ * Applies the Quantum Fourier Transform to each qubit in @p qureg. Alternatively,
+ * applies the Inverse Quantum Fourier Transform according to @p inverse.
+ * 
+ * @formulae
+ * 
+ * The Quantum Fourier Transform maps each computational basis state @f$ \ket{j} @f$
+ * in an @f$ N @f$ qubit @p qureg according to
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ * Similarly the Inverse Quantum Fourier Transform maps each basis state like
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{-2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ *
+ * @equivalences
+ *
+ * - This function wraps applyQuantumFourierTransform() with all the qubits in the @p qureg as @p targets.
+ *
+ * @param[in,out] qureg      the state to modify.
+ * @param[in]     inverse    whether to apply the inverse QFT or forward QFT
+ * @throws @validationerror
+ * - if @p qureg is uninitialised.
+ * @see
+ * - applyQuantumFourierTransform()
+ * @author Vasco Ferreira
+ */
+void applyFullQuantumFourierTransform(Qureg qureg, bool inverse);
 
 
 // end de-mangler
@@ -2429,7 +2487,7 @@ void applyFullQuantumFourierTransform(Qureg qureg);
 /// @notyetdoced
 /// @cppvectoroverload
 /// @see applyQuantumFourierTransform()
-void applyQuantumFourierTransform(Qureg qureg, std::vector<int> targets);
+void applyQuantumFourierTransform(Qureg qureg, std::vector<int> targets, bool inverse);
 
 
 #endif // __cplusplus

quest/src/api/operations.cpp

@@ -1736,41 +1736,58 @@ qreal applyForcedMultiQubitMeasurement(Qureg qureg, vector<int> qubits, vector<i
 
 extern "C" {
 
-void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets) {
+void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets, bool inverse) {
     validate_quregFields(qureg, __func__);
     validate_targets(qureg, targets, numTargets, __func__);
 
     /// @todo
     /// change this placeholder implementation to the bespoke, optimised routine,
     /// wherein each contiguous controlled-phase gate is merged
 
-    for (int n=numTargets-1; n>=0; n--) {
-        applyHadamard(qureg, targets[n]);
-        for (int m=0; m<n; m++) {
-            qreal arg = const_PI / powerOf2(m+1);
-            applyTwoQubitPhaseShift(qureg, targets[n], targets[n-m-1], arg);
+    int mid = numTargets/2; // floors
+
+    // don't think it is worth overcomplicating this to avoid slight repetition
+    if (inverse) {
+        for (int n=0; n<mid; n++)
+            applySwap(qureg, targets[n], targets[numTargets-1-n]);
+
+        for (int n=0; n<numTargets; n++) {
+            for (int m=0; m<n; m++) {
+                qreal arg = - const_PI / powerOf2(m+1);
+                applyTwoQubitPhaseShift(qureg, targets[n], targets[n-m-1], arg);
+            }
+            applyHadamard(qureg, targets[n]);
         }
-    }
 
-    int mid = numTargets/2; // floors
-    for (int n=0; n<mid; n++)
-        applySwap(qureg, targets[n], targets[numTargets-1-n]);
+    } else {
+        for (int n=numTargets-1; n>=0; n--) {
+            applyHadamard(qureg, targets[n]);
+
+            for (int m=0; m<n; m++) {
+                qreal arg = const_PI / powerOf2(m+1);
+                applyTwoQubitPhaseShift(qureg, targets[n], targets[n-m-1], arg);
+            }
+        }
+
+        for (int n=0; n<mid; n++)
+            applySwap(qureg, targets[n], targets[numTargets-1-n]);
+    }
 }
 
-void applyFullQuantumFourierTransform(Qureg qureg) {
+void applyFullQuantumFourierTransform(Qureg qureg, bool inverse) {
     validate_quregFields(qureg, __func__);
 
     // tiny; no need to validate alloc
     vector<int> targets(qureg.numQubits);
     for (size_t i=0; i<targets.size(); i++)
         targets[i] = i;
 
-    applyQuantumFourierTransform(qureg, targets.data(), targets.size());
+    applyQuantumFourierTransform(qureg, targets.data(), targets.size(), inverse);
 }
 
 } // end de-mangler
 
-void applyQuantumFourierTransform(Qureg qureg, vector<int> targets) {
+void applyQuantumFourierTransform(Qureg qureg, vector<int> targets, bool inverse) {
 
-    applyQuantumFourierTransform(qureg, targets.data(), targets.size());
+    applyQuantumFourierTransform(qureg, targets.data(), targets.size(), inverse);
 }

tests/unit/operations.cpp

@@ -1372,16 +1372,18 @@ TEST_CASE( "applyQuantumFourierTransform", TEST_CATEGORY_OPS ) {
 
         int numTargs = GENERATE_COPY( range(1,numQubits+1) );
         auto targs = GENERATE_TARGS( numQubits, numTargs );
+        bool inverse = GENERATE(false, true);
 
         CAPTURE( targs );
 
         SECTION( LABEL_STATEVEC ) { 
 
             auto testFunc = [&](Qureg qureg, qvector& ref) {
-                applyQuantumFourierTransform(qureg, targs.data(), targs.size());
-                ref = getDisceteFourierTransform(ref, targs);
+                applyQuantumFourierTransform(qureg, targs.data(), targs.size(), inverse);
+                ref = getDiscreteFourierTransform(ref, targs, inverse);
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc);
         }
 
@@ -1395,15 +1397,16 @@ TEST_CASE( "applyQuantumFourierTransform", TEST_CATEGORY_OPS ) {
 
                 // overwrite the Qureg debug state set by caller to above mixture
                 setQuregToReference(qureg, getMixture(states, probs));
-                applyQuantumFourierTransform(qureg, targs.data(), targs.size());
+                applyQuantumFourierTransform(qureg, targs.data(), targs.size(), inverse);
 
                 ref = getZeroMatrix(ref.size());
                 for (size_t i=0; i<states.size(); i++) {
-                    qvector vec = getDisceteFourierTransform(states[i], targs);
+                    qvector vec = getDiscreteFourierTransform(states[i], targs, inverse);
                     ref += probs[i] * getOuterProduct(vec, vec);
                 }
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc);
         }
     }
@@ -1419,14 +1422,16 @@ TEST_CASE( "applyFullQuantumFourierTransform", TEST_CATEGORY_OPS ) {
     SECTION( LABEL_CORRECTNESS ) {
 
         GENERATE( range(0,10) );
+        bool inverse = GENERATE(false, true);
 
         SECTION( LABEL_STATEVEC ) { 
 
             auto testFunc = [&](Qureg qureg, qvector& ref) {
-                applyFullQuantumFourierTransform(qureg);
-                ref = getDisceteFourierTransform(ref);
+                applyFullQuantumFourierTransform(qureg, inverse);
+                ref = getDiscreteFourierTransform(ref, inverse);
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc);
         }
 
@@ -1440,15 +1445,16 @@ TEST_CASE( "applyFullQuantumFourierTransform", TEST_CATEGORY_OPS ) {
 
                 // overwrite the Qureg debug state set by caller to above mixture
                 setQuregToReference(qureg, getMixture(states, probs));
-                applyFullQuantumFourierTransform(qureg);
+                applyFullQuantumFourierTransform(qureg, inverse);
 
                 ref = getZeroMatrix(ref.size());
                 for (size_t i=0; i<states.size(); i++) {
-                    qvector vec = getDisceteFourierTransform(states[i]);
+                    qvector vec = getDiscreteFourierTransform(states[i], inverse);
                     ref += probs[i] * getOuterProduct(vec, vec);
                 }
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc);
         }
     }

tests/utils/linalg.cpp

@@ -157,7 +157,7 @@ qvector getNormalised(qvector vec) {
 }
 
 
-qvector getDisceteFourierTransform(qvector in) {
+qvector getDiscreteFourierTransform(qvector in, bool inverse) {
     DEMAND( in.size() > 0 );
 
     size_t dim = in.size();
@@ -166,7 +166,7 @@ qvector getDisceteFourierTransform(qvector in) {
     // PI must be accurate here
     qreal pi = 3.14159265358979323846;
     qreal a = 1 / std::sqrt(dim);
-    qreal b = 2 * pi / dim;
+    qreal b = (inverse ? -1 : 1) * 2 * pi / dim;
 
     for (size_t x=0; x<dim; x++)
         for (size_t y=0; y<dim; y++)
@@ -176,7 +176,7 @@ qvector getDisceteFourierTransform(qvector in) {
 }
 
 
-qvector getDisceteFourierTransform(qvector in, vector<int> targs) {
+qvector getDiscreteFourierTransform(qvector in, vector<int> targs, bool inverse) {
     DEMAND( in.size() > 0 );
 
     size_t dim = in.size();
@@ -185,7 +185,7 @@ qvector getDisceteFourierTransform(qvector in, vector<int> targs) {
     qindex len = getPow2(targs.size());
     qreal pi = 3.14159265358979323846;
     qreal a = 1 / std::sqrt(len);
-    qreal b = 2 * pi / len;
+    qreal b = (inverse ? -1 : 1) * 2 * pi / len;
 
     for (size_t i=0; i<dim; i++) {
         size_t x = getBitsAt(i, targs);

tests/utils/linalg.hpp

@@ -32,8 +32,8 @@ qindex getPow2(int);
 qreal getSum(vector<qreal> vec);
 qcomp getSum(qvector);
 qvector getNormalised(qvector);
-qvector getDisceteFourierTransform(qvector);
-qvector getDisceteFourierTransform(qvector in, vector<int> targs);
+qvector getDiscreteFourierTransform(qvector in, bool inverse);
+qvector getDiscreteFourierTransform(qvector in, vector<int> targs, bool inverse);
 
 qcomp   getInnerProduct(qvector bra, qvector ket);
 qmatrix getOuterProduct(qvector ket, qvector bra);