Setup

The input consists of two registers (namely, two parts): the upper ${\displaystyle n}$ qubits comprise the first register, and the lower ${\displaystyle m}$ qubits are the second register.

The initial state of the system is:

${\displaystyle |0\rangle ^{\otimes n}|\psi \rangle .}$

After applying n-bit Hadamard gate operation ${\displaystyle H^{\otimes n}}$ on the first register, the state becomes:

${\displaystyle {\frac {1}{2^{\frac {n}{2}}}}(|0\rangle +|1\rangle )^{\otimes n}|\psi \rangle ={\frac {1}{2^{n/2}}}\sum _{j=0}^{2^{n}-1}|j\rangle |\psi \rangle }$.

Let ${\displaystyle U}$ be a unitary operator with eigenvector ${\displaystyle |\psi \rangle }$ such that ${\displaystyle U|\psi \rangle =e^{2\pi i\theta }|\psi \rangle }$. Thus,

${\displaystyle U^{2^{j}}|\psi \rangle =e^{2\pi i2^{j}\theta }|\psi \rangle }$.

Overall, the transformation implemented on the two registers by the controlled gates applying ${\displaystyle U,U^{2},U^{2^{2}},\ldots ,U^{2^{n-1}}}$ is

$${\displaystyle |k\rangle |\psi \rangle \mapsto |k\rangle U^{k}|\psi \rangle }$$

This can be seen by the decomposition of ${\displaystyle k}$ into its bitstring ${\displaystyle k_{n-1}k_{n-2}\ldots k_{1}k_{0}}$ and binary representation ${\displaystyle 2^{n-1}k_{n-1}+2^{n-2}k_{n-2}+\ldots +2k_{1}+k_{0}}$, where ${\displaystyle k_{n-1},\ldots ,k_{1},k_{0}\in \{0,1\}}$. Clearly, ${\displaystyle U^{k}}$ becomes

$${\displaystyle U^{2^{n-1}k_{n-1}+\ldots +2^{2}k_{2}+2k_{1}+k_{0}}=U^{2^{n-1}k_{n-1}}\ldots U^{2^{2}k_{2}}U^{2k_{1}}U^{k_{0}}}$$

Each ${\displaystyle U^{2^{j}k_{j}}}$ will only apply if the qubit ${\displaystyle k_{j}}$ is ${\displaystyle 1}$, implying that it is controlled by that bit. Therefore the overall transformation to ${\displaystyle |k\rangle U^{k}|\psi \rangle }$ is equivalent to the controlled ${\displaystyle U^{2^{j}}}$ gates from each ${\displaystyle j}$-th qubit. Therefore, the state will be transformed by the controlled ${\displaystyle U^{2^{j}}}$ gates like so:

$${\displaystyle {\frac {1}{\sqrt {2^{n}}}}\sum _{j=0}^{2^{n}-1}|j\rangle |\psi \rangle \mapsto {\frac {1}{\sqrt {2^{n}}}}\sum _{j=0}^{2^{n}-1}e^{2\pi ij\theta }|j\rangle |\psi \rangle }$$

At this point, the second register with the eigenvector is not needed. It can be reused again in another run of phase estimation. The state without ${\displaystyle |\psi \rangle }$ is

${\displaystyle {\frac {1}{\sqrt {2^{n}}}}\sum _{j=0}^{2^{n}-1}e^{2\pi ij\theta }|j\rangle }$

Apply inverse quantum Fourier transform

Applying the inverse quantum Fourier transform on

${\displaystyle {\frac {1}{2^{\frac {n}{2}}}}\sum _{k=0}^{2^{n}-1}e^{2\pi i\theta k}|k\rangle }$

yields

${\displaystyle {\frac {1}{2^{\frac {n}{2}}}}\sum _{k=0}^{2^{n}-1}e^{2\pi i\theta k}\left({\frac {1}{2^{\frac {n}{2}}}}\sum _{x=0}^{2^{n}-1}e^{\frac {-2\pi ikx}{2^{n}}}|x\rangle \right)={\frac {1}{2^{n}}}\sum _{x=0}^{2^{n}-1}\sum _{k=0}^{2^{n}-1}e^{-{\frac {2\pi ik}{2^{n}}}\left(x-2^{n}\theta \right)}|x\rangle .}$

