Quantum Phase Estimation in Q#

Quantum Phase Estimation (QPE) is one of the most important subroutines in quantum computing. It appears as a core component in Shor’s factoring algorithm, the HHL algorithm for linear systems, and quantum chemistry simulations. Given a unitary operator U and one of its eigenstates |psi>, QPE estimates the phase phi in the eigenvalue equation U|psi> = e^(2pii*phi)|psi>.

This tutorial walks through the mathematical foundations, builds the algorithm from scratch in Q#, analyzes its error behavior, and connects QPE to its major applications in factoring and chemistry.

Mathematical Foundations

Eigenvalues of Unitary Operators

Every unitary operator U satisfies U*U^dagger = I, which constrains its eigenvalues. If U|psi> = lambda|psi>, then unitarity requires |lambda|^2 = 1, so every eigenvalue has magnitude 1 and can be written as:

lambda = e^(2pii*phi)

where phi is a real number in [0, 1). We call phi the phase of the eigenvalue. QPE takes as input a unitary U and a state close to one of its eigenstates, and it outputs an n-bit approximation of phi.

The Phase Kickback Mechanism

The core idea behind QPE is phase kickback. When a controlled-U gate acts on a control qubit in state |1> and a target register in eigenstate |psi>, the eigenvalue appears as a phase on the control qubit:

CU(|1>|psi>) = |1>(U|psi>) = |1>(e^(2piiphi)|psi>) = e^(2piiphi)|1>|psi>

The target register is unchanged, but the control qubit picks up the phase factor. This is the mechanism that transfers information about phi from the target register into the counting register.

Derivation of the QPE State

Consider n counting qubits and a target register prepared in eigenstate |psi>. After applying Hadamard gates to all counting qubits, the state is:

|state_1> = (1/sqrt(2^n)) * sum_{k=0}^{2^n - 1} |k> tensor |psi>

Each counting qubit j (for j = 0, 1, …, n-1) controls the application of U^(2^j). The controlled-U^(2^j) gate on qubit j acts as:

|0>|psi> -> |0>|psi> |1>|psi> -> e^(2piiphi2^j)|1>|psi>

After all controlled applications, the state of the counting register becomes:

|state_2> = (1/sqrt(2^n)) * sum_{k=0}^{2^n - 1} e^(2piiphik) |k> tensor |psi>

To see why, note that the integer k has binary representation k = k_{n-1} * 2^{n-1} + … + k_1 * 2 + k_0. Each bit k_j contributes a phase factor of e^(2piiphi2^j) when k_j = 1, and these multiply together to give e^(2piiphik).

The Inverse QFT Extracts the Phase

The Quantum Fourier Transform maps computational basis states to frequency states:

QFT|m> = (1/sqrt(2^n)) * sum_{k=0}^{2^n - 1} e^(2piimk/2^n) |k>

Comparing this with our state after the controlled unitaries:

|state_2> = (1/sqrt(2^n)) * sum_{k=0}^{2^n - 1} e^(2piiphik) |k>

If phi = m/2^n for some integer m, then |state_2> is exactly QFT|m>. Applying the inverse QFT recovers |m>, and measuring gives m. The estimated phase is then phi = m/2^n.

Precision and the Non-Exact Case

The counting register has n qubits, so it can represent 2^n distinct values. The phase grid consists of the points {0, 1/2^n, 2/2^n, …, (2^n - 1)/2^n}, giving a precision of 1/2^n.

When the true phase phi does not equal any grid point exactly, QPE returns a probabilistic result. Let phi = (m + delta)/2^n where m is the nearest integer below and 0 < delta < 1. The probability of measuring outcome j is:

P(j) = (1/2^n) * |sin(pi*(phi2^n - j)) / sin(pi(phi*2^n - j)/2^n)|^2

This is a peaked distribution centered near the true value. The probability of obtaining one of the two nearest grid points (m or m+1) satisfies:

P(m) + P(m+1) >= 1 - 1/(2*(2^n - 1))

For n = 3 this gives a combined probability of at least approximately 0.96. The estimate is therefore accurate to within 1/2^n with high probability.

The Quantum Fourier Transform

Before implementing QPE, we need the QFT and its inverse. The QFT is a quantum analogue of the discrete Fourier transform, operating on quantum amplitudes.

QFT Circuit Structure

For n qubits, the QFT applies the following sequence of gates:

1. For qubit 0 (most significant): apply H, then controlled-R_2 (controlled by qubit 1), controlled-R_3 (controlled by qubit 2), and so on through controlled-R_n (controlled by qubit n-1).
2. For qubit 1: apply H, then controlled-R_2 through controlled-R_{n-1}.
3. Continue this pattern for all qubits.
4. Finally, swap the qubit ordering to match the standard convention.

Here R_k = diag(1, e^(2pii/2^k)) is a phase rotation gate. In Q#, R1(theta, qubit) applies diag(1, e^(itheta)), so R_k corresponds to R1(2pi/2^k, qubit).

QFT Matrix for n = 2

For two qubits, the QFT matrix is:

QFT_2 = (1/2) * [[1, 1, 1, 1], [1, i, -1, -i], [1, -1, 1, -1], [1, -i, -1, i]]

Each entry (j, k) equals (1/2) * e^(2piijk/4). The columns are orthonormal, confirming that the QFT is unitary.

Why Swaps Are Needed

The QFT circuit naturally produces output in bit-reversed order. The controlled rotation structure processes the most significant bit first but leaves it in the first qubit position, while the standard QFT definition places the most significant bit last. The swap operations at the end (or beginning, depending on convention) fix this ordering.

