El algoritmo de Shor

Ahora dirigiremos nuestra atención al problema de la factorización de enteros y veremos cómo puede resolverse eficientemente en una computadora cuántica mediante la estimación de fase. El algoritmo que obtendremos es el algoritmo de Shor para la factorización de enteros. Shor no describió su algoritmo específicamente en términos de estimación de fase, pero es una forma natural e intuitiva de explicar cómo funciona.

Comenzaremos analizando un problema intermedio conocido como el problema de búsqueda de orden, y veremos cómo la estimación de fase proporciona una solución a este problema. Luego veremos cómo una solución eficiente al problema de búsqueda de orden nos da una solución eficiente al problema de factorización de enteros. (Cuando la solución de un problema proporciona la solución de otro problema de este modo, decimos que el segundo problema se reduce al primero — así que en este caso estamos reduciendo la factorización de enteros a la búsqueda de orden.) Esta segunda parte del algoritmo de Shor no hace uso de la computación cuántica en absoluto; es completamente clásica. La computación cuántica solo es necesaria para resolver la búsqueda de orden.

El problema de búsqueda de orden

Algo de teoría básica de números

Para explicar el problema de búsqueda de orden y cómo puede resolverse mediante estimación de fase, será útil comenzar con un par de conceptos básicos de teoría de números, e introducir algo de notación conveniente en el camino.

Para empezar, para cualquier entero positivo $N$ dado, definimos el conjunto $\mathbb{Z}_N$ de la siguiente manera.

$$\mathbb{Z}_N = \{0,1,\ldots,N-1\}$$

Por ejemplo, $\mathbb{Z}_1 = \{0\},\;$ $\mathbb{Z}_2 = \{0,1\},\;$ $\mathbb{Z}_3 = \{0,1,2\},\;$ y así sucesivamente.

Estos son conjuntos de números, pero podemos pensar en ellos como algo más que simples conjuntos. En particular, podemos pensar en operaciones aritméticas sobre $\mathbb{Z}_N$, como la suma y la multiplicación — y si acordamos tomar siempre nuestros resultados módulo $N$ (es decir, dividir entre $N$ y tomar el resto como resultado), siempre permaneceremos dentro de este conjunto al realizar estas operaciones. Las dos operaciones específicas de suma y multiplicación, ambas tomadas módulo $N$, convierten $\mathbb{Z}_N$ en un anillo, que es un tipo de objeto fundamentalmente importante en álgebra.

Por ejemplo, $3$ y $5$ son elementos de $\mathbb{Z}_7$, y si los multiplicamos obtenemos $3\cdot 5 = 15$, que deja un resto de $1$ al dividirlo entre $7$. A veces esto se expresa de la siguiente forma.

$$3 \cdot 5 \equiv 1 \; (\textrm{mod } 7)$$

Pero también podemos simplemente escribir $3 \cdot 5 = 1$, siempre que haya quedado claro que estamos trabajando en $\mathbb{Z}_7$, para mantener la notación lo más sencilla posible.

Como ejemplo, aquí están las tablas de suma y multiplicación para $\mathbb{Z}_6.$

$$\begin{array}{c|cccccc} + & 0 & 1 & 2 & 3 & 4 & 5 \\\hline 0 & 0 & 1 & 2 & 3 & 4 & 5 \\ 1 & 1 & 2 & 3 & 4 & 5 & 0 \\ 2 & 2 & 3 & 4 & 5 & 0 & 1 \\ 3 & 3 & 4 & 5 & 0 & 1 & 2 \\ 4 & 4 & 5 & 0 & 1 & 2 & 3 \\ 5 & 5 & 0 & 1 & 2 & 3 & 4 \\ \end{array} \qquad \begin{array}{c|cccccc} \cdot & 0 & 1 & 2 & 3 & 4 & 5 \\\hline 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 2 & 3 & 4 & 5 \\ 2 & 0 & 2 & 4 & 0 & 2 & 4 \\ 3 & 0 & 3 & 0 & 3 & 0 & 3 \\ 4 & 0 & 4 & 2 & 0 & 4 & 2 \\ 5 & 0 & 5 & 4 & 3 & 2 & 1 \\ \end{array}$$

Entre los $N$ elementos de $\mathbb{Z}_N$, los elementos $a\in\mathbb{Z}_N$ que satisfacen $\gcd(a,N) = 1$ son especiales. Frecuentemente, el conjunto que contiene estos elementos se denota con un asterisco de la siguiente manera.

$$\mathbb{Z}_N^{\ast} = \{a\in \mathbb{Z}_N : \gcd(a,N) = 1\}$$

Si centramos nuestra atención en la operación de multiplicación, el conjunto $\mathbb{Z}_N^{\ast}$ forma un grupo — específicamente un grupo abeliano — que es otro tipo importante de objeto en álgebra. Es un hecho básico sobre estos conjuntos (y los grupos finitos en general) que si tomamos cualquier elemento $a\in\mathbb{Z}_N^{\ast}$ y multiplicamos $a$ por sí mismo repetidamente, siempre acabaremos obteniendo el número $1$.

