Added ApplyQPE and ApplyOperationPowerCA to Canon namespace (#2473)

- Added ApplyQPE and ApplyOperationPowerCA to Canon namespace.
- Added doc comments
- Added tests

ApplyQPE along with the previously implemented ApplyQFT are useful
primitives for quantum algorithms. To keep these implementations simple
and avoid performance penalty no swaps are performed on arguments or
results. In most cases such results can be used without problems. Note,
that ApplyQFT and ApplyQPE work with little-endian registers.

---------

Co-authored-by: Dmitry Vasilevsky <dmitryv@microsoft.com>

Author: DmitryVasilevsky

compiler/qsc_partial_eval/src/tests/qubits.rs

@@ -377,6 +377,6 @@ fn qubit_relabel_in_dynamic_block_triggers_capability_error() {
 
     assert_error(
         &error,
-        &expect!["CapabilityError(UseOfDynamicQubit(Span { lo: 60173, hi: 60186 }))"],
+        &expect!["CapabilityError(UseOfDynamicQubit(Span { lo: 62782, hi: 62795 }))"],
     );
 }

library/src/tests/canon.rs

@@ -504,6 +504,44 @@ fn check_qft_le_sample_4() {
     );
 }
 
+const QPE_TEST_LIB: &str = include_str!("resources/src/qpe.qs");
+
+#[test]
+fn check_qpe_z() {
+    test_expression_with_lib(
+        "Test.TestQPE_Z()",
+        QPE_TEST_LIB,
+        &Value::Tuple(vec![].into()),
+    );
+}
+
+#[test]
+fn check_qpe_s() {
+    test_expression_with_lib(
+        "Test.TestQPE_S()",
+        QPE_TEST_LIB,
+        &Value::Tuple(vec![].into()),
+    );
+}
+
+#[test]
+fn check_qpe_t() {
+    test_expression_with_lib(
+        "Test.TestQPE_T()",
+        QPE_TEST_LIB,
+        &Value::Tuple(vec![].into()),
+    );
+}
+
+#[test]
+fn check_qpe_p() {
+    test_expression_with_lib(
+        "Test.TestQPE_P()",
+        QPE_TEST_LIB,
+        &Value::Tuple(vec![].into()),
+    );
+}
+
 #[test]
 fn check_swap_reverse_register() {
     test_expression(
@@ -571,3 +609,18 @@ fn check_apply_operation_power_a() {
         &Value::RESULT_ZERO,
     );
 }
+
+#[test]
+fn check_apply_operation_power_ca() {
+    test_expression(
+        {
+            "{
+            use q = Qubit();
+            ApplyOperationPowerCA(12, Rx(Std.Math.PI()/16.0, _), q);
+            ApplyOperationPowerCA(-3, Rx(Std.Math.PI()/4.0, _), q);
+            M(q)
+        }"
+        },
+        &Value::RESULT_ZERO,
+    );
+}

library/src/tests/resources/src/qpe.qs

@@ -0,0 +1,50 @@
+namespace Test {
+    import Std.Diagnostics.*;
+
+    operation TestQPE_Phase(phaseGate : Qubit => Unit is Ctl + Adj, power : Int) : Unit {
+        use state = Qubit();
+        use phase = Qubit[4];
+
+        // Estimate eigenvalue of a phase gate on eigenvector of |0⟩
+        ApplyQPE(ApplyOperationPowerCA(_, qs => phaseGate(qs[0]), _), [state], phase);
+        // Eigenvalue should be 1 = exp(i * 2π * 0.0000), so the estimaped phase
+        // should be 0.0000 fraction of 2π.
+        Fact(CheckAllZero(phase), "Expected state |0000⟩");
+
+        // Estimate eigenvalue of a phase gate on eigenvector of |1⟩
+        X(state);
+        ApplyQPE(ApplyOperationPowerCA(_, qs => phaseGate(qs[0]), _), [state], phase);
+        // The eigenvalue of a phase gate on eigenvector |1⟩ is exp(i * 2π / 2^power),
+        // so the eigenvalue phase is the binary fraction 0.0…01 of 2π,
+        // where the 1 is at the position `power` after the point (counting from 1).
+        // So the a little-endian register `phase` should be
+        // phase[0] = 0, phase[1] = 0, …, phase[N-power] = 1, phase[N-power+1] = 0, …
+        X(phase[Length(phase)-power]);
+        Fact(CheckAllZero(phase), $"Incorrect phase for exp(i * 2π / 2^{power}).");
+
+        Reset(state);
+    }
+
+    operation TestQPE_Z() : Unit {
+        TestQPE_Phase(Z, 1); // eigenvalue = exp(i * 2π / 2^1)
+    }
+
+    operation TestQPE_S() : Unit {
+        TestQPE_Phase(S, 2); // eigenvalue = exp(i * 2π / 2^2)
+
+    }
+
+    operation TestQPE_T() : Unit {
+        TestQPE_Phase(T, 3); // eigenvalue = exp(i * 2π / 2^3)
+    }
+
+    // Phase gate is a rotation around the Z axis and an ajustment for the global phase.
+    operation P(phase : Double, q : Qubit) : Unit is Ctl + Adj {
+        Rz(phase, q);
+        Exp([], phase / 2.0, []);
+    }
+
+    operation TestQPE_P() : Unit {
+        TestQPE_Phase(P(Std.Math.PI() / 8.0, _), 4); // eigenvalue = exp(i * 2π / 2^4)
+    }
+}