Why Adjoint QFT = Inverse QFT

Because the QFT is unitary, its inverse equals its adjoint (conjugate transpose). In Q#, any operation annotated with is Adj automatically supports the Adjoint functor, which reverses the gate sequence and conjugates all rotation angles. This means we can write a forward QFT with the Adj annotation and obtain the inverse QFT for free by calling Adjoint QFT(register).

QFT Implementation in Q#

QFT.qs

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Applies the Quantum Fourier Transform to a register.
    /// Annotated with Adj + Ctl so Adjoint QFT gives the inverse QFT.
    operation QFT(register : Qubit[]) : Unit is Adj + Ctl {
        let n = Length(register);

        for i in 0..n - 1 {
            // Apply Hadamard to qubit i
            H(register[i]);

            // Apply controlled R_k rotations
            for j in i + 1..n - 1 {
                let k = j - i + 1;
                let angle = 2.0 * PI() / IntAsDouble(1 <<< k);
                Controlled R1([register[j]], (angle, register[i]));
            }
        }

        // Swap qubits to match the standard QFT output ordering
        for i in 0..n / 2 - 1 {
            SWAP(register[i], register[n - 1 - i]);
        }
    }
}

With this definition, Adjoint QFT(register) applies the inverse QFT. The Q# compiler reverses the gate order, negates all rotation angles, and reverses the swap pattern automatically.

How the Algorithm Works

QPE uses two registers: a counting register of n qubits (which determines the precision of the estimate) and a target register holding the eigenstate |psi>. The algorithm proceeds in three stages:

1. Superposition: Apply Hadamard gates to every qubit in the counting register, creating an equal superposition of all 2^n basis states.
2. Controlled unitary applications: Apply controlled-U^(2^k) operations, where counting qubit k controls the application of U raised to the power 2^k on the target register. This encodes the phase into the relative phases of the counting register.
3. Inverse Quantum Fourier Transform: Apply the inverse QFT to the counting register. This converts the phase information from the frequency domain into a computational basis state that can be measured directly.

After measurement, the counting register contains an n-bit binary approximation of phi. With n counting qubits, the estimate has precision 1/2^n.

Worked Example: Phase of the T Gate

The T gate applies a phase of pi/4 to the |1> state:

T|1> = e^(ipi/4)|1> = e^(2pii(1/8))|1>

So the phase we expect QPE to find is phi = 1/8 = 0.125. With 3 counting qubits, we can represent this exactly since 1/8 = 1/2^3 = 0.001 in binary.

Worked Example: Phase of the S Gate

The S gate (also called the phase gate or P gate) applies a phase of pi/2 to the |1> state:

S|1> = e^(ipi/2)|1> = e^(2pii(1/4))|1>

The phase is phi = 1/4 = 0.25. In binary, 1/4 = 0.01, which requires exactly 2 bits.

S Gate with 2 Counting Qubits (Exact Representation)

With n = 2 counting qubits, the grid points are {0, 1/4, 2/4, 3/4}. The phase 1/4 falls exactly on the grid point m = 1. QPE returns m = 1 with certainty, and the estimated phase is 1/4 = 0.25.

S Gate with 1 Counting Qubit (Insufficient Precision)

With only n = 1 counting qubit, the grid points are {0, 1/2}. The phase 1/4 sits exactly halfway between grid points 0 and 1/2. In this case, the QPE circuit reduces to:

1. Apply H to the single counting qubit, producing (|0> + |1>)/sqrt(2).
2. Apply controlled-S to the target. The counting qubit picks up phase: (|0> + e^(2pii/4)|1>)/sqrt(2) = (|0> + i|1>)/sqrt(2).
3. Apply inverse QFT (which is just H for one qubit): H maps this to ((1+i)|0> + (1-i)|1>)/(2).

The probability of measuring 0 is |1+i|^2/4 = 2/4 = 1/2, and the probability of measuring 1 is |1-i|^2/4 = 2/4 = 1/2. QPE gives 0 or 1 with equal probability, corresponding to phase estimates of 0 or 0.5, neither of which is correct. This illustrates concretely why sufficient counting qubits are essential.

ControlledSPower Operation

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Applies controlled-S^(2^power) where S has phase pi/2.
    /// S^(2^k) applies a phase of e^(i * pi/2 * 2^k) to |1>.
    operation ControlledSPower(
        control : Qubit,
        target : Qubit,
        power : Int
    ) : Unit {
        // S gate has angle pi/2. S^(2^k) has angle pi/2 * 2^k.
        let angle = PI() / 2.0 * IntAsDouble(1 <<< power);
        Controlled R1([control], (angle, target));
    }

    /// QPE for the S gate with configurable counting qubits.
    operation RunQPE_SGate(nCounting : Int) : Double {
        use countingRegister = Qubit[nCounting];
        use target = Qubit();

        // Prepare eigenstate |1>
        X(target);

        // Hadamard on counting qubits
        for q in countingRegister {
            H(q);
        }

        // Controlled-S^(2^k) applications
        for k in 0..nCounting - 1 {
            ControlledSPower(countingRegister[k], target, k);
        }

        // Inverse QFT
        InverseQFT(countingRegister);

        // Measure counting register
        mutable result = 0;
        for k in 0..nCounting - 1 {
            let bit = M(countingRegister[k]);
            if bit == One {
                set result = result + (1 <<< k);
            }
        }

        ResetAll(countingRegister);
        Reset(target);

        return IntAsDouble(result) / IntAsDouble(1 <<< nCounting);
    }
}

