In this section, it will be shown that the problem of finding the order \(r\) satisfying \(x^r=1(mod\hspace{1mm} N)\) for some coprime integers \(x\) and \(N\), can be solved through classical means provided with random samples of integer multiples \(k/r\) of \(1/r\). Furthermore, a quantum algorithm will be constructed based off of the eigenvalue estimation algorithm which can efficiently sample estimates of these random integer multiples of \(1/r\). The results obtained from these sub-routines can then be used to factor \(N\) via the reductions described in the previous section.

Let \(r\) be some positive integer and consider the set \(\{k/r\hspace{1mm}|\hspace{1mm}k\in\mathbb{Z}\}\) consisting of integer multiples of \(1/r\). Now let \(k_1/r\) and \(k_2/r\) be two random elements from this set, and suppose the integer \(r\) is unknown. Given the values of the \(k_i/r\) it may not be possible to determine the value of \(r\) because the fractions may or may not be simplified into lower terms leaving a certain ambiguity in determining \(r\). However, if two other fractions can be found of the form \(c_1/r_1\) and \(c_2/r_2\) such that \(gcm(c_1,r_1)=gcm(c_2,r_2)=1\) with \(c_1/r_1=k_1/r\) and \(c_2/r_2=k_2/r\), then it is possible to determine the value of \(r\). With these conditions, the fractions \(c_1/r_1\) and \(c_2/r_2\) are both simplified to lowest terms so \(r_1\) and \(r_2\) can be uniquely determined. This would in turn determine \(r\) since \(lcm(r_1,r_2)=r\), where \(lcm(r_1,r_2)\) denotes the least common multiple of \(r_1\) and \(r_2\).

The classical

The problem remains to have the ability to randomly sample integer multiples \(k/r\) of \(1/r\). Actually, as a consequence of the continued fractions algorithm it will suffice to obtain an estimate of \(k/r\) to a certain degree of precision in order to arrive at the fraction \(c_1/r_1\). This motivates the following problem.

The distance constraint required in the problem statement ensures that the estimate meets a certain degree of precision in order to yield appropriate results. In addition, there is a probability constrained present to make sure such a sampling is faithful in practice.

To aid in solving the sampling problem, define the operator given by

\[U_x: \left|s\right>\mapsto \left|sx (mod\hspace{1mm}N)\right>.\]

with the assumption that \(x\) is coprime to \(N\) such a map always has a well-defined inverse so that \(U_x\) in indeed a unitary operator. The order \(r\) of \(x\) is reflected in the action of the operator \(U_x\) since \(x^r=1(mod\hspace{1mm} N)\) and

\[U_x^r\left|s\right>=\left|sx^r (mod\hspace{1mm}N)\right>=\left|s\right>.\]

This implies that the eigenvalues of \(\left|s\right>\) will actually be \(r^{th}\) roots of unity of the form \(e^{i2\pi k/r}\) for integers \(k\in\{0,1,\dots, r-1\}\).

Now, consider the state

\[\left|u_k\right>=\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\left|x^s (mod\hspace{1mm}N)\right>.\]

Then observe that

\[\begin{array}{r l}

U_x\left|u_k\right>=&\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}U_x\left|x^s (mod\hspace{1mm}N)\right> \\

=& \displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\left|x^{s+1} (mod\hspace{1mm}N)\right> \\

=& \displaystyle\frac{e^{i2\pi \frac{k}{r}}}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} (s+1)}\left|x^{s+1} (mod\hspace{1mm}N)\right>\\

=&e^{i2\pi \frac{k}{r}}\left|u_k\right>,

\end{array}\]

where the last identity follows from \[e^{i2\pi \frac{k}{r}r}\left|x^r (mod\hspace{1mm}N)\right>=e^{i2\pi \frac{k}{r}0}\left|x^0 (mod\hspace{1mm}N)\right>=\left|1 (mod\hspace{1mm}N)\right>.\] Hence, the \(\left|u_k\right>\) is an eigenstate of \(U_x\) with eigenvalue \(e^{i2\pi \frac{k}{r}}\).

