Algoritmo de Shor

Estimación de uso: Tres segundos en un procesador Eagle r3 (NOTA: Esto es sólo una estimación. Su tiempo de ejecución puede variar)

Resultados del aprendizaje

Una vez completado este tutorial, los usuarios deberían comprender:

• Los fundamentos matemáticos del algoritmo de Shor para la factorización de números enteros
• Cómo ejecutar una instancia de ejemplo de este algoritmo en un equipo

Requisitos previos

Recomendamos a los usuarios que se familiaricen con los siguientes temas antes de seguir este tutorial:

• Fundamentos de los algoritmos cuánticos.
• Estimación de fase y factorización. En este tutorial tratamos parte de este tema.

En segundo plano

El algoritmo de Shor, desarrollado por Peter Shor en 1994, es un algoritmo cuántico revolucionario que permite factorizar números enteros en tiempo polinomial. Su importancia radica en su capacidad para factorizar números enteros grandes de forma exponencialmente más rápida que cualquier algoritmo clásico conocido, lo que pone en peligro la seguridad de sistemas criptográficos de uso generalizado como el RSA, que se basan en la dificultad de factorizar números grandes. Si se resolviera este problema de manera eficiente en un ordenador cuántico lo suficientemente potente, el algoritmo de Shor podría revolucionar campos como la criptografía, la ciberseguridad y las matemáticas computacionales, lo que pondría de relieve el poder transformador de la computación cuántica.

Este tutorial se centra en la demostración del algoritmo de Shor mediante la factorización de 15 en un ordenador cuántico.

En primer lugar, definimos el problema de búsqueda de orden y construimos los circuitos correspondientes a partir del protocolo cuántico de estimación de fase. A continuación, ejecutamos los circuitos de búsqueda de órdenes en hardware real utilizando circuitos de menor profundidad que podamos transpilar. La última sección completa el algoritmo de Shor conectando el problema de búsqueda de órdenes con la factorización entera.

Terminamos el tutorial con una discusión sobre otras demostraciones del algoritmo de Shor en hardware real, centrándonos tanto en implementaciones genéricas como en aquellas adaptadas a la factorización de enteros específicos como 15 y 21.

Nota: Este tutorial se centra más en la implementación y demostración de los circuitos relativos al algoritmo de Shor. Para profundizar en el material, consulte el curso Fundamentos de los algoritmos cuánticos, del Dr. John Watrous, y los artículos de la sección Referencias.

Requisitos

Antes de empezar este tutorial, asegúrate de que tienes instalado lo siguiente:

• Qiskit SDK v2.0 o posterior, con soporte de visualización
• Qiskit Runtime v0.40 o posterior (pip install qiskit-ibm-runtime)

Configuración

import numpy as np
import pandas as pd
from fractions import Fraction
from math import floor, gcd, log

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit.transpiler import CouplingMap, generate_preset_pass_manager
from qiskit.visualization import plot_histogram

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Paso 1: Asignar entradas clásicas a un problema cuántico

El algoritmo de Shor para la factorización de enteros utiliza un problema intermedio conocido como problema de búsqueda de órdenes. En esta sección, demostramos cómo resolver el problema de búsqueda de orden utilizando la estimación cuántica de fase.

Problema de estimación de fase

En el problema de estimación de fase, se nos da un estado cuántico $\ket{\psi}$ de $n$ qubits, junto con un circuito cuántico unitario que actúa sobre $n$ qubits. Nos han prometido que $\ket{\psi}$ es un vector propio de la matriz unitaria $U$ que describe la acción del circuito, y nuestro objetivo es calcular o aproximar el valor propio $\lambda = e^{2 \pi i \theta}$ al que corresponde $\ket{\psi}$. En otras palabras, el circuito debe dar como salida una aproximación al número $\theta \in [0, 1)$ que satisfaga $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ El objetivo del circuito de estimación de fase es aproximar $\theta$ en $m$ bits. Matemáticamente hablando, nos gustaría encontrar $y$ tal que $\theta \approx y / 2^m$, donde $y \in {0, 1, 2, \dots, 2^{m-1}}$. La siguiente imagen muestra el circuito cuántico que estima $y$ en $m$ bits realizando una medida en $m$ qubits.