Como primer ejemplo, tomemos $N=6$. Tenemos que $5\in\mathbb{Z}_6^{\ast}$ porque $\gcd(5,6) = 1$, y si multiplicamos $5$ por sí mismo obtenemos $1$, como confirma la tabla anterior.

$$5^2 = 1 \quad \text{(trabajando dentro de $\mathbb{Z}_6$)}$$

Como segundo ejemplo, tomemos $N = 21$. Si recorremos los números del $0$ al $20$, los que tienen MCD igual a $1$ con $21$ son los siguientes.

$$\mathbb{Z}_{21}^{\ast} = \{1,2,4,5,8,10,11,13,16,17,19,20\}$$

Para cada uno de estos elementos, es posible elevar ese número a una potencia entera positiva para obtener $1$. A continuación se muestran las potencias más pequeñas para las que esto funciona:

$$\begin{array}{ccc} 1^{1} = 1 \quad & 8^{2} = 1 \quad & 16^{3} = 1 \\[1mm] 2^{6} = 1 \quad & 10^{6} = 1 \quad & 17^{6} = 1 \\[1mm] 4^{3} = 1 \quad & 11^{6} = 1 \quad & 19^{6} = 1 \\[1mm] 5^{6} = 1 \quad & 13^{2} = 1 \quad & 20^{2} = 1 \end{array}$$

Naturalmente estamos trabajando dentro de $\mathbb{Z}_{21}$ en todas estas ecuaciones, lo cual no hemos escrito explícitamente — lo damos por implícito para no recargar la notación. Continuaremos haciendo esto a lo largo del resto de la lección.

Enunciado del problema y conexión con la estimación de fase

Ahora podemos enunciar el problema de búsqueda de orden.

Búsqueda de orden

Entrada: enteros positivos $N$ y $a$ que satisfacen $\gcd(N,a) = 1$
Salida: el entero positivo más pequeño $r$ tal que $a^r \equiv 1$ $(\textrm{mod } N)$

Alternativamente, en términos de la notación que acabamos de introducir, se nos da $a \in \mathbb{Z}_N^{\ast}$, y buscamos el entero positivo más pequeño $r$ tal que $a^r = 1$. Este número $r$ se llama el orden de $a$ módulo $N$.

Para conectar el problema de búsqueda de orden con la estimación de fase, pensemos en la operación definida sobre un sistema cuyos estados clásicos corresponden a $\mathbb{Z}_N$, donde multiplicamos por un elemento fijo $a\in\mathbb{Z}_N^{\ast}$.

$$M_a \vert x\rangle = \vert ax \rangle \qquad \text{(para cada $x\in\mathbb{Z}_N$)}$$

Para ser precisos, estamos realizando la multiplicación en $\mathbb{Z}_N$, por lo que se sobreentiende que tomamos el producto módulo $N$ dentro del ket en el lado derecho de la ecuación.

Por ejemplo, si tomamos $N = 15$ y $a=2$, entonces la acción de $M_2$ sobre la base estándar $\{\vert 0\rangle,\ldots,\vert 14\rangle\}$ es la siguiente.

$$\begin{array}{ccc} M_{2} \vert 0 \rangle = \vert 0\rangle \quad & M_{2} \vert 5 \rangle = \vert 10\rangle \quad & M_{2} \vert 10 \rangle = \vert 5\rangle \\[1mm] M_{2} \vert 1 \rangle = \vert 2\rangle \quad & M_{2} \vert 6 \rangle = \vert 12\rangle \quad & M_{2} \vert 11 \rangle = \vert 7\rangle \\[1mm] M_{2} \vert 2 \rangle = \vert 4\rangle \quad & M_{2} \vert 7 \rangle = \vert 14\rangle \quad & M_{2} \vert 12 \rangle = \vert 9\rangle \\[1mm] M_{2} \vert 3 \rangle = \vert 6\rangle \quad & M_{2} \vert 8 \rangle = \vert 1\rangle \quad & M_{2} \vert 13 \rangle = \vert 11\rangle \\[1mm] M_{2} \vert 4 \rangle = \vert 8\rangle \quad & M_{2} \vert 9 \rangle = \vert 3\rangle \quad & M_{2} \vert 14 \rangle = \vert 13\rangle \end{array}$$

Esta es una operación unitaria siempre que $\gcd(a,N)=1$; permuta los elementos de la base estándar $\{\vert 0\rangle,\ldots,\vert N-1\rangle\}$, de modo que como matriz es una matriz de permutación. Es evidente por su definición que esta operación es determinista, y una forma sencilla de ver que es invertible consiste en pensar en el orden $r$ de $a$ módulo $N$, y reconocer que la inversa de $M_a$ es $M_a^{r-1}$.

$$M_a^{r-1} M_a = M_a^r = M_{a^r} = M_1 = \mathbb{I}$$