If such a state \(\left|u_k\right>\) was provided, then the eigenvalue estimation algorithm can be implemented on the input state \(\left|0\right>\left|u_k\right>\) to yield the output \(\left|\tilde{k/r}\right>\left|u_k\right>\). Then upon measuring the first register an estimate \(\tilde{k/r}\) of some integer multiple of \(1/r\) could be determined. The issue here is that, since \(r\) is originally unknown, the state \(\left|u_k\right>\) cannot be prepared directly because it is defined in terms of \(r\). Fortunately, such a state does not need to be constructed directly in order to effectively sample estimate to an integer multiple of \(1/r\). For consider the following equally weighted superposition of all \(r\) states \(\left|u_k\right>\)

\[\begin{array}{r l}

\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|u_k\right>=& \displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\left|x^s (mod\hspace{1mm}N)\right> \\

=& \displaystyle\frac{1}{r}\displaystyle\sum\limits_{s=0}^{r-1}\left(\displaystyle\sum\limits_{k=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\right)\left|x^s (mod\hspace{1mm}N)\right>

\end{array}\]

The amplitude \(\alpha_1\) of the state \(\left|1\right>=\left|x^0 (mod\hspace{1mm}N)\right>\) in the above superposition. Notice this occurs when \(s=0\) since \(x^0=1\). So the amplitude is therefore

\[\alpha_1=\displaystyle\frac{1}{r}\displaystyle\sum\limits_{k=0}^{r-1}e^{-i2\pi \frac{k}{r} 0}=\displaystyle\frac{1}{r}\displaystyle\sum\limits_{k=0}^{r-1}1=1.\]

Since the superposition satisfies the normalization constraint, the amplitudes of all other states \(\alpha_s\) for \(s \neq 0\) must be zero. Hence,

\[\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|u_k\right>=\left|1\right>\] and this superposition is really just the computational basis state \(\left|1\right>\) in disguise.

If instead this superposition was input into the eigenvalue estimation algorithm, then some superposition of states corresponding to estimates of integer multiples of \(1/r\) will be output. This is ultimately a consequence of the linearity of the operator implementing the eigenvalue estimation algorithm. More explicitly, the action of the eigenvalue algorithm when the target register is prepared in the superposition consisting of equally weighted eigenstates, outputs an equally weighted superposition of estimates of integer multiples of \(1/r\):

\[\left|0\right>\left|1\right>=\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|0\right>\left|u_k\right> \mapsto \displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|\tilde{k/r}\right>\left|u_k\right>.\]

By measuring the control register some state \(\left|\tilde{k/r}\right>\) corresponding to an integer \(x\) will be observed with equal probability such that \(x/2^n\) is an estimate of an integer multiple \(\frac{k}{r}\). Therefore, a uniformly random integer multiple of \(1/r\) can be sampled. A circuit depicting the algorithm uses to sample estimates of a random integer multiple of \(1/r\) is shown in the figure below. Having found such samples, allows the order \(r\) to be determined using the continued fractions algorithm. Then, this order \(r\) can be used to find a non-trivial factor of the integer \(N\) as described through the classical methods involved in the reduction of the problem.

(A circuit for the problem of sampling random integer multiples of \(1/r\). This circuit is analogous to the general circuit for the eigenvalue estimation problem shown here, except here the controlled \(-U_{x}^{x}\) operation is defined differently. In this circuit, the input in the target register is prepared to be in the state \(\left|1\right>\) instead of an eigenstate of \(U_x\). Consequently the output in the control register is a superposition of estimates \(\left|\tilde{k/r}\right>\) of integer multiples of \(1/r\).)

This circuit is nearly identical to the circuit presented for the general eigenvalue estimation algorithm. Here, the controlled\(-U_{x}^{x}\) gate is supposed to signify an analogous sequence of controlled \(-U_{x}^{2^s}\) operations, where \(s\) ranges from \(0\leq s \leq r-1\). Since the action of the \(U_x\) gate is defined to multiply a basis state by \(x(mod\hspace{1mm}N)\), a natural way to compute the \(c-U_x^{2^s}\) gate for some \(s\) would be to simply apply the \(c-U_x\) gate \(2^s\) times. This, however, would not be efficient as it requires an exponential number of operations to implement. Instead, observe that by applying a similarly defined gate \(c-U_{x^{2^s}}\) accomplishes an equivalent operation and only needs to be applied a single time. Moreover, calculating the value of \(x^{2^s}\) by repeatedly squaring \(x^2\) modulo \(N\) can be done classically with only \(s\)-many multiplications.

To be more rigorous, it can be shown that the \(c-U_{x}^{x}\) operation can be computed with time in \(O((logN)^2loglogNlogloglogN)\).\cite{mosca} Recall that the \(QFT\) and \(QFT^{-1}\) can each be computed in time \(O((logN)^2)\). Therefore, the whole sampling problem can be implemented efficiently in this case with time \(O((logN)^2loglogNlogloglogN)\). When compared to the classical complexity of the same problem the best known bounds are in \(e^{O((logN)^{1/2}(loglogN)^{1/2})}\). Thus, the quantum algorithm can perform exponentially faster than the classical when trying to factor integers.

It is worth mentioning that the approach taken by Shor in his factoring algorithm \cite{shor} involved a slightly different analysis. However, when comparing his method to the one outlined here it can be seen that the main different results from analyzing the problem in a different basis.

Let \(r\) be some positive integer and consider the set \(\{k/r\hspace{1mm}|\hspace{1mm}k\in\mathbb{Z}\}\) consisting of integer multiples of \(1/r\). Now let \(k_1/r\) and \(k_2/r\) be two random elements from this set, and suppose the integer \(r\) is unknown. Given the values of the \(k_i/r\) it may not be possible to determine the value of \(r\) because the fractions may or may not be simplified into lower terms leaving a certain ambiguity in determining \(r\). However, if two other fractions can be found of the form \(c_1/r_1\) and \(c_2/r_2\) such that \(gcm(c_1,r_1)=gcm(c_2,r_2)=1\) with \(c_1/r_1=k_1/r\) and \(c_2/r_2=k_2/r\), then it is possible to determine the value of \(r\). With these conditions, the fractions \(c_1/r_1\) and \(c_2/r_2\) are both simplified to lowest terms so \(r_1\) and \(r_2\) can be uniquely determined. This would in turn determine \(r\) since \(lcm(r_1,r_2)=r\), where \(lcm(r_1,r_2)\) denotes the least common multiple of \(r_1\) and \(r_2\).

The classical

*continued fractions algorithm*\cite{hardy} (although not discussed here for the sake of brevity), can be used to find fractions of the type \(c_1/r_1\) and \(c_2/r_2\) provided with some random multiples \(k_1/r\) and \(k_2/r\), and can do so efficiently with a time complexity in \(O(log^2r)\). Making use of this technique will be essential in the post-processing required to determine orders after the quantum part of the algorithm is applied.The problem remains to have the ability to randomly sample integer multiples \(k/r\) of \(1/r\). Actually, as a consequence of the continued fractions algorithm it will suffice to obtain an estimate of \(k/r\) to a certain degree of precision in order to arrive at the fraction \(c_1/r_1\). This motivates the following problem.

**Sampling Estimates to random integer multiples of 1/r****Input:**Integers \(x\) and \(N\) such that \(gcd(x,N)=1\)**Problem:**Output a number \(y\in\{0,1,\dots,2^{n}-1\}\) such that for each \(k\in\{0,1,\dots,r-1\}\), the probability \(p\) such that \(\left| \frac{y}{2^n}- \frac{k}{r}\right|\leq \frac{1}{2r^2}\) is \(p\geq\frac{c}{r}\) for some constant \(c>0\), where \(r\) denotes the unknown order of \(x\).The distance constraint required in the problem statement ensures that the estimate meets a certain degree of precision in order to yield appropriate results. In addition, there is a probability constrained present to make sure such a sampling is faithful in practice.

To aid in solving the sampling problem, define the operator given by

\[U_x: \left|s\right>\mapsto \left|sx (mod\hspace{1mm}N)\right>.\]

with the assumption that \(x\) is coprime to \(N\) such a map always has a well-defined inverse so that \(U_x\) in indeed a unitary operator. The order \(r\) of \(x\) is reflected in the action of the operator \(U_x\) since \(x^r=1(mod\hspace{1mm} N)\) and

\[U_x^r\left|s\right>=\left|sx^r (mod\hspace{1mm}N)\right>=\left|s\right>.\]

This implies that the eigenvalues of \(\left|s\right>\) will actually be \(r^{th}\) roots of unity of the form \(e^{i2\pi k/r}\) for integers \(k\in\{0,1,\dots, r-1\}\).

Now, consider the state

\[\left|u_k\right>=\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\left|x^s (mod\hspace{1mm}N)\right>.\]

Then observe that

\[\begin{array}{r l}

U_x\left|u_k\right>=&\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}U_x\left|x^s (mod\hspace{1mm}N)\right> \\

=& \displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\left|x^{s+1} (mod\hspace{1mm}N)\right> \\

=& \displaystyle\frac{e^{i2\pi \frac{k}{r}}}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} (s+1)}\left|x^{s+1} (mod\hspace{1mm}N)\right>\\