En el circuito anterior, los qubits superiores $m$ se inician en el estado $\ket{0^m}$, y los qubits inferiores $n$ se inician en $\ket{\psi}$, que se promete que es un eigenvector de $U$. El primer ingrediente en el circuito de estimación de fase son las operaciones unitarias controladas que son responsables de realizar un retroceso de fase a su qubit de control correspondiente. Estos unitarios controlados se exponencian en función de la posición del qubit de control, desde el bit menos significativo hasta el bit más significativo. Dado que $\ket{\psi}$ es un eigenvector de $U$, el estado de los qubits inferiores $n$ no se ve afectado por esta operación, pero la información de fase del eigenvalor se propaga a los qubits superiores $m$.

Resulta que después de la operación de retroceso de fase mediante unidades controladas, todos los estados posibles de los qubits superiores $m$ son ortonormales entre sí para cada vector propio $\ket{\psi}$ del unitario $U$. Por lo tanto, estos estados son perfectamente distinguibles, y podemos rotar la base que forman de nuevo a la base computacional para hacer una medición. Un análisis matemático muestra que esta matriz de rotación corresponde a la transformada cuántica de Fourier (QFT) inversa en el espacio de Hilbert $2^m$ -dimensional. La intuición detrás de esto es que la estructura periódica de los operadores de exponenciación modular está codificada en el estado cuántico, y la QFT convierte esta periodicidad en picos medibles en el dominio de la frecuencia.

Para una comprensión más profunda de por qué se emplea el circuito QFT en el algoritmo de Shor, remitimos al lector al curso Fundamentos de algoritmos cuánticos.

Ahora estamos preparados para utilizar el circuito de estimación de fase para la búsqueda de órdenes.

Problema con el pedido

Para definir el problema de búsqueda de órdenes, empezamos con algunos conceptos de teoría de números. En primer lugar, para cualquier número entero positivo dado $N$, defina el conjunto $\mathbb{Z}_N$ como $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Todas las operaciones aritméticas en $\mathbb{Z}_N$ se realizan modulo $N$. En particular, todos los elementos $a \in \mathbb{Z}_n$ que son coprimos con $N$ son especiales y constituyen $\mathbb{Z}^*_N$ como $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Para un elemento $a \in \mathbb{Z}^*_N$, el menor número entero positivo $r$ tal que $a^r \equiv 1 \; (\mathrm{mod} \; N)$ se define como el orden de $a$ módulo $N$. Como veremos más adelante, encontrar el orden de un $a \in \mathbb{Z}^*_N$ nos permitirá factorizar $N$.

Para construir el circuito de búsqueda de orden a partir del circuito de estimación de fase, necesitamos dos consideraciones. En primer lugar, necesitamos definir el unitario $U$ que nos permitirá encontrar el orden $r$, y en segundo lugar, necesitamos definir un vector propio $\ket{\psi}$ de $U$ para preparar el estado inicial del circuito de estimación de fase.

Para conectar el problema de búsqueda de orden con la estimación de fase, consideramos 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$. En particular, definimos este operador de multiplicación $M_a$ tal que $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ para cada $x \in \mathbb{Z}_N$. Nótese que está implícito que estamos tomando el producto modulo $N$ dentro del ket en el lado derecho de la ecuación. Un análisis matemático muestra que $M_a$ es un operador unitario. Además, resulta que $M_a$ tiene pares de vectores y valores propios que nos permiten conectar el orden $r$ de $a$ con el problema de estimación de fase. En concreto, para cualquier elección de $j \in \{0, \dots, r-1\}$, tenemos que $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ es un vector propio de $M_a$ cuyo valor propio correspondiente es $\omega^{j}_{r}$, donde $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$