Hay otra forma de pensar en la inversa que no requiere ningún conocimiento de $r$ (que, al fin y al cabo, es lo que estamos tratando de calcular). Para todo elemento $a\in\mathbb{Z}_N^{\ast}$ siempre existe un único elemento $b\in\mathbb{Z}_N^{\ast}$ que satisface $ab=1$. Denotamos este elemento $b$ como $a^{-1}$, y puede calcularse eficientemente; una extensión del algoritmo MCD de Euclides lo hace con un coste cuadrático en $\operatorname{lg}(N)$. Y así

$$M_{a^{-1}} M_a = M_{a^{-1}a} = M_1 = \mathbb{I}.$$

Por tanto, la operación $M_a$ es a la vez determinista e invertible. Eso implica que está descrita por una matriz de permutación y, por lo tanto, es unitaria.

Ahora pensemos en los vectores propios y los valores propios de la operación $M_a$, suponiendo que $a\in\mathbb{Z}_N^{\ast}$. Como acabamos de argumentar, esta suposición nos dice que $M_a$ es unitaria.

Hay $N$ valores propios de $M_a$, posiblemente incluyendo el mismo valor propio repetido varias veces, y en general hay cierta libertad al seleccionar los vectores propios correspondientes — pero no necesitamos preocuparnos por todas las posibilidades. Empecemos de manera sencilla e identifiquemos solo un vector propio de $M_a$.

$$\vert \psi_0 \rangle = \frac{\vert 1 \rangle + \vert a \rangle + \cdots + \vert a^{r-1} \rangle}{\sqrt{r}}$$

El número $r$ es el orden de $a$ módulo $N$, aquí y en el resto de la lección. El valor propio asociado a este vector propio es $1$ porque no cambia cuando multiplicamos por $a$.

$$M_a \vert \psi_0 \rangle = \frac{\vert a \rangle + \cdots + \vert a^{r-1} \rangle + \vert a^r \rangle}{\sqrt{r}} = \frac{\vert a \rangle + \cdots + \vert a^{r-1} \rangle + \vert 1 \rangle}{\sqrt{r}} = \vert \psi_0 \rangle$$

Esto ocurre porque $a^r = 1$, de modo que cada estado de la base estándar $\vert a^k \rangle$ se desplaza a $\vert a^{k+1} \rangle$ para $k\leq r-1$, y $\vert a^{r-1} \rangle$ vuelve a desplazarse a $\vert 1\rangle$. Informalmente hablando, es como si estuviéramos agitando lentamente $\vert \psi_0 \rangle$, pero ya está completamente mezclado y nada cambia.

Aquí hay otro ejemplo de vector propio de $M_a$. Este resulta ser más interesante en el contexto de la búsqueda de orden y la estimación de fase.

$$\vert \psi_1 \rangle = \frac{\vert 1 \rangle + \omega_r^{-1} \vert a \rangle + \cdots + \omega_r^{-(r-1)}\vert a^{r-1} \rangle}{\sqrt{r}}$$

Alternativamente, podemos escribir este vector usando una suma de la siguiente manera.

$$\vert \psi_1 \rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^k \rangle$$

Aquí vemos aparecer de forma natural el número complejo $\omega_r = e^{2\pi i/r}$, debido a la manera en que funciona la multiplicación por $a$ módulo $N$. Esta vez el valor propio correspondiente es $\omega_r$. Para verlo, podemos calcular primero de la siguiente manera.

$$M_a \vert \psi_1 \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} M_a\vert a^k \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^{k+1} \rangle = \frac{1}{\sqrt{r}}\sum_{k = 1}^{r} \omega_r^{-(k - 1)} \vert a^{k} \rangle = \frac{1}{\sqrt{r}}\omega_r \sum_{k = 1}^{r} \omega_r^{-k} \vert a^{k} \rangle$$

Luego, como $\omega_r^{-r} = 1 = \omega_r^0$ y $\vert a^r \rangle = \vert 1\rangle = \vert a^0\rangle$, vemos que

$$\frac{1}{\sqrt{r}}\sum_{k = 1}^{r} \omega_r^{-k} \vert a^{k} \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^k \rangle = \vert\psi_1\rangle,$$

por lo que $M_a \vert\psi_1\rangle = \omega_r \vert\psi_1\rangle.$

Usando el mismo razonamiento, podemos identificar pares adicionales de vector propio/valor propio para $M_a$. Para cualquier elección de $j\in\{0,\ldots,r-1\}$ tenemos que

$$\vert \psi_j \rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \omega_r^{-jk} \vert a^k \rangle$$

es un vector propio de $M_a$ cuyo valor propio correspondiente es $\omega_r^j$.

$$M_a \vert \psi_j \rangle = \omega_r^j \vert \psi_j \rangle$$

Existen otros vectores propios de $M_a$, pero no necesitamos preocuparnos por ellos — nos centraremos únicamente en los vectores propios $\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle$ que acabamos de identificar.

Búsqueda de orden mediante estimación de fase

Para resolver el problema de búsqueda de orden para una elección dada de $a\in\mathbb{Z}_N^{\ast},$ podemos aplicar el procedimiento de estimación de fase a la operación $M_a.$