We can approximate the value of ${\displaystyle \theta \in [0,1]}$ by rounding ${\displaystyle 2^{n}\theta }$ to the nearest integer. This means that ${\displaystyle 2^{n}\theta =a+2^{n}\delta ,}$ where ${\displaystyle a}$ is the nearest integer to ${\displaystyle 2^{n}\theta ,}$ and the difference ${\displaystyle 2^{n}\delta }$ satisfies ${\displaystyle 0\leqslant |2^{n}\delta |\leqslant {\tfrac {1}{2}}}$.

Using this decomposition we can rewrite the state as ${\textstyle \sum _{x=0}^{2^{n}-1}c_{x}|x\rangle ,}$ where

${\displaystyle c_{x}\equiv {\frac {1}{2^{n}}}\sum _{k=0}^{2^{n}-1}e^{-{\frac {2\pi ik}{2^{n}}}(x-2^{n}\theta )}={\frac {1}{2^{n}}}\sum _{k=0}^{2^{n}-1}e^{-{\frac {2\pi ik}{2^{n}}}\left(x-a\right)}e^{2\pi i\delta k}.}$

Measurement

Performing a measurement in the computational basis on the first register yields the outcome ${\displaystyle |y\rangle }$ with probability

$${\displaystyle \Pr(y)=|c_{y}|^{2}=\left|{\frac {1}{2^{n}}}\sum _{k=0}^{2^{n}-1}e^{{\frac {-2\pi ik}{2^{n}}}(y-a)}e^{2\pi i\delta k}\right|^{2}.}$$

It follows that ${\displaystyle \operatorname {Pr} (a)=1}$ if ${\displaystyle \delta =0}$, that is, when ${\displaystyle \theta }$ can be written as ${\displaystyle \theta =a/2^{n}}$, one always finds the outcome ${\displaystyle y=a}$. On the other hand, if ${\displaystyle \delta \neq 0}$, the probability reads

$${\displaystyle \operatorname {Pr} (a)={\frac {1}{2^{2n}}}\left|\sum _{k=0}^{2^{n}-1}e^{2\pi i\delta k}\right|^{2}={\frac {1}{2^{2n}}}\left|{\frac {1-{e^{2\pi i2^{n}\delta }}}{1-{e^{2\pi i\delta }}}}\right|^{2}.}$$

From this expression we can see that ${\displaystyle \Pr(a)\geqslant {\frac {4}{\pi ^{2}}}\approx 0.405}$ when ${\displaystyle \delta \neq 0}$. To see this, we observe that from the definition of ${\displaystyle \delta }$ we have the inequality ${\displaystyle |\delta |\leqslant {\tfrac {1}{2^{n+1}}}}$, and thus:^{[4]: 157 [5]: 348}

$${\displaystyle {\begin{aligned}\Pr(a)&={\frac {1}{2^{2n}}}\left|{\frac {1-{e^{2\pi i2^{n}\delta }}}{1-{e^{2\pi i\delta }}}}\right|^{2}&&{\text{for }}\delta \neq 0\\[6pt]&={\frac {1}{2^{2n}}}\left|{\frac {2\sin \left(\pi 2^{n}\delta \right)}{2\sin(\pi \delta )}}\right|^{2}&&\left|1-e^{2ix}\right|^{2}=4\left|\sin(x)\right|^{2}\\[6pt]&={\frac {1}{2^{2n}}}{\frac {\left|\sin \left(\pi 2^{n}\delta \right)\right|^{2}}{|\sin(\pi \delta )|^{2}}}\\[6pt]&\geqslant {\frac {1}{2^{2n}}}{\frac {\left|\sin \left(\pi 2^{n}\delta \right)\right|^{2}}{|\pi \delta |^{2}}}&&|\sin(\pi \delta )|\leqslant |\pi \delta |\\[6pt]&\geqslant {\frac {1}{2^{2n}}}{\frac {|2\cdot 2^{n}\delta |^{2}}{|\pi \delta |^{2}}}&&|2\cdot 2^{n}\delta |\leqslant |\sin(\pi 2^{n}\delta )|{\text{ for }}|\delta |\leqslant {\frac {1}{2^{n+1}}}\\[6pt]&\geqslant {\frac {4}{\pi ^{2}}}.\end{aligned}}}$$

We conclude that the algorithm provides the best ${\displaystyle n}$-bit estimate (i.e., one that is within ${\displaystyle 1/2^{n}}$ of the correct answer) of ${\displaystyle \theta }$ with probability at least ${\displaystyle 4/\pi ^{2}}$. By adding a number of extra qubits on the order of ${\displaystyle O(\log(1/\epsilon ))}$ and truncating the extra qubits the probability can increase to ${\displaystyle 1-\epsilon }$.^[5]