Por observación, vemos que un par eigenvector/eigenvalor conveniente es el estado $\ket{\psi_1}$ con $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Por lo tanto, si pudiéramos encontrar el vector propio $\ket{\psi_1}$, podríamos estimar la fase $\theta=1/r$ con nuestro circuito cuántico y, por lo tanto, obtener una estimación del orden $r$. Sin embargo, no es fácil hacerlo, y tenemos que considerar una alternativa.

Consideremos cuál sería el resultado del circuito si preparamos el estado computacional $\ket{1}$ como estado inicial. No se trata de un estado propio de $M_a$, sino de la superposición uniforme de los estados propios que acabamos de describir. En otras palabras, se cumple la siguiente relación. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ La implicación de la ecuación anterior es que si fijamos el estado inicial en $\ket{1}$, obtendremos precisamente el mismo resultado de medición que si hubiéramos elegido $k \in \{ 0, \dots, r-1\}$ uniformemente al azar y utilizado $\ket{\psi_k}$ como un vector propio en el circuito de estimación de fase. En otras palabras, una medición de los $m$ qubits superiores produce una aproximación $y / 2^m$ al valor $k / r$ donde $k \in \{ 0, \dots, r-1\}$ se elige uniformemente al azar. Esto nos permite aprender $r$ con un alto grado de confianza tras varias ejecuciones independientes, que era nuestro objetivo.

Operadores de exponenciación modular

Hasta ahora, hemos vinculado el problema de estimación de fase al problema de búsqueda de orden definiendo $U = M_a$ y $\ket{\psi} = \ket{1}$ en nuestro circuito cuántico. Por lo tanto, el último ingrediente que queda es encontrar una manera eficiente de definir exponenciales modulares de $M_a$ como $M_a^k$ para $k = 1, 2, 4, \dots, 2^{m-1}$. Para realizar este cálculo, encontramos que 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 \; \mathrm{mod} \; N$ y luego usando el circuito para $M_b$. Dado que sólo necesitamos las potencias que son potencias de 2, podemos hacer esto de forma clásica y eficiente utilizando el cuadrado iterativo.

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

Ejemplo específico con $N = 15$ y $a=2$

Podemos detenernos aquí para discutir un ejemplo específico y construir el circuito de búsqueda de orden para $N=15$. Obsérvese que los posibles $a \in \mathbb{Z}_N^*$ no triviales para $N=15$ son $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Para este ejemplo, elegimos $a=2$. Construiremos el operador $M_2$ y los operadores de exponenciación modular $M_2^k$.

La acción de $M_2$ sobre los estados de la base de cálculo es la siguiente. $M_2 \ket{0} = \ket{0} \quad M_2 \ket{5} = \ket{10} \quad M_2 \ket{10} = \ket{5}$ $M_2 \ket{1} = \ket{2} \quad M_2 \ket{6} = \ket{12} \quad M_2 \ket{11} = \ket{7}$ $M_2 \ket{2} = \ket{4} \quad M_2 \ket{7} = \ket{14} \quad M_2 \ket{12} = \ket{9}$ $M_2 \ket{3} = \ket{6} \quad M_2 \ket{8} = \ket{1} \quad M_2 \ket{13} = \ket{11}$ $M_2 \ket{4} = \ket{8} \quad M_2 \ket{9} = \ket{3} \quad M_2 \ket{14} = \ket{13}$ Por observación, podemos ver que los estados base se barajan, por lo que tenemos una matriz de permutación. Podemos construir esta operación en cuatro qubits con puertas de intercambio. A continuación, construimos las operaciones $M_2$ y $M_2$ controlada.

def M2mod15():
    """
    M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M2 operator
M2 = M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

def controlled_M2mod15():
    """
    Controlled M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Output:

Las puertas que actúan sobre más de dos qubits se descompondrán en puertas de dos qubits.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

Ahora necesitamos construir los operadores de exponenciación modular. Para obtener suficiente precisión en la estimación de fase, utilizaremos ocho qubits para la medición de la estimación. Por lo tanto, necesitamos construir $M_b$ con $b = a^{2^k} \; (\mathrm{mod} \; N)$ para cada $k = 0, 1, \dots, 7$.