Para hacer esto, necesitamos implementar eficientemente con un circuito cuántico no solo $M_a,$ sino también $M_a^2,$ $M_a^4,$ $M_a^8,$ y así sucesivamente, avanzando tanto como sea necesario para obtener una estimación suficientemente precisa del procedimiento de estimación de fase. Aquí explicaremos cómo puede hacerse esto, y determinaremos exactamente cuánta precisión se necesita más adelante.

Empecemos con la operación $M_a$ por sí sola. Naturalmente, como estamos trabajando con el modelo de circuitos cuánticos, usaremos notación binaria para codificar los números entre $0$ y $N-1.$ El número más grande que necesitamos codificar es $N-1,$ por lo que el número de bits que necesitamos es

$$n = \operatorname{lg}(N-1) = \lfloor \log(N-1) \rfloor + 1.$$

Por ejemplo, si $N = 21$ tenemos $n = \operatorname{lg}(N-1) = 5.$ Así es como se ve la codificación de los elementos de $\mathbb{Z}_{21}$ como cadenas binarias de longitud $5.$

$$\begin{gathered} 0 \mapsto 00000\\[1mm] 1 \mapsto 00001\\[1mm] \vdots\\[1mm] 20 \mapsto 10100 \end{gathered}$$

Y ahora, aquí está la definición precisa de cómo se define $M_a$ como una operación de $n$ cúbits.

$$M_a \vert x\rangle = \begin{cases} \vert ax \; (\textrm{mod}\;N)\rangle & 0\leq x < N\\[1mm] \vert x\rangle & N\leq x < 2^n \end{cases}$$

El punto es que aunque solo nos importa cómo funciona $M_a$ para $\vert 0\rangle,\ldots,\vert N-1\rangle,$ sí tenemos que especificar cómo funciona para los restantes $2^n - N$ estados de la base estándar — y necesitamos hacerlo de una forma que nos siga dando una operación unitaria. Definir $M_a$ de manera que no haga nada a los estados de la base estándar restantes logra esto.

Usando los algoritmos para la multiplicación entera y la división discutidos en la lección anterior, junto con la metodología para implementaciones reversibles y libres de basura de estos, podemos construir un circuito cuántico que realice $M_a,$ para cualquier elección de $a\in\mathbb{Z}_N^{\ast},$ con un costo de $O(n^2).$ Una forma de hacer esto es la siguiente.

1. Construir un circuito para realizar la operación

   $$\vert x \rangle \vert y \rangle \mapsto \vert x \rangle \vert y \oplus f_a(x)\rangle$$

   donde

   $$f_a(x) = \begin{cases} ax \; (\textrm{mod}\;N) & 0\leq x < N\\[1mm] x & N\leq x < 2^n \end{cases}$$

   usando el método descrito en la lección anterior. Esto nos da un circuito de tamaño $O(n^2).$

2. Intercambiar los dos sistemas de $n$ cúbits usando $n$ compuertas swap para intercambiar los cúbits individualmente.

3. De forma similar al primer paso, construir un circuito para la operación

   $$\vert x \rangle \vert y \rangle \mapsto \vert x \rangle \bigl\vert y \oplus f_{a^{-1}}(x)\bigr\rangle$$

   donde $a^{-1}$ es el inverso de $a$ en $\mathbb{Z}_N^{\ast}.$

Al inicializar los $n$ cúbits inferiores y componer los tres pasos, obtenemos esta transformación:

$$\vert x \rangle \vert 0^n \rangle \stackrel{\text{paso 1}}{\mapsto} \vert x \rangle \vert f_a(x)\rangle \stackrel{\text{paso 2}}{\mapsto} \vert f_a(x)\rangle \vert x \rangle \stackrel{\text{paso 3}}{\mapsto} \vert f_a(x)\rangle \bigl\vert x \oplus f_{a^{-1}}(f_a(x)) \bigr\rangle = \vert f_a(x)\rangle\vert 0^n \rangle$$

El método requiere cúbits de espacio de trabajo, pero estos se devuelven a su estado inicializado al final, lo que nos permite usar estos circuitos para la estimación de fase. El costo total del circuito que obtenemos es $O(n^2).$

Para realizar $M_a^2,$ $M_a^4,$ $M_a^8,$ y así sucesivamente, podemos usar exactamente el mismo método, excepto que reemplazamos $a$ por $a^2,$ $a^4,$ $a^8,$ y así sucesivamente, como elementos de $\mathbb{Z}_N^{\ast}.$ Es decir, para cualquier potencia $k$ que elijamos, podemos crear un circuito para $M_a^k$ no iterando $k$ veces el circuito para $M_a,$ sino calculando $b = a^k \in \mathbb{Z}_N^{\ast}$ y luego usando el circuito para $M_b.$