library/std/src/Std/Canon.qs

@@ -599,6 +599,72 @@ operation ApplyOperationPowerA<'T>(
     }
 }
 
+/// # Summary
+/// Applies operation `op` to the `target` `power` times.
+/// If `power` is negative, the adjoint of `op` is used.
+/// If `power` is 0, the operation `op` is not applied.
+operation ApplyOperationPowerCA<'T>(
+    power : Int,
+    op : 'T => Unit is Ctl + Adj,
+    target : 'T
+) : Unit is Ctl + Adj {
+    let u = if power >= 0 { op } else { Adjoint op };
+    for _ in 1..AbsI(power) {
+        u(target);
+    }
+}
+
+/// # Summary
+/// Performs the quantum phase estimation algorithm for a unitary U represented
+/// by `applyPowerOfU`, and a `targetState`. Returns the phase estimation
+/// in the range of [0, 2π) as a fraction of 2π in the little-endian register `phase`.
+///
+/// # Input
+/// ## applyPowerOfU
+/// An oracle implementing Uᵐ for a unitary U and a given integer power m.
+/// `ApplyOperationPowerCA(_, U, _)` can be used to implement this oracle,
+/// if no better performing implementation is available.
+/// ## targetState
+/// A quantum register acted on by U. When this `targetState` is the eigenstate
+/// of the operator U, the corresponding eigenvalue is estimated.
+/// ## phase
+/// A little-endian quantum register containing the phase estimation result. The phase is
+/// a fraction of 2π represented in binary with q[0] containing 2π/2=π bit, q[1] containing 2π/4=π/2 bit,
+/// q[2] containing 2π/8=π/4 bit, and so on. The length of the register indicates the desired precision.
+/// The phase register is assumed to be in the state |0...0⟩ when the operation is invoked.
+///
+/// # Remarks
+/// All eigenvalues of a unitary operator are of magnitude 1 and therefore can be represented as
+/// $e^{i\phi}$ for some $\phi \in [0, 2\pi)$. Finding the phase of an eigenvalue is therefore
+/// equivalent to finding the eigenvalue itself. Passing the eigenvector as `targetState`
+/// to the phase estimation operation allows finding the eigenvalue to the desired precision
+/// determined by the length of the `phase` register.
+///
+/// # Reference
+/// - [Quantum phase estimation algorithm](https://wikipedia.org/wiki/Quantum_phase_estimation_algorithm)
+///
+/// # See Also
+/// - [Std.Canon.ApplyOperationPowerCA](xref:Qdk.Std.Canon.ApplyOperationPowerCA)
+operation ApplyQPE(
+    applyPowerOfU : (Int, Qubit[]) => Unit is Adj + Ctl,
+    targetState : Qubit[],
+    phase : Qubit[]
+) : Unit is Adj + Ctl {
+
+    let nQubits = Length(phase);
+    ApplyToEachCA(H, phase);
+
+    for i in 0..nQubits - 1 {
+        let power = 2^((nQubits - i) - 1);
+        Controlled applyPowerOfU(
+            phase[i..i],
+            (power, targetState)
+        );
+    }
+
+    Adjoint ApplyQFT(phase);
+}
+
 /// # Summary
 /// Relabels the qubits in the `current` array with the qubits in the `updated` array. The `updated` array
 /// must be a valid permutation of the `current` array.
@@ -660,8 +726,10 @@ export
     ApplyControlledOnInt,
     ApplyControlledOnBitString,
     ApplyQFT,
+    ApplyQPE,
     SwapReverseRegister,
     ApplyXorInPlace,
     ApplyXorInPlaceL,
     ApplyOperationPowerA,
+    ApplyOperationPowerCA,
     Relabel;