Running RunQPE_SGate(2) returns 0.25 every time. Running RunQPE_SGate(1) returns either 0.0 or 0.5 with equal probability.

Using EstimateFrequencyA

Q# provides built-in operations for frequency estimation. The EstimateFrequencyA operation from the Microsoft.Quantum.Characterization namespace repeatedly prepares a state, applies an operation, and measures to estimate the probability of a given outcome. While not a full QPE implementation, it demonstrates the estimation concept.

EstimatePhaseSimple.qs

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Measurement;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Estimates the phase of the T gate by repeated sampling.
    /// Prepares |1>, applies T, then measures in a rotated basis.
    operation EstimateTPhaseByFrequency(nTrials : Int) : Double {
        mutable countOne = 0;

        for _ in 1..nTrials {
            use q = Qubit();
            X(q);               // Prepare |1> (eigenstate of T)
            T(q);               // Apply T: introduces phase e^(i*pi/4)
            H(q);               // Rotate to measure phase information
            let result = M(q);
            if result == One {
                set countOne += 1;
            }
            Reset(q);
        }

        // The probability of measuring One encodes the phase
        let probOne = IntAsDouble(countOne) / IntAsDouble(nTrials);
        return probOne;
    }
}

This sampling approach gives a statistical estimate. For an exact circuit-based approach, we need to build the full QPE circuit.

Building QPE from Scratch in Q#

The following implementation constructs the complete QPE circuit with controlled power-of-two applications of a unitary and an inverse QFT.

QPE.qs

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Measurement;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;
    open Microsoft.Quantum.Arrays;

    /// Applies the inverse Quantum Fourier Transform to a register.
    operation InverseQFT(register : Qubit[]) : Unit is Adj + Ctl {
        let n = Length(register);

        // Reverse qubit order for the QFT convention
        for i in 0..n/2 - 1 {
            SWAP(register[i], register[n - 1 - i]);
        }

        for i in 0..n - 1 {
            // Apply controlled rotations from previous qubits
            for j in 0..i - 1 {
                let angle = -PI() / IntAsDouble(1 <<< (i - j));
                Controlled R1([register[j]], (angle, register[i]));
            }
            H(register[i]);
        }
    }

    /// Applies controlled-U^(2^power) where U is the T gate.
    /// T^(2^k) applies a phase of e^(i * pi/4 * 2^k) to |1>.
    operation ControlledTPower(
        control : Qubit,
        target : Qubit,
        power : Int
    ) : Unit {
        // T gate has angle pi/4. T^(2^k) has angle pi/4 * 2^k.
        let angle = PI() / 4.0 * IntAsDouble(1 <<< power);
        Controlled R1([control], (angle, target));
    }

    /// Full Quantum Phase Estimation for the T gate.
    /// Uses nCounting qubits for precision.
    /// Returns the estimated phase as a fraction of 2*pi.
    operation RunQPE(nCounting : Int) : Double {
        use countingRegister = Qubit[nCounting];
        use target = Qubit();

        // Step 1: Prepare eigenstate |1> in target register
        X(target);

        // Step 2: Apply Hadamard to all counting qubits
        for q in countingRegister {
            H(q);
        }

        // Step 3: Apply controlled-U^(2^k) operations
        // Counting qubit k controls U^(2^k)
        for k in 0..nCounting - 1 {
            ControlledTPower(countingRegister[k], target, k);
        }

        // Step 4: Apply inverse QFT to counting register
        InverseQFT(countingRegister);

        // Step 5: Measure the counting register
        mutable result = 0;
        for k in 0..nCounting - 1 {
            let bit = M(countingRegister[k]);
            if bit == One {
                set result = result + (1 <<< k);
            }
        }

        // Clean up
        ResetAll(countingRegister);
        Reset(target);

        // Convert measured integer to phase estimate
        // phase = measured_value / 2^n
        return IntAsDouble(result) / IntAsDouble(1 <<< nCounting);
    }
}

Python host (run_qpe.py)

import qsharp
from PhaseEstimation import RunQPE

# Run QPE with 3 counting qubits (precision = 1/8)
phase = RunQPE.simulate(nCounting=3)
print(f"Estimated phase: {phase}")

# Run with higher precision
phase_5 = RunQPE.simulate(nCounting=5)
print(f"Estimated phase (5 qubits): {phase_5}")

Expected Output

Estimated phase: 0.125
Estimated phase (5 qubits): 0.125

With 3 counting qubits, the T gate phase of 1/8 is represented exactly in binary (0.001), so QPE returns the exact value 0.125 every time. This is not a coincidence; when the true phase aligns with one of the 2^n grid points, QPE succeeds with certainty.

Interpreting Results

The measured integer from the counting register encodes the phase as a binary fraction. If you measure the integer m using n counting qubits, the estimated phase is:

phi = m / 2^n

For the T gate example with n = 3, measuring m = 1 gives phi = 1/8 = 0.125, which corresponds to the eigenvalue e^(2pii/8) = e^(i*pi/4).

When the true phase does not fall exactly on the grid, the measurement outcome is probabilistic. The probability concentrates on the two nearest grid points, and the estimate is accurate to within 1/2^n. Increasing the number of counting qubits improves precision exponentially.

Phase Error Analysis

When the true phase does not land on a grid point, QPE returns a probability distribution over all possible outcomes. Understanding this distribution is essential for reasoning about QPE accuracy in practice.