El cálculo de potencias $a^k \in \mathbb{Z}_N$ es el problema de exponenciación modular mencionado en la lección anterior. Este cálculo puede realizarse clásicamente, usando el algoritmo de exponenciación modular mencionado en la lección anterior (frecuentemente llamado el algoritmo de potenciación en la teoría computacional de números). De hecho, solo necesitamos potencias de $a$ que sean potencias de 2, en particular $a^2, a^4, \ldots a^{2^{m-1}} \in \mathbb{Z}_N^{\ast},$ y podemos obtener estas potencias elevando al cuadrado iterativamente $m-1$ veces. Cada elevación al cuadrado puede realizarse mediante un circuito booleano de tamaño $O(n^2).$

En esencia, lo que estamos haciendo aquí es delegar el problema de iterar $M_a$ hasta $2^{m-1}$ veces a un cálculo clásico eficiente. ¡Y es una gran suerte que esto sea posible! Para una elección arbitraria de un circuito cuántico en el problema de estimación de fase, es poco probable que esto sea posible — y en ese caso el costo resultante para la estimación de fase crece exponencialmente en el número de cúbits de control $m.$

Solución dado un autovector conveniente

Para entender cómo podemos resolver el problema de búsqueda de orden usando estimación de fase, empecemos suponiendo que ejecutamos el procedimiento de estimación de fase sobre la operación $M_a$ usando el autovector $\vert\psi_1\rangle.$ Conseguir este autovector no es fácil, como veremos, así que esta no será la historia completa — pero es útil comenzar aquí.

El autovalor de $M_a$ correspondiente al autovector $\vert \psi_1\rangle$ es

$$\omega_r = e^{2\pi i \frac{1}{r}}.$$

Es decir, $\omega_r = e^{2\pi i \theta}$ para $\theta = 1/r.$ Entonces, si ejecutamos el procedimiento de estimación de fase sobre $M_a$ usando el autovector $\vert\psi_1\rangle,$ obtendremos una aproximación a $1/r.$ Calculando el recíproco podremos aprender $r$ — siempre que nuestra aproximación sea lo suficientemente buena.

Con más detalle, cuando ejecutamos el procedimiento de estimación de fase usando $m$ cúbits de control, lo que obtenemos es un número $y\in\{0,\ldots,2^m-1\}.$ Luego tomamos $y/2^m$ como estimación para $\theta,$ que en este caso es $1/r.$ Para determinar qué es $r$ a partir de esta aproximación, lo natural es calcular el recíproco de nuestra aproximación y redondear al entero más cercano.

$$\left\lfloor \frac{2^m}{y} + \frac{1}{2} \right\rfloor$$

Por ejemplo, supongamos que $r = 6$ y realizamos la estimación de fase sobre $M_a$ con el autovector $\vert\psi_1\rangle$ usando $m = 5$ bits de control. La mejor aproximación de $5$ bits a $1/r = 1/6$ es $5/32,$ y tenemos una buena probabilidad (alrededor del $68\%$ en este caso) de obtener el resultado $y=5$ de la estimación de fase. Tenemos

$$\frac{2^m}{y} = \frac{32}{5} = 6.4,$$

y redondeando al entero más cercano obtenemos $6,$ que es la respuesta correcta.

Por otro lado, si no usamos suficiente precisión, podríamos no obtener la respuesta correcta. Por ejemplo, si tomamos $m = 4$ cúbits de control en la estimación de fase, podríamos obtener la mejor aproximación de $4$ bits a $1/r = 1/6,$ que es $3/16.$ Tomando el recíproco obtenemos

$$\frac{2^m}{y} = \frac{16}{3} = 5.333 \cdots$$

y redondeando al entero más cercano obtenemos la respuesta incorrecta de $5.$

Entonces, ¿cuánta precisión necesitamos para obtener la respuesta correcta? Sabemos que el orden $r$ es un entero, e intuitivamente lo que necesitamos es suficiente precisión para distinguir $1/r$ de posibilidades cercanas, incluyendo $1/(r+1)$ y $1/(r-1).$ El número más cercano a $1/r$ del que debemos preocuparnos es $1/(r+1),$ y la distancia entre estos dos números es

$$\frac{1}{r} - \frac{1}{r+1} = \frac{1}{r(r+1)}.$$

Entonces, si queremos asegurarnos de no confundir $1/r$ con $1/(r+1),$ es suficiente usar suficiente precisión para garantizar que la mejor aproximación $y/2^m$ a $1/r$ sea más cercana a $1/r$ que a $1/(r+1).$ Si usamos suficiente precisión para garantizar que

$$\left\vert \frac{y}{2^m} - \frac{1}{r} \right\vert < \frac{1}{2 r (r+1)},$$

de modo que el error sea menor que la mitad de la distancia entre $1/r$ y $1/(r+1),$ entonces $y/2^m$ estará más cerca de $1/r$ que de cualquier otra posibilidad, incluyendo $1/(r+1)$ y $1/(r-1).$

Podemos verificar esto de la siguiente manera. Supongamos que

$$\frac{y}{2^m} = \frac{1}{r} + \varepsilon$$

para $\varepsilon$ satisfaciendo