def a2kmodN(a, k, N):
    """Compute a^{2^k} (mod N) by repeated squaring"""
    for _ in range(k):
        a = int(np.mod(a**2, N))
    return a

k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]

print(b_list)

Output:

[2, 4, 1, 1, 1, 1, 1, 1]

Como podemos ver en la lista de valores $b$, además de $M_2$ que construimos previamente, también necesitamos construir $M_4$ y $M_1$. Nótese que $M_1$ actúa trivialmente sobre los estados de la base computacional, por lo que es simplemente el operador identidad.

$M_4$ actúa sobre los estados base computacionales de la siguiente manera. $M_4 \ket{0} = \ket{0} \quad M_4 \ket{5} = \ket{5} \quad M_4 \ket{10} = \ket{10}$ $M_4 \ket{1} = \ket{4} \quad M_4 \ket{6} = \ket{9} \quad M_4 \ket{11} = \ket{14}$ $M_4 \ket{2} = \ket{8} \quad M_4 \ket{7} = \ket{13} \quad M_4 \ket{12} = \ket{3}$ $M_4 \ket{3} = \ket{12} \quad M_4 \ket{8} = \ket{2} \quad M_4 \ket{13} = \ket{7}$ $M_4 \ket{4} = \ket{1} \quad M_4 \ket{9} = \ket{6} \quad M_4 \ket{14} = \ket{11}$

Por lo tanto, esta permutación se puede construir con la siguiente operación de intercambio.

def M4mod15():
    """
    M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M4 operator
M4 = M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

def controlled_M4mod15():
    """
    Controlled M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Output:

Las puertas que actúan sobre más de dos qubits se descompondrán en puertas de dos qubits.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

Hemos visto que los operadores $M_b$ para un $b \in \mathbb{Z}^*_N$ dado son operaciones de permutación. Debido al tamaño relativamente pequeño del problema de permutación que tenemos aquí, ya que $N=15$ requiere sólo cuatro qubits, pudimos sintetizar estas operaciones directamente con SWAP gates por inspección. En general, éste podría no ser un enfoque escalable. En su lugar, podríamos necesitar construir la matriz de permutación explícitamente, y utilizar la clase UnitaryGate de Qiskit y los métodos de transpilación para sintetizar esta matriz de permutación. Sin embargo, esto puede dar lugar a circuitos mucho más profundos. A continuación se proporciona un ejemplo.

def mod_mult_gate(b, N):
    """
    Modular multiplication gate from permutation matrix.
    """
    if gcd(b, N) > 1:
        print(f"Error: gcd({b},{N}) > 1")
    else:
        n = floor(log(N - 1, 2)) + 1
        U = np.full((2**n, 2**n), 0)
        for x in range(N):
            U[b * x % N][x] = 1
        for x in range(N, 2**n):
            U[x][x] = 1
        G = UnitaryGate(U)
        G.name = f"M_{b}"
        return G

# Let's build M2 using the permutation matrix definition
M2_other = mod_mult_gate(2, 15)

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2_other, inplace=True)
circ = circ.decompose()

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