Analytical Example: phi = 0.1 with n = 3

With n = 3 counting qubits, the grid has 8 points: {0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8}. Suppose the true phase is phi = 0.1, which falls between grid points 0 (= 0.0) and 1/8 (= 0.125).

The quantity phi * 2^n = 0.1 * 8 = 0.8. The probability of measuring outcome j is:

P(j) = (1/(2^n)^2) * |sin(pi * (0.8 - j)) / sin(pi * (0.8 - j) / 2^n)|^2

Computing each probability:

Outcome j	Phase estimate j/8	P(j)
0	0.000	0.057
1	0.125	0.877
2	0.250	0.026
3	0.375	0.0048
4	0.500	0.0025
5	0.625	0.0025
6	0.750	0.0048
7	0.875	0.0048

The two nearest grid points (j = 0 and j = 1) together capture about 93.4% of the probability. The distribution is strongly peaked, confirming that QPE gives a useful estimate even when the phase is not exactly representable.

Q# Code: Histogramming QPE Results

The following code runs QPE for a custom phase angle and collects statistics over many shots.

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Applies controlled-Rz^(2^power) where Rz has a specified angle.
    /// This simulates QPE on a unitary with an arbitrary phase.
    operation ControlledPhasePower(
        control : Qubit,
        target : Qubit,
        phaseAngle : Double,
        power : Int
    ) : Unit {
        let angle = phaseAngle * IntAsDouble(1 <<< power);
        Controlled R1([control], (angle, target));
    }

    /// Runs QPE once for an arbitrary phase and returns the raw integer.
    operation RunQPESingleShot(
        nCounting : Int,
        phaseAngle : Double
    ) : Int {
        use countingRegister = Qubit[nCounting];
        use target = Qubit();

        // Prepare eigenstate |1>
        X(target);

        // Hadamard on all counting qubits
        for q in countingRegister {
            H(q);
        }

        // Controlled phase applications
        for k in 0..nCounting - 1 {
            ControlledPhasePower(
                countingRegister[k], target, phaseAngle, k
            );
        }

        // Inverse QFT
        Adjoint QFT(countingRegister);

        // Measure
        mutable result = 0;
        for k in 0..nCounting - 1 {
            let bit = M(countingRegister[k]);
            if bit == One {
                set result = result + (1 <<< k);
            }
        }

        ResetAll(countingRegister);
        Reset(target);

        return result;
    }

    /// Runs QPE many times and returns an array of measured integers.
    operation RunQPEHistogram(
        nCounting : Int,
        phaseAngle : Double,
        nShots : Int
    ) : Int[] {
        mutable results = [0, size = nShots];
        for i in 0..nShots - 1 {
            set results w/= i <- RunQPESingleShot(nCounting, phaseAngle);
        }
        return results;
    }
}

Python host for histogram analysis:

import qsharp
from PhaseEstimation import RunQPEHistogram
from collections import Counter
import math

# Phase of 0.1 * 2*pi
phi = 0.1
phase_angle = 2.0 * math.pi * phi
n_counting = 3
n_shots = 1000

results = RunQPEHistogram.simulate(
    nCounting=n_counting,
    phaseAngle=phase_angle,
    nShots=n_shots
)

counts = Counter(results)
print(f"QPE results for phi = {phi}, n = {n_counting}:")
for outcome in sorted(counts.keys()):
    phase_est = outcome / (2**n_counting)
    freq = counts[outcome] / n_shots
    print(f"  outcome {outcome} (phase {phase_est:.4f}): {freq:.3f}")

A typical run produces output showing that outcomes 0 and 1 dominate, consistent with the analytical distribution above.

Generalized QPE for Arbitrary Unitaries

The examples so far hardcode a specific gate (T or S). A more flexible implementation accepts any single-qubit unitary as a parameter.

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Generalized QPE that accepts any single-qubit unitary.
    /// The unitary is passed as an operation parameter.
    /// The caller must also provide a way to apply U^(2^k)
    /// via the unitaryPower parameter.
    operation GeneralizedQPE(
        nCounting : Int,
        prepareEigenstate : (Qubit => Unit),
        unitaryPower : ((Qubit, Qubit, Int) => Unit)
    ) : Double {
        use countingRegister = Qubit[nCounting];
        use target = Qubit();

        // Prepare the eigenstate in the target register
        prepareEigenstate(target);

        // Hadamard on all counting qubits
        for q in countingRegister {
            H(q);
        }

        // Controlled-U^(2^k) for each counting qubit
        for k in 0..nCounting - 1 {
            unitaryPower(countingRegister[k], target, k);
        }

        // Inverse QFT
        InverseQFT(countingRegister);

        // Measure
        mutable result = 0;
        for k in 0..nCounting - 1 {
            let bit = M(countingRegister[k]);
            if bit == One {
                set result = result + (1 <<< k);
            }
        }

        ResetAll(countingRegister);
        Reset(target);

        return IntAsDouble(result) / IntAsDouble(1 <<< nCounting);
    }

    // --- Example callers ---

    /// QPE for the T gate using the generalized implementation.
    operation RunGeneralizedQPE_T() : Double {
        return GeneralizedQPE(
            3,
            X,   // Prepare |1>
            ControlledTPower
        );
    }

    /// QPE for the S gate using the generalized implementation.
    operation RunGeneralizedQPE_S() : Double {
        return GeneralizedQPE(
            3,
            X,
            ControlledSPower
        );
    }

    /// Controlled-Rz(pi/3)^(2^power).
    /// Rz(theta)|1> = e^(-i*theta/2)|1>, so the phase is -theta/(4*pi)
    /// mod 1. For theta = pi/3, phi = -1/12 mod 1 = 11/12.
    operation ControlledRzPiOver3Power(
        control : Qubit,
        target : Qubit,
        power : Int
    ) : Unit {
        let angle = -PI() / 3.0 * IntAsDouble(1 <<< power);
        Controlled R1([control], (angle, target));
    }

    /// QPE for Rz(pi/3).
    /// The expected phase is 11/12 = 0.91666...
    operation RunGeneralizedQPE_Rz() : Double {
        return GeneralizedQPE(
            5,
            X,
            ControlledRzPiOver3Power
        );
    }
}