=&e^{i2\pi \frac{k}{r}}\left|u_k\right>,

\end{array}\]

where the last identity follows from \[e^{i2\pi \frac{k}{r}r}\left|x^r (mod\hspace{1mm}N)\right>=e^{i2\pi \frac{k}{r}0}\left|x^0 (mod\hspace{1mm}N)\right>=\left|1 (mod\hspace{1mm}N)\right>.\] Hence, the \(\left|u_k\right>\) is an eigenstate of \(U_x\) with eigenvalue \(e^{i2\pi \frac{k}{r}}\).

If such a state \(\left|u_k\right>\) was provided, then the eigenvalue estimation algorithm can be implemented on the input state \(\left|0\right>\left|u_k\right>\) to yield the output \(\left|\tilde{k/r}\right>\left|u_k\right>\). Then upon measuring the first register an estimate \(\tilde{k/r}\) of some integer multiple of \(1/r\) could be determined. The issue here is that, since \(r\) is originally unknown, the state \(\left|u_k\right>\) cannot be prepared directly because it is defined in terms of \(r\). Fortunately, such a state does not need to be constructed directly in order to effectively sample estimate to an integer multiple of \(1/r\). For consider the following equally weighted superposition of all \(r\) states \(\left|u_k\right>\)

\[\begin{array}{r l}

\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|u_k\right>=& \displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{s=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\left|x^s (mod\hspace{1mm}N)\right> \\

=& \displaystyle\frac{1}{r}\displaystyle\sum\limits_{s=0}^{r-1}\left(\displaystyle\sum\limits_{k=0}^{r-1}e^{-i2\pi \frac{k}{r} s}\right)\left|x^s (mod\hspace{1mm}N)\right>

\end{array}\]

The amplitude \(\alpha_1\) of the state \(\left|1\right>=\left|x^0 (mod\hspace{1mm}N)\right>\) in the above superposition. Notice this occurs when \(s=0\) since \(x^0=1\). So the amplitude is therefore

\[\alpha_1=\displaystyle\frac{1}{r}\displaystyle\sum\limits_{k=0}^{r-1}e^{-i2\pi \frac{k}{r} 0}=\displaystyle\frac{1}{r}\displaystyle\sum\limits_{k=0}^{r-1}1=1.\]

Since the superposition satisfies the normalization constraint, the amplitudes of all other states \(\alpha_s\) for \(s \neq 0\) must be zero. Hence,

\[\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|u_k\right>=\left|1\right>\] and this superposition is really just the computational basis state \(\left|1\right>\) in disguise.

If instead this superposition was input into the eigenvalue estimation algorithm, then some superposition of states corresponding to estimates of integer multiples of \(1/r\) will be output. This is ultimately a consequence of the linearity of the operator implementing the eigenvalue estimation algorithm. More explicitly, the action of the eigenvalue algorithm when the target register is prepared in the superposition consisting of equally weighted eigenstates, outputs an equally weighted superposition of estimates of integer multiples of \(1/r\):

\[\left|0\right>\left|1\right>=\displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|0\right>\left|u_k\right> \mapsto \displaystyle\frac{1}{\sqrt{r}}\displaystyle\sum\limits_{k=0}^{r-1}\left|\tilde{k/r}\right>\left|u_k\right>.\]

By measuring the control register some state \(\left|\tilde{k/r}\right>\) corresponding to an integer \(x\) will be observed with equal probability such that \(x/2^n\) is an estimate of an integer multiple \(\frac{k}{r}\). Therefore, a uniformly random integer multiple of \(1/r\) can be sampled. A circuit depicting the algorithm uses to sample estimates of a random integer multiple of \(1/r\) is shown in the figure below. Having found such samples, allows the order \(r\) to be determined using the continued fractions algorithm. Then, this order \(r\) can be used to find a non-trivial factor of the integer \(N\) as described through the classical methods involved in the reduction of the problem.

(A circuit for the problem of sampling random integer multiples of \(1/r\). This circuit is analogous to the general circuit for the eigenvalue estimation problem shown here, except here the controlled \(-U_{x}^{x}\) operation is defined differently. In this circuit, the input in the target register is prepared to be in the state \(\left|1\right>\) instead of an eigenstate of \(U_x\). Consequently the output in the control register is a superposition of estimates \(\left|\tilde{k/r}\right>\) of integer multiples of \(1/r\).)

This circuit is nearly identical to the circuit presented for the general eigenvalue estimation algorithm. Here, the controlled\(-U_{x}^{x}\) gate is supposed to signify an analogous sequence of controlled \(-U_{x}^{2^s}\) operations, where \(s\) ranges from \(0\leq s \leq r-1\). Since the action of the \(U_x\) gate is defined to multiply a basis state by \(x(mod\hspace{1mm}N)\), a natural way to compute the \(c-U_x^{2^s}\) gate for some \(s\) would be to simply apply the \(c-U_x\) gate \(2^s\) times. This, however, would not be efficient as it requires an exponential number of operations to implement. Instead, observe that by applying a similarly defined gate \(c-U_{x^{2^s}}\) accomplishes an equivalent operation and only needs to be applied a single time. Moreover, calculating the value of \(x^{2^s}\) by repeatedly squaring \(x^2\) modulo \(N\) can be done classically with only \(s\)-many multiplications.

To be more rigorous, it can be shown that the \(c-U_{x}^{x}\) operation can be computed with time in \(O((logN)^2loglogNlogloglogN)\).\cite{mosca} Recall that the \(QFT\) and \(QFT^{-1}\) can each be computed in time \(O((logN)^2)\). Therefore, the whole sampling problem can be implemented efficiently in this case with time \(O((logN)^2loglogNlogloglogN)\). When compared to the classical complexity of the same problem the best known bounds are in \(e^{O((logN)^{1/2}(loglogN)^{1/2})}\). Thus, the quantum algorithm can perform exponentially faster than the classical when trying to factor integers.

It is worth mentioning that the approach taken by Shor in his factoring algorithm \cite{shor} involved a slightly different analysis. However, when comparing his method to the one outlined here it can be seen that the main different results from analyzing the problem in a different basis.