$$\vert\varepsilon\vert < \frac{1}{2 r (r+1)}.$$

Cuando tomamos el recíproco obtenemos

$$\frac{2^m}{y} = \frac{1}{\frac{1}{r} + \varepsilon} = \frac{r}{1+\varepsilon r} = r - \frac{\varepsilon r^2}{1+\varepsilon r}.$$

Maximizando en el numerador y minimizando en el denominador, podemos acotar qué tan lejos estamos de $r$ de la siguiente manera.

$$\left\vert \frac{\varepsilon r^2}{1+\varepsilon r} \right\vert \leq \frac{ \frac{r^2}{2 r(r+1)}}{1 - \frac{r}{2r(r+1)}} = \frac{r}{2 r + 1} < \frac{1}{2}$$

Estamos a menos de $1/2$ de distancia de $r,$ por lo que como se esperaba obtendremos $r$ al redondear.

Desafortunadamente, como todavía no sabemos qué es $r,$ no podemos usarlo para indicarnos cuánta precisión necesitamos. Lo que podemos hacer en cambio es usar el hecho de que $r$ debe ser menor que $N$ para asegurarnos de usar suficiente precisión. En particular, si usamos suficiente precisión para garantizar que la mejor aproximación $y/2^m$ a $1/r$ satisfaga

$$\left\vert \frac{y}{2^m} - \frac{1}{r} \right\vert \leq \frac{1}{2N^2},$$

entonces tendremos suficiente precisión para determinar correctamente $r$ cuando tomemos el recíproco. Tomar $m = 2\operatorname{lg}(N)+1$ garantiza que tenemos una alta probabilidad de obtener una estimación con esta precisión usando el método descrito anteriormente. (Tomar $m = 2\operatorname{lg}(N)$ es suficiente si estamos cómodos con un límite inferior del 40% en la probabilidad de éxito.)

Solución general

Como acabamos de ver, si tenemos el autovector $\vert \psi_1 \rangle$ de $M_a,$ podemos aprender $r$ mediante estimación de fase, siempre y cuando usemos suficientes cúbits de control para hacerlo con la precisión suficiente. Desafortunadamente, no es fácil conseguir el autovector $\vert\psi_1\rangle,$ por lo que necesitamos determinar cómo proceder.

Supongamos momentáneamente que procedemos igual que antes, pero con el autovector $\vert\psi_k\rangle$ en lugar de $\vert\psi_1\rangle,$ para cualquier elección de $k\in\{0,\ldots,r-1\}$ que queramos considerar. El resultado que obtenemos del procedimiento de estimación de fase será una aproximación

$$\frac{y}{2^m} \approx \frac{k}{r}.$$

Bajo el supuesto de que no conocemos ni $k$ ni $r,$ esto puede o no permitirnos identificar $r.$ Por ejemplo, si $k = 0$ obtendremos una aproximación $y/2^m$ a $0,$ lo que desafortunadamente no nos dice nada. Sin embargo, este es un caso inusual; para otros valores de $k,$ al menos podremos aprender algo sobre $r.$

Podemos usar un algoritmo conocido como el algoritmo de fracciones continuas para convertir nuestra aproximación $y/2^m$ en fracciones cercanas — incluyendo $k/r$ si la aproximación es lo suficientemente buena. No explicaremos el algoritmo de fracciones continuas aquí. En cambio, aquí está el enunciado de un hecho conocido sobre este algoritmo.

Hecho

Dados un entero $N\geq 2$ y un número real $\alpha\in(0,1),$ existe a lo sumo una elección de enteros $u,v\in\{0,\ldots,N-1\}$ con $v\neq 0$ y $\gcd(u,v)=1$ que satisface $\vert \alpha - u/v\vert < \frac{1}{2N^2}.$ Dados $\alpha$ y $N,$ el algoritmo de fracciones continuas encuentra $u$ y $v,$ o informa que no existen. Este algoritmo puede implementarse como un circuito booleano de tamaño $O((\operatorname{lg}(N))^3).$

Si tenemos una aproximación muy cercana $y/2^m$ a $k/r,$ y ejecutamos el algoritmo de fracciones continuas para $N$ y $\alpha = y/2^m,$ obtendremos $u$ y $v,$ tal como se describen en el hecho. Un análisis del hecho nos permite concluir que

$$\frac{u}{v} = \frac{k}{r}.$$

Nótese en particular que no necesariamente aprendemos $k$ y $r,$ solo aprendemos $k/r$ en su expresión irreducible.

Por ejemplo, y como ya hemos notado, no vamos a aprender nada de $k=0.$ Pero ese es el único valor de $k$ donde eso ocurre. Cuando $k$ es distinto de cero, puede tener factores comunes con $r,$ pero el número $v$ que obtenemos del algoritmo de fracciones continuas debe al menos dividir a $r.$