The key design challenge for generalized QPE is implementing U^(2^k) efficiently. For phase gates, raising to a power simply multiplies the angle. For an arbitrary unitary, computing U^(2^k) might require applying U sequentially 2^k times, which is exponentially expensive. Efficient implementations exploit the structure of U to avoid this blowup. For example, modular exponentiation in Shor’s algorithm uses repeated squaring at the circuit level.

Iterative Phase Estimation

Standard QPE uses n counting qubits simultaneously. Iterative Phase Estimation (IPE) replaces these with a single counting qubit that is reused across n rounds, extracting one bit of the phase per round. IPE requires mid-circuit measurement and classical feedback, but it uses far fewer qubits.

How IPE Works

IPE extracts the phase bits from least significant to most significant. In round j (for j = n-1, n-2, …, 0):

1. Prepare the counting qubit in |0> and apply H.
2. Apply controlled-U^(2^j) between the counting qubit and the target.
3. Apply a phase correction based on the bits already measured in previous rounds. The correction angle is theta_j = -pi * sum_{l=1}^{n-1-j} (bit_{j+l} / 2^l), where bit_{j+l} is the phase bit extracted in the round that measured bit position j+l.
4. Apply R1(theta_j) to the counting qubit, then apply H and measure.
5. The measurement result is the j-th bit of the phase.

This classical feedback loop is what distinguishes IPE from standard QPE. Each measurement collapses the counting qubit, and the correction angle steers subsequent rounds to extract the next bit.

Q# Implementation of IPE

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Iterative Phase Estimation using a single counting qubit.
    /// Extracts n bits of the phase of the T gate.
    operation IterativePhaseEstimation(n : Int) : Double {
        use counting = Qubit();
        use target = Qubit();

        // Prepare eigenstate |1>
        X(target);

        // Array to store measured bits (index 0 = least significant)
        mutable phaseBits = [false, size = n];

        // Iterate from least significant bit to most significant
        // Round j extracts bit position (n - 1 - j) when counting
        // from the MSB, but we iterate to fill phaseBits from MSB down.
        for round in 0..n - 1 {
            let bitIndex = n - 1 - round;

            // Prepare counting qubit in |+> state
            H(counting);

            // Apply controlled-U^(2^bitIndex)
            ControlledTPower(counting, target, bitIndex);

            // Apply phase correction from previously measured bits
            mutable correctionAngle = 0.0;
            for l in 1..round {
                if phaseBits[n - 1 - (round - l)] {
                    set correctionAngle =
                        correctionAngle - PI() / IntAsDouble(1 <<< l);
                }
            }
            R1(correctionAngle, counting);

            // Measure in X basis (H then measure in Z)
            H(counting);
            let bit = M(counting);
            set phaseBits w/= bitIndex <- (bit == One);

            // Reset for next round
            Reset(counting);
        }

        Reset(target);

        // Convert bit array to phase value
        mutable phase = 0.0;
        for k in 0..n - 1 {
            if phaseBits[k] {
                set phase = phase +
                    1.0 / IntAsDouble(1 <<< (k + 1));
            }
        }

        return phase;
    }
}

When to Use IPE Over Standard QPE

IPE is preferred in several scenarios:

• Qubit-limited hardware: When the device has very few qubits, using a single counting qubit instead of n is a significant advantage.
• Mid-circuit measurement support: IPE requires measuring and resetting qubits during the circuit. Hardware that supports this (such as trapped-ion systems and some superconducting processors) can run IPE natively.
• Shallow circuits: Each round of IPE has depth proportional to 2^j for the controlled-U application, but the rounds run sequentially rather than simultaneously. The maximum circuit depth in any single round is 2^(n-1), the same as in standard QPE.

The tradeoff is that IPE requires n sequential rounds of quantum operations with classical processing in between, making it slower in wall-clock time than standard QPE, which runs all counting qubits in parallel.

Connection to Shor’s Algorithm

QPE is the core subroutine inside Shor’s algorithm for integer factoring. Understanding this connection reveals why QPE has such profound implications for cryptography.

Order Finding via QPE

To factor an integer N, Shor’s algorithm reduces the problem to order finding: given a random integer a coprime to N, find the smallest positive integer r such that a^r = 1 (mod N). This r is called the order of a modulo N.

Shor’s algorithm constructs a unitary operator U_a that acts on a register of ceil(log2(N)) qubits:

U_a|x> = |a*x mod N> for x in {0, 1, …, N-1}

Eigenstates and Eigenvalues of U_a

The eigenstates of U_a are the “period states”:

|phi_s> = (1/sqrt(r)) * sum_{k=0}^{r-1} e^(-2piisk/r) |a^k mod N>

for s = 0, 1, …, r-1. The corresponding eigenvalues are:

U_a|phi_s> = e^(2pii*s/r) |phi_s>

QPE on U_a with eigenstate |phi_s> estimates the phase s/r.

From Phase Estimate to Period

After QPE returns an approximation of s/r, the continued fraction algorithm recovers r. The continued fraction expansion of a rational number p/q (where p/q approximates s/r) converges to s/r when the QPE estimate has sufficient precision. Specifically, using n = 2*ceil(log2(N)) counting qubits provides enough precision to distinguish different fractions with denominator at most N.

The Hard Part: Modular Exponentiation

The controlled-U_a^(2^k) gates require computing a^(2^k) * x mod N in superposition. This is modular exponentiation, and it is the most expensive part of Shor’s algorithm. Each modular multiplication requires O(n^2) elementary gates (using schoolbook multiplication) or O(n * log(n) * log(log(n))) gates (using more advanced techniques), where n = ceil(log2(N)).

For a 2048-bit RSA modulus, the modular exponentiation circuit requires billions of T-gates, making this far beyond the reach of current hardware. Nonetheless, the polynomial scaling (in the number of bits of N) represents an exponential speedup over the best known classical factoring algorithms.

A Subtle Point About Eigenstates

In practice, we do not know the eigenstates |phi_s> because computing them requires knowing r (which is what we are trying to find). Fortunately, the uniform superposition state |1> can be decomposed as:

|1> = (1/sqrt(r)) * sum_{s=0}^{r-1} |phi_s>

Running QPE with |1> as the target produces a superposition of estimates s/r for random s. Each measurement gives a different s, but the continued fraction algorithm extracts r from any of them (with high probability).

Connection to Quantum Chemistry

QPE is the foundation of quantum algorithms for computing molecular energies, making it central to the promise of quantum advantage in chemistry and materials science.

Energy Estimation via Phase Estimation

A molecular Hamiltonian H has eigenstates |E_k> with energies E_k:

H|E_k> = E_k|E_k>

The time evolution operator U(t) = e^(-iHt) is a unitary whose eigenstates are the same |E_k>, with eigenvalues:

U(t)|E_k> = e^(-iE_kt)|E_k>

Comparing with the standard QPE eigenvalue form e^(2piiphi), we identify phi = -E_kt/(2pi). QPE on U(t) estimates phi, from which we recover E_k = -2pi*phi/t.

The QPE Workflow for Chemistry

The chemistry QPE workflow proceeds as follows:

1. Prepare a trial state: Start with a classically computable approximation to the ground state, such as a Hartree-Fock state or a state prepared via a Unitary Coupled Cluster (UCC) ansatz. The trial state should have significant overlap with the true ground state.

2. Implement U(t) via Trotterization: The Hamiltonian H is a sum of terms H = H_1 + H_2 + … + H_L. The Trotter-Suzuki decomposition approximates e^(-iHt) as a product of simpler exponentials:

   e^(-iHt) approx (e^(-iH_1dt) * e^(-iH_2dt) * … * e^(-iH_Ldt))^(t/dt)

   Each e^(-iH_jdt) can be implemented with elementary quantum gates.

3. Run QPE: Use the trial state as the target register and the Trotterized U(t) as the unitary in QPE.

4. Read out the energy: Convert the measured phase to an energy using E = -2piphi/t.

The Precision-Depth Tradeoff

Larger values of t give better energy resolution (since the energy maps to a larger phase), but they require deeper Trotter circuits that accumulate more gate errors. On near-term hardware, this tradeoff severely limits the achievable precision. For a target energy accuracy of epsilon, the required evolution time scales as t ~ 1/epsilon, and the Trotter circuit depth scales as t^2/epsilon_trotter, where epsilon_trotter is the Trotter error tolerance.

Q# Structure for Chemistry QPE

While a full chemistry QPE implementation requires a Hamiltonian simulation library, the structural outline in Q# looks like this:

namespace ChemistryQPE {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Placeholder: prepare a Hartree-Fock initial state.
    operation PrepareTrialState(register : Qubit[]) : Unit is Adj + Ctl {
        // In practice, this encodes the Hartree-Fock determinant
        // as a computational basis state. For H2 in minimal basis,
        // the Hartree-Fock state is |0011>.
        X(register[0]);
        X(register[1]);
    }

    /// Placeholder: apply one Trotter step of e^(-i*H*dt).
    /// A real implementation decomposes H into Pauli terms
    /// and applies Rxx, Ryy, Rzz rotations.
    operation TrotterStep(
        register : Qubit[],
        dt : Double
    ) : Unit is Adj + Ctl {
        // Placeholder: apply rotations corresponding to
        // each term in the Hamiltonian
        let nQubits = Length(register);
        for i in 0..nQubits - 2 {
            CNOT(register[i], register[i + 1]);
            Rz(dt, register[i + 1]);
            CNOT(register[i], register[i + 1]);
        }
    }

    /// Apply U(t) = e^(-i*H*t) via nSteps Trotter steps.
    operation TimeEvolution(
        register : Qubit[],
        t : Double,
        nSteps : Int
    ) : Unit is Adj + Ctl {
        let dt = t / IntAsDouble(nSteps);
        for _ in 1..nSteps {
            TrotterStep(register, dt);
        }
    }
}

The placeholder TrotterStep would be replaced by a real Hamiltonian simulation for a specific molecule. Libraries like Microsoft’s quantum chemistry library provide these decompositions for molecular Hamiltonians encoded in second quantization.