print(f"qubits: {circ.num_qubits}")
print(
    f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.decompose().draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

Output:

qubits: 4
2q-depth: 94
2q-size: 96
Operator counts: OrderedDict({'cx': 45, 'swap': 32, 'u': 24, 'u1': 7, 'u3': 4, 'unitary': 3, 'circuit-335': 1, 'circuit-338': 1, 'circuit-341': 1, 'circuit-344': 1, 'circuit-347': 1, 'circuit-350': 1, 'circuit-353': 1, 'circuit-356': 1, 'circuit-359': 1, 'circuit-362': 1, 'circuit-365': 1, 'circuit-368': 1, 'circuit-371': 1, 'circuit-374': 1, 'circuit-377': 1, 'circuit-380': 1})

Comparemos estos recuentos con la profundidad del circuito compilado de nuestra implementación manual de la puerta $M_2$.

# Get the M2 operator from our manual construction
M2 = M2mod15()

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ = circ.decompose(reps=3)

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

print(f"qubits: {circ.num_qubits}")
print(
    f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

Output:

qubits: 4
2q-depth: 9
2q-size: 9
Operator counts: OrderedDict({'cx': 9})

Como podemos ver, el enfoque de la matriz de permutación resultó en un circuito significativamente profundo incluso para una única puerta $M_2$ en comparación con nuestra implementación manual del mismo. Por lo tanto, continuaremos con nuestra aplicación anterior de las operaciones $M_b$.

Ahora, estamos listos para construir el circuito de búsqueda de orden completo utilizando nuestros operadores de exponenciación modular controlada previamente definidos. En el siguiente código, también importamos el circuito QFT de la librería Qiskit Circuit, que utiliza puertas Hadamard en cada qubit, una serie de puertas controlled-U1 (o Z, dependiendo de la fase), y una capa de puertas swap.

# Order finding problem for N = 15 with a = 2
N = 15
a = 2

# Number of qubits
num_target = floor(log(N - 1, 2)) + 1  # for modular exponentiation operators
num_control = 2 * num_target  # for enough precision of estimation

# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]

# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)

# Initialize the target register to the state |1>
circuit.x(num_control)

# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
    circuit.h(k)
    b = b_list[k]
    if b == 2:
        circuit.compose(
            M2mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    elif b == 4:
        circuit.compose(
            M4mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    else:
        continue  # M1 is the identity operator

# Apply the inverse QFT to the control register
circuit.compose(QFT(num_control, inverse=True), qubits=control, inplace=True)

# Measure the control register
circuit.measure(control, output)

circuit.draw("mpl", fold=-1)

Output:

Nótese que omitimos las operaciones de exponenciación modular controlada de los qubits de control restantes porque $M_1$ es el operador identidad.

Tenga en cuenta que más adelante en este tutorial, ejecutaremos este circuito en el backend ibm_marrakesh . Para ello, transpilamos el circuito según este backend específico e informamos de la profundidad del circuito y el recuento de puertas.

service = QiskitRuntimeService()
backend = service.backend("ibm_marrakesh")
pm = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuit = pm.run(circuit)

print(
    f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}"
)
print(
    f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}"
)
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

Output:

2q-depth: 187
2q-size: 260
Operator counts: OrderedDict({'sx': 521, 'rz': 354, 'cz': 260, 'measure': 8, 'x': 4})

Paso 3: Ejecutar utilizando Qiskit primitives

En primer lugar, analizamos lo que obtendríamos teóricamente si ejecutáramos este circuito en un simulador ideal. A continuación, tenemos una serie de resultados de simulación del circuito anterior utilizando 1024 tomas. Como vemos, obtenemos una distribución aproximadamente uniforme en cuatro cadenas de bits sobre los qubits de control.

# Obtained from the simulator
counts = {"00000000": 264, "01000000": 268, "10000000": 249, "11000000": 243}

plot_histogram(counts)

Output:

Midiendo los qubits de control, obtenemos una estimación de fase de ocho bits del operador $M_a$. Podemos convertir esta representación binaria a decimal para hallar la fase medida. Como podemos ver en el histograma anterior, se midieron cuatro cadenas de bits diferentes, y cada una de ellas corresponde a un valor de fase, como se indica a continuación.

# Rows to be displayed in table
rows = []
# Corresponding phase of each bitstring
measured_phases = []

for output in counts:
    decimal = int(output, 2)  # Convert bitstring to decimal
    phase = decimal / (2**num_control)  # Find corresponding eigenvalue
    measured_phases.append(phase)
    # Add these values to the rows in our table:
    rows.append(
        [
            f"{output}(bin) = {decimal:>3}(dec)",
            f"{decimal}/{2 ** num_control} = {phase:.2f}",
        ]
    )

# Print the rows in a table
headers = ["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Output:

            Register Output           Phase
0  00000000(bin) =   0(dec)    0/256 = 0.00
1  01000000(bin) =  64(dec)   64/256 = 0.25
2  10000000(bin) = 128(dec)  128/256 = 0.50
3  11000000(bin) = 192(dec)  192/256 = 0.75

Recordemos que la fase medida cualquiera corresponde a $\theta = k / r$ donde $k$ se muestrea uniformemente al azar a partir de $\{0, 1, \dots, r-1 \}$. Por lo tanto, podemos utilizar el algoritmo de fracciones continuas para intentar encontrar $k$ y el orden $r$. Python tiene esta funcionalidad incorporada. Podemos utilizar el módulo fractions para convertir un flotador en un objeto Fraction , por ejemplo:

Fraction(0.666)

Output:

Fraction(5998794703657501, 9007199254740992)

Debido a que esto da fracciones que devuelven el resultado exactamente (en este caso, 0.6660000...), esto puede dar resultados retorcidos como el de arriba. Podemos utilizar el método .limit_denominator() para obtener la fracción que más se parezca a nuestro flotante, con un denominador inferior a un determinado valor:

# Get fraction that most closely resembles 0.666
# with denominator < 15
Fraction(0.666).limit_denominator(15)

Output:

Fraction(2, 3)

Esto es mucho más bonito. El orden (r) debe ser menor que N, por lo que fijaremos el denominador máximo en 15:

# Rows to be displayed in a table
rows = []

for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append(
        [phase, f"{frac.numerator}/{frac.denominator}", frac.denominator]
    )

# Print the rows in a table
headers = ["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Output:

   Phase Fraction  Guess for r
0   0.00      0/1            1
1   0.25      1/4            4
2   0.50      1/2            2
3   0.75      3/4            4

Podemos ver que dos de los valores propios medidos nos dieron el resultado correcto: $r=4$, y podemos ver que el algoritmo de Shor para encontrar el orden tiene una posibilidad de fallar. Estos malos resultados se deben a que $k = 0$, o porque $k$ y $r$ no son coprimos - y en lugar de $r$, se nos da un factor de $r$. La solución más fácil para esto es simplemente repetir el experimento hasta que obtengamos un resultado satisfactorio para $r$.

Hasta ahora, hemos implementado el problema de búsqueda de orden para $N=15$ con $a=2$ utilizando el circuito de estimación de fase en un simulador. El último paso del algoritmo de Shor consistirá en relacionar el problema de búsqueda de órdenes con el problema de factorización de enteros. Esta última parte del algoritmo es puramente clásica y puede resolverse en un ordenador clásico una vez obtenidas las medidas de fase en un ordenador cuántico. Por lo tanto, aplazamos la última parte del algoritmo hasta después de demostrar cómo podemos ejecutar el circuito de búsqueda de órdenes en hardware real.

Ejecución de hardware

Ahora podemos ejecutar el circuito de búsqueda de órdenes que transpilamos anteriormente para ibm_marrakesh. Aquí nos centramos en el desacoplamiento dinámico (DD) para la supresión de errores, y en el giro de compuertas para la mitigación de errores. La DD consiste en aplicar secuencias de pulsos de control sincronizados con precisión a un dispositivo cuántico, lo que elimina las interacciones ambientales no deseadas y la decoherencia. El giro de puertas, por su parte, aleatoriza puertas cuánticas específicas para transformar los errores coherentes en errores Pauli, que se acumulan linealmente en lugar de cuadráticamente. Ambas técnicas se combinan a menudo para mejorar la coherencia y fidelidad de los cálculos cuánticos.

# Sampler primitive to obtain the probability distribution
sampler = Sampler(backend)

# Turn on dynamical decoupling with sequence XpXm
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XpXm"
# Enable gate twirling
sampler.options.twirling.enable_gates = True

# Assign tags before executing
sampler.options.environment.job_tags = ["TUT_SA"]

pub = transpiled_circuit
job = sampler.run([pub], shots=1024)

result = job.result()[0]
counts = result.data["out"].get_counts()

plot_histogram(counts, figsize=(35, 5))

Output:

Como podemos ver, obtuvimos las mismas cadenas de bits con los recuentos más altos. Como el hardware cuántico tiene ruido, hay algunas fugas a otras cadenas de bits, que podemos filtrar estadísticamente.

# Dictionary of bitstrings and their counts to keep
counts_keep = {}
# Threshold to filter
threshold = np.max(list(counts.values())) / 2

for key, value in counts.items():
    if value > threshold:
        counts_keep[key] = value

print(counts_keep)

Output:

{'00000000': 58, '01000000': 41, '11000000': 42, '10000000': 40}

Paso 4: Procesamiento posterior y devolución del resultado en el formato clásico deseado

Factorización de números enteros

Hasta ahora, hemos discutido cómo podemos implementar el problema de búsqueda de orden utilizando un circuito de estimación de fase. Ahora, conectamos el problema de búsqueda de orden a la factorización entera, lo que completa el algoritmo de Shor. Tenga en cuenta que esta parte del algoritmo es clásica.

Ahora lo demostraremos con nuestro ejemplo de $N = 15$ y $a = 2$. Recordemos que la fase que medimos es $k / r$, donde $a^r \; (\textrm{mod} \; N) = 1$ y $k$ es un número entero aleatorio entre $0$ y $r - 1$. A partir de esta ecuación, tenemos $(a^r - 1) \; (\textrm{mod} \; N) = 0,$, lo que significa que $N$ debe dividir a $a^r-1$. Si $r$ también es par, entonces podemos escribir $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Si $r$ no es par, no podemos ir más allá y debemos intentarlo de nuevo con un valor diferente para $a$; de lo contrario, hay una alta probabilidad de que el máximo común divisor de $N$ y $a^{r/2}-1$, o $a^{r/2}+1$ sea un factor propio de $N$.

Dado que algunas ejecuciones del algoritmo fallarán estadísticamente, repetiremos este algoritmo hasta que se encuentre al menos un factor de $N$.

La celda siguiente repite el algoritmo hasta que se encuentra al menos un factor de $N=15$. Utilizaremos los resultados de la ejecución de hardware anterior para adivinar la fase y el factor correspondiente en cada iteración.

a = 2
N = 15

FACTOR_FOUND = False
num_attempt = 0

while not FACTOR_FOUND:
    print(f"\nATTEMPT {num_attempt}:")
    # Here, we get the bitstring by iterating over outcomes
    # of a previous hardware run with multiple shots.
    # Instead, we can also perform a single-shot measurement
    # here in the loop.
    bitstring = list(counts_keep.keys())[num_attempt]
    num_attempt += 1
    # Find the phase from measurement
    decimal = int(bitstring, 2)
    phase = decimal / (2**num_control)  # phase = k / r
    print(f"Phase: theta = {phase}")

    # Guess the order from phase
    frac = Fraction(phase).limit_denominator(N)
    r = frac.denominator  # order = r
    print(f"Order of {a} modulo {N} estimated as: r = {r}")

    if phase != 0:
        # Guesses for factors are gcd(a^{r / 2} ± 1, 15)
        if r % 2 == 0:
            x = pow(a, r // 2, N) - 1
            d = gcd(x, N)
            if d > 1:
                FACTOR_FOUND = True
                print(f"*** Non-trivial factor found: {x} ***")

Output:

ATTEMPT 0:
Phase: theta = 0.0
Order of 2 modulo 15 estimated as: r = 1

ATTEMPT 1:
Phase: theta = 0.25
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

Debate

Trabajos relacionados

En esta sección, analizamos otros trabajos que han demostrado el algoritmo de Shor en hardware real.

El trabajo seminal [3] de IBM® demostró por primera vez el algoritmo de Shor, factorizando el número 15 en sus factores primos 3 y 5 utilizando un ordenador cuántico de resonancia magnética nuclear (RMN) de siete qubits. Otro experimento [4] factorizó 15 utilizando qubits fotónicos. Empleando un único qubit reciclado varias veces y codificando el registro de trabajo en estados de mayor dimensión, los investigadores redujeron el número necesario de qubits a un tercio del del protocolo estándar, utilizando un algoritmo compilado de dos fotones. Un artículo significativo en la demostración del algoritmo de Shor es [5], que utiliza la técnica de estimación de fase iterativa de Kitaev [8] para reducir el requisito de qubits del algoritmo. Los autores utilizaron siete qubits de control y cuatro qubits de caché, junto con la implementación de multiplicadores modulares. Esta implementación, sin embargo, requiere mediciones en mitad del circuito con operaciones feed-forward y reciclado de qubits con operaciones de reinicio. Esta demostración se realizó en un ordenador cuántico de trampa de iones.

Un trabajo más reciente [6] se centró en la factorización de 15, 21 y 35 en el hardware IBM Quantum®. Al igual que en trabajos anteriores, los investigadores utilizaron una versión compilada del algoritmo que empleaba una transformada cuántica de Fourier semiclásica como la propuesta por Kitaev para minimizar el número de qubits físicos y puertas. Un trabajo más reciente [7] también realizó una demostración de prueba de concepto para factorizar el entero 21. Esta demostración también implicó el uso de una versión compilada de la rutina de estimación de fase cuántica, y se basó en la demostración anterior de [4]. Los autores fueron más allá de este trabajo utilizando una configuración de puertas Toffoli aproximadas con desplazamientos de fase residuales. El algoritmo se implementó en procesadores cuánticos IBM utilizando sólo cinco qubits, y se verificó con éxito la presencia de entrelazamiento entre los qubits de control y registro.

Escalado del algoritmo

El cifrado RSA suele requerir claves del orden de 2048 a 4096 bits. Intentar factorizar un número de 2048 bits con el algoritmo de Shor dará como resultado un circuito cuántico con millones de qubits, incluida la sobrecarga de corrección de errores y una profundidad de circuito del orden de mil millones, que está más allá de los límites de ejecución del hardware cuántico actual. Por lo tanto, el algoritmo de Shor necesitará métodos de construcción de circuitos optimizados o una corrección de errores cuántica robusta para ser viable en la práctica a la hora de romper sistemas criptográficos modernos. Le remitimos a [9] para una discusión más detallada sobre la estimación de recursos para el algoritmo de Shor.

Reto

¡Enhorabuena por terminar el tutorial! Ahora es un buen momento para poner a prueba tus conocimientos. ¿Podrías intentar construir el circuito para factorizar 21? Puede seleccionar un $a$ de su elección. Tendrás que decidir la precisión de bits del algoritmo para elegir el número de qubits, y tendrás que diseñar los operadores de exponenciación modular $M_a$. Te animamos a que pruebes esto por ti mismo, y luego leas sobre las metodologías mostradas en la Fig. 9 de [6] y la Fig. 2 de [7].

def M_a_mod21():
    """
    M_a (mod 21)
    """

    # Your code here
    pass

Referencias

1. Shor, Peter W. "Algoritmos de tiempo polinómico para la factorización de primos y logaritmos discretos en un ordenador cuántico " Revisión SIAM 41.2 (1999): 303-332.
2. IBM Quantum Curso «Fundamentos de los algoritmos cuánticos», impartido por el Dr. John Watrous.
3. Vandersypen, Lieven MK, et al. "Experimental realization of Shor's quantum factoring algorithm using nuclear magnetic resonance " Nature 414.6866 (2001): 883-887.
4. Martín-López, Enrique, et al. "Experimental realization of Shor's quantum factoring algorithm using qubit recycling" Nature photonics 6.11 (2012): 773-776.
5. Monz, Thomas, et al. "Realización de un algoritmo Shor escalable " Science 351.6277 (2016): 1068-1070.
6. Amico, Mirko, Zain H. Saleem y Muir Kumph. "Estudio experimental del algoritmo de factorización de Shor utilizando la experiencia IBM Q " Physical Review A 100.1 (2019): 012305.
7. Skosana, Unathi y Mark Tame. "Demostración del algoritmo de factorización de Shor para N=21 en procesadores cuánticos IBM " Informes científicos 11.1 (2021): 16599.
8. Kitaev, A. Yu. "Medidas cuánticas y el problema del estabilizador abeliano " arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig, y Martin Ekerå. "Cómo factorizar enteros RSA de 2048 bits en 8 horas usando 20 millones de qubits ruidosos " Quantum 5 (2021): 433.