Aunque no es evidente, es cierto que si tenemos la capacidad de aprender $u$ y $v$ para $u/v = k/r$ con $k\in\{0,\ldots,r-1\}$ elegido uniformemente al azar, entonces es muy probable que podamos recuperar $r$ después de solo unas pocas muestras. En particular, si nuestra estimación para $r$ es el mínimo común múltiplo de todos los valores del denominador $v$ que observamos, tendremos razón con alta probabilidad. Intuitivamente, algunos valores de $k$ no son buenos porque comparten factores comunes con $r,$ y esos factores comunes quedan ocultos cuando aprendemos $u$ y $v.$ Pero las elecciones aleatorias de $k$ no tienden a ocultar factores de $r$ por mucho tiempo, y la probabilidad de no adivinar $r$ correctamente tomando el mínimo común múltiplo de los denominadores observados decrece exponencialmente en el número de muestras.

Queda por abordar la cuestión de cómo conseguimos un autovector $\vert\psi_k\rangle$ de $M_a$ sobre el cual ejecutar el procedimiento de estimación de fase. Resulta que ¡en realidad no necesitamos crearlos!

Lo que haremos en su lugar es ejecutar el procedimiento de estimación de fase sobre el estado $\vert 1\rangle,$ con lo que nos referimos a la codificación binaria de $n$ bits del número $1,$ en lugar de un autovector $\vert\psi\rangle$ de $M_a.$ Hasta ahora solo hemos hablado de ejecutar el procedimiento de estimación de fase sobre un autovector particular, pero nada nos impide ejecutar el procedimiento sobre un estado de entrada que no sea un autovector de $M_a,$ y eso es lo que hacemos aquí con el estado $\vert 1\rangle.$ (Esto no es un autovector de $M_a$ a menos que $a=1,$ que no es una elección que nos interesará.)

La justificación para elegir el estado $\vert 1\rangle$ en lugar de un autovector de $M_a$ es que la siguiente ecuación es verdadera.

$$\vert 1\rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \vert \psi_k\rangle$$

Una forma de verificar esta ecuación es comparar los productos internos de ambos lados con cada estado de la base estándar, usando fórmulas mencionadas anteriormente en la lección para ayudar a evaluar los resultados del lado derecho. Como consecuencia, obtendremos exactamente los mismos resultados de medición que si hubiéramos elegido $k\in\{0,\ldots,r-1\}$ uniformemente al azar y usado $\vert\psi_k\rangle$ como autovector.

Con mayor detalle, imaginemos que ejecutamos el procedimiento de estimación de fase con el estado $\vert 1\rangle$ en lugar de uno de los autovectores $\vert\psi_k\rangle.$ Después de realizar la transformada cuántica de Fourier inversa, esto nos deja con el estado

$$\frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \vert \psi_k\rangle \vert \gamma_k\rangle,$$

donde

$$\vert\gamma_k\rangle = \frac{1}{2^m} \sum_{y=0}^{2^m - 1} \sum_{x=0}^{2^m-1} e^{2\pi i x (k/r - y/2^m)} \vert y\rangle.$$

El vector $\vert\gamma_k\rangle$ representa el estado de los $m$ cúbits superiores después de que se ha aplicado la inversa de la transformada cuántica de Fourier sobre ellos.

Entonces, en virtud de que $\{\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle\}$ es un conjunto ortonormal, encontramos que una medición de los $m$ cúbits superiores produce una aproximación $y/2^m$ al valor $k/r$ donde $k\in\{0,\ldots,r-1\}$ se elige uniformemente al azar. Como ya hemos discutido, esto nos permite aprender $r$ con un alto grado de confianza después de varias ejecuciones independientes, que era nuestro objetivo.

Costo total

El costo de implementar cada unitaria controlada $M_a^k$ es $O(n^2).$ Hay $m$ operaciones unitarias controladas, y tenemos $m = O(n),$ por lo que el costo total de las operaciones unitarias controladas es $O(n^3).$ Además, tenemos $m$ compuertas Hadamard (que contribuyen $O(n)$ al costo), y la transformada cuántica de Fourier inversa contribuye $O(n^2)$ al costo. Por lo tanto, el costo de las operaciones unitarias controladas domina el costo de todo el procedimiento — que es entonces $O(n^3).$

Además del circuito cuántico en sí, hay algunos cálculos clásicos que deben realizarse en el camino. Esto incluye calcular las potencias $a^k$ en $\mathbb{Z}_N$ para $k = 2, 4, 8, \ldots, 2^{m-1},$ que se necesitan para crear las compuertas unitarias controladas, así como el algoritmo de fracciones continuas que convierte aproximaciones de $\theta$ en fracciones. Estos cálculos pueden realizarse mediante circuitos booleanos con un costo total de $O(n^3).$

Como es habitual, todas estas cotas pueden mejorarse usando algoritmos asintóticamente rápidos; estas cotas asumen que estamos usando algoritmos estándar para las operaciones aritméticas básicas.

Factorización mediante búsqueda de orden

Lo último que necesitamos discutir es cómo resolver el problema de búsqueda de orden nos ayuda a factorizar. Esta parte es completamente clásica — no tiene nada que ver específicamente con la computación cuántica.