Resource Analysis

Understanding the gate cost of QPE is critical for determining when it becomes feasible on real hardware.

Gate Count Breakdown

For QPE on a single-qubit unitary with n counting qubits, the circuit contains:

1. Hadamard gates: n (one per counting qubit for the initial superposition).
2. Controlled unitary applications: Counting qubit k requires U^(2^k), which may need 2^k sequential applications of controlled-U. The total number of controlled-U applications is sum_{k=0}^{n-1} 2^k = 2^n - 1.
3. Inverse QFT gates: The QFT on n qubits uses n Hadamard gates and n*(n-1)/2 controlled rotation gates, totaling n + n*(n-1)/2 = n*(n+1)/2 gates.

For the T gate specifically, each controlled-T^(2^k) reduces to a single controlled-R1 rotation, so the “2^k applications” collapse to a single gate. The unitary application count is then just n (one per counting qubit).

Cost Table for T Gate QPE

n (counting qubits)	Controlled-U gates	QFT gates	Total gates	Phase precision
3	3	6	12	1/8 = 0.125
5	5	15	25	1/32 = 0.03125
10	10	55	75	1/1024
20	20	210	250	~10^(-6)

Cost Table for General Unitary QPE

When U^(2^k) cannot be simplified and requires 2^k applications of U:

n (counting qubits)	Controlled-U gates	QFT gates	Total gates	Phase precision
3	7	6	16	1/8
5	31	15	51	1/32
10	1023	55	1088	1/1024
20	~10^6	210	~10^6	~10^(-6)

The exponential growth in unitary applications is the primary cost driver. This is why efficient implementations of U^(2^k) (using repeated squaring or algebraic structure) are critical for practical QPE.

Phase Estimation on Real Hardware

On a noiseless simulator, QPE for the T gate with 3 counting qubits returns phi = 0.125 with certainty. Real hardware introduces gate errors, readout errors, and decoherence, which cause incorrect outcomes.

Majority Voting for Noise Robustness

A straightforward strategy for handling noise is to run QPE multiple times and take the most frequent result (the mode) as the estimate. This majority voting approach works well when the per-shot error probability is below 50%.

namespace PhaseEstimation {
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;
    open Microsoft.Quantum.Convert;

    /// Runs QPE nShots times and returns the mode (most common result).
    operation RunQPEWithMajorityVote(
        nCounting : Int,
        nShots : Int
    ) : Double {
        // Count occurrences of each outcome
        mutable counts = [0, size = 1 <<< nCounting];

        for _ in 1..nShots {
            use countingRegister = Qubit[nCounting];
            use target = Qubit();

            X(target);
            for q in countingRegister {
                H(q);
            }
            for k in 0..nCounting - 1 {
                ControlledTPower(countingRegister[k], target, k);
            }
            InverseQFT(countingRegister);

            mutable result = 0;
            for k in 0..nCounting - 1 {
                let bit = M(countingRegister[k]);
                if bit == One {
                    set result = result + (1 <<< k);
                }
            }

            set counts w/= result <- counts[result] + 1;

            ResetAll(countingRegister);
            Reset(target);
        }

        // Find the mode
        mutable maxCount = 0;
        mutable mode = 0;
        let nOutcomes = 1 <<< nCounting;
        for i in 0..nOutcomes - 1 {
            if counts[i] > maxCount {
                set maxCount = counts[i];
                set mode = i;
            }
        }

        return IntAsDouble(mode) / IntAsDouble(1 <<< nCounting);
    }
}

Hardware Limitations

For current quantum hardware, the practical limits of QPE are:

• Gate error rates: Superconducting devices have two-qubit gate error rates around 0.5-1%. A QPE circuit with n = 3 has roughly 12 gates, giving a total error probability of about 6-12%. Majority voting over 100 shots handles this reliably.
• Circuit depth: For n = 10, the circuit depth grows substantially. With QFT contributing O(n^2) = 100 gates and the unitary applications adding more, the total depth can exceed the coherence time of the qubits.
• Beyond n = 10: Circuits with more than roughly 10 counting qubits generally require error correction to produce reliable results on today’s hardware. This is the domain of fault-tolerant quantum computing.

Connecting to Azure Quantum Resource Estimator

The Azure Quantum Resource Estimator provides detailed physical resource estimates for quantum algorithms running on fault-tolerant hardware with surface code error correction. This helps answer the question: “When will QPE become practical for real problems?”

Estimating Resources for QPE

To use the resource estimator with a Q# program, you configure the target hardware model and error budget, then run the estimator against your Q# operation:

import qsharp
from qsharp.estimator import EstimatorParams, QubitParams

# Configure the resource estimator
params = EstimatorParams()
params.error_budget = 0.01  # 1% total error budget
params.qubit_params.name = QubitParams.GATE_NS_E3  # Superconducting, 10^-3 error rate

# Run the estimator on the QPE operation
result = qsharp.estimate(
    "PhaseEstimation.RunQPE(nCounting=10)",
    params=params
)

print(f"Logical qubits: {result['logicalCounts']['numQubits']}")
print(f"T-gate count: {result['logicalCounts']['tCount']}")
print(f"Physical qubits: {result['physicalCounts']['physicalQubits']}")
print(f"Runtime: {result['physicalCounts']['runtime']}")

Typical Resource Estimates

For different QPE configurations with surface code error correction:

Application	Logical qubits	T-gates	Physical qubits (est.)	Runtime (est.)
T gate QPE, n=10	11	~100	~5,000	milliseconds
Simple chemistry, n=20	~50	~10^6	~100,000	seconds
Shor’s (2048-bit RSA)	~4,000	~10^10	~20 million	hours

These estimates depend heavily on the assumed physical error rate, the error correction code distance, and the magic state distillation factory layout. The key takeaway is that while QPE for simple cases is feasible in the near-term fault-tolerant era, cryptographically relevant applications require hardware that is still years away.

Iterative vs. Standard QPE: A Comparison

Feature	Standard QPE	Iterative QPE
Counting qubits	n	1
Circuit depth (max)	O(2^n + n^2)	O(2^(n-1))
Measurement	Once at the end	n mid-circuit measures
Classical feedback	None	Required
Total shots	1 per estimate	1 per estimate
Hardware requirement	Many qubits	Mid-circuit measurement

Both variants achieve the same 1/2^n precision. The choice depends on hardware capabilities and qubit budgets.

Common Mistakes

Even experienced quantum programmers encounter pitfalls with QPE. Here are the most frequent errors and how to avoid them.

1. Preparing the Wrong Eigenstate

QPE assumes the target register holds an eigenstate of U. If you prepare a superposition of eigenstates, QPE returns a superposition of the corresponding phases. For example, if the target is (|phi_1> + |phi_2>)/sqrt(2), measurement randomly gives the phase of |phi_1> or |phi_2> with equal probability. This is sometimes useful (as in Shor’s algorithm) but is a common source of confusion when debugging.

2. Off-by-One in Binary Fraction Conversion

The measured integer m with n counting qubits gives phase phi = m/2^n, not m/2^(n+1) or m/(2^n - 1). A frequent mistake is dividing by the wrong power of 2. Always verify: with n = 3, measuring m = 1 gives phi = 1/8, not 1/16 or 1/7.

3. QFT Direction (Forward vs. Inverse)

QPE requires the inverse QFT, not the forward QFT. Applying the forward QFT instead gives the conjugate phase (1 - phi instead of phi, or equivalently, negated rotation angles). In Q#, if you define a forward QFT operation, use Adjoint QFT(register) in your QPE circuit. Swapping the direction is a subtle error that produces results that look plausible (they are valid phases) but are incorrect.

4. Insufficient Counting Qubits

For a desired precision epsilon, you need at least n > log2(1/epsilon) counting qubits. Common mistakes include using n = 3 (precision 1/8) when you need precision 0.01 (requiring n >= 7) or not accounting for the probabilistic nature of QPE when the phase is not on the grid. For high-confidence estimates, add 1-2 extra counting qubits beyond the minimum.

5. Forgetting to Reset Qubits

Q# requires that qubits be in the |0> state before deallocation. Forgetting Reset or ResetAll after measurement causes a runtime error. Always reset all qubits at the end of your operation:

// Correct: reset after measurement
ResetAll(countingRegister);
Reset(target);

// Wrong: forgetting reset causes runtime error
// (qubits may be in |1> after measurement)

6. Incorrect Power Computation

When implementing ControlledTPower, the power of 2 should use the bitshift operator <<<, not multiplication by 2 in a loop. In Q#:

// Correct: use bitshift for 2^k
let scaleFactor = 1 <<< power;

// Wrong: multiplying in a loop is error-prone and slow
mutable scaleFactor = 1;
for _ in 1..power {
    set scaleFactor = scaleFactor * 2;  // Works but fragile
}

The bitshift approach is both clearer and avoids potential integer overflow issues for large powers.

7. Confusing R1 and Rz

In Q#, R1(theta, qubit) applies diag(1, e^(itheta)), while Rz(theta, qubit) applies diag(e^(-itheta/2), e^(itheta/2)). These differ by a global phase of e^(-itheta/2). For controlled operations, global phase becomes relative phase, so Controlled R1 and Controlled Rz produce different results. QPE implementations should consistently use R1 for phase rotations to avoid sign errors.

Practical Considerations

QPE requires implementing controlled-U^(2^k) efficiently. For simple gates like T, this is straightforward because T^(2^k) is just another phase rotation. For more complex unitaries (such as those arising in chemistry simulations), decomposing controlled powers into elementary gates can be expensive and is often the bottleneck in practical applications.

The inverse QFT itself uses O(n^2) gates for n counting qubits, but approximate versions using only O(n log n) gates exist and introduce negligible error. The approximate QFT drops controlled rotations with angles below a threshold, since these contribute exponentially small phases that are masked by measurement uncertainty.

On near-term hardware, QPE circuits can be deep and sensitive to noise. Iterative Phase Estimation (IPE) offers an alternative that uses only a single counting qubit, repeating the circuit multiple times with classical feedback to extract one bit of phase per iteration. Q# supports this pattern through its classical control flow, as demonstrated in the IPE section above.

Next Steps

After mastering QPE, consider exploring these related topics:

• Variational Quantum Eigensolver (VQE): A hybrid classical-quantum algorithm that estimates ground state energies without the deep circuits required by QPE. VQE trades circuit depth for more classical optimization iterations.
• Quantum Amplitude Estimation (QAE): Uses QPE as a subroutine to estimate the amplitude of a marked state, providing a quadratic speedup over classical Monte Carlo methods.
• Robust Phase Estimation: Modified versions of QPE that use Bayesian inference or signal processing techniques to extract phases with fewer measurements.
• Hamiltonian Simulation: The techniques for implementing e^(-iHt) are the enabling technology for chemistry QPE and have applications beyond eigenvalue estimation.