Aquí está la idea básica. Queremos factorizar el número $N,$ y podemos hacer esto recursivamente. Específicamente, podemos enfocarnos en la tarea de dividir $N,$ que significa encontrar dos enteros cualesquiera $b,c\geq 2$ para los cuales $N = bc.$ Esto no es posible si $N$ es un número primo, pero podemos probar eficientemente si $N$ es primo usando un algoritmo de prueba de primalidad primero, y si $N$ no es primo intentaremos dividirlo. Una vez que dividimos $N,$ podemos simplemente recurrir sobre $b$ y $c$ hasta que todos nuestros factores sean primos y obtengamos la factorización prima de $N.$

Dividir enteros pares es fácil: simplemente damos como resultado $2$ y $N/2.$

También es fácil dividir potencias perfectas, es decir, números de la forma $N = s^j$ para enteros $s,j\geq 2,$ simplemente aproximando las raíces $N^{1/2},$ $N^{1/3},$ $N^{1/4},$ y así sucesivamente, y comprobando los enteros cercanos como candidatos para $s.$ No necesitamos ir más allá de $\log(N)$ pasos en esta secuencia, porque en ese punto la raíz cae por debajo de $2$ y no revelará candidatos adicionales.

Es bueno que podamos hacer ambas cosas porque la búsqueda de orden no nos ayudará a factorizar números pares ni potencias primas, donde el número $s$ resulta ser primo. Si $N$ es impar y no es una potencia prima, sin embargo, la búsqueda de orden nos permite dividir $N.$

Algoritmo probabilístico para dividir un entero compuesto impar N que no es una potencia prima

1. Elige aleatoriamente $a\in\{2,\ldots,N-1\}.$

2. Calcula $d=\gcd(a,N).$

3. Si $d > 1$ entonces da como resultado $b = d$ y $c = N/d$ y detente. De lo contrario continúa con el siguiente paso sabiendo que $a\in\mathbb{Z}_N^{\ast}.$

4. Sea $r$ el orden de $a$ módulo $N.$ (Aquí es donde necesitamos la búsqueda de orden.)

5. Si $r$ es par:

   5.1 Calcula $x = a^{r/2} - 1$ módulo $N$
   5.2 Calcula $d = \gcd(x,N).$
   5.3 Si $d>1$ entonces da como resultado $b=d$ y $c = N/d$ y detente.

6. Si se llega a este punto, el algoritmo no ha podido encontrar un factor de $N.$

Una ejecución de este algoritmo puede no encontrar un factor de $N.$ Específicamente, esto ocurre en dos situaciones:

• El orden de $a$ módulo $N$ es impar.
• El orden de $a$ módulo $N$ es par y $\gcd\bigl(a^{r/2} - 1, N\bigr) = 1.$

Usando teoría básica de números se puede probar que, para una elección aleatoria de $a,$ con probabilidad al menos $1/2$ ninguno de estos eventos ocurre. De hecho, la probabilidad de que ocurra alguno de estos eventos es a lo sumo $2^{-(m-1)}$ siendo $m$ el número de factores primos distintos de $N,$ razón por la cual se necesita el supuesto de que $N$ no es una potencia prima. (El supuesto de que $N$ es impar también es necesario para que este hecho sea verdadero.)

Esto significa que cada ejecución tiene al menos un 50% de probabilidad de dividir $N.$ Por lo tanto, si ejecutamos el algoritmo $t$ veces, eligiendo $a$ aleatoriamente cada vez, tendremos éxito en dividir $N$ con probabilidad al menos $1 - 2^{-t}.$

La idea básica detrás del algoritmo es la siguiente. Si tenemos una elección de $a$ para la cual el orden $r$ de $a$ módulo $N$ es par, entonces $r/2$ es un entero y podemos considerar los números

$$a^{r/2} - 1\; (\textrm{mod}\; N) \quad \text{y} \quad a^{r/2} + 1\; (\textrm{mod}\; N).$$

Usando la fórmula $Z^2 - 1 = (Z+1)(Z-1),$ concluimos que

$$\bigl(a^{r/2} - 1\bigr) \bigl(a^{r/2} + 1\bigr) = a^r - 1.$$

Ahora, sabemos que $a^r \; (\textrm{mod}\; N) = 1$ por la definición del orden — lo que es otra forma de decir que $N$ divide exactamente a $a^r - 1.$ Eso significa que $N$ divide exactamente al producto

$$\bigl(a^{r/2} - 1\bigr) \bigl(a^{r/2} + 1\bigr).$$

Para que esto sea verdad, todos los factores primos de $N$ deben ser también factores primos de $a^{r/2} - 1$ o de $a^{r/2} + 1$ (o de ambos) — y para una selección aleatoria de $a$ resulta ser poco probable que todos los factores primos de $N$ dividan a uno de los términos sin que ninguno divida al otro. De lo contrario, siempre que algunos de los factores primos de $N$ dividan al primer término y algunos dividan al segundo término, podremos encontrar un factor no trivial de $N$ calculando el MCD con el primer término.