Algoritmo de Shor

Estimación de uso: Tres segundos en un procesador Eagle r3 (NOTA: Esto es solo una estimación. Su tiempo de ejecución podría variar.)

Resultados de aprendizaje

Tras completar este tutorial, los usuarios deberán comprender:

• El contexto matemático del algoritmo de Shor para la factorización de enteros
• Cómo ejecutar una instancia de ejemplo de este algoritmo en hardware

Requisitos previos

Sugerimos que los usuarios estén familiarizados con los siguientes temas antes de completar este tutorial:

• Fundamentos de algoritmos cuánticos.
• Estimación de fase y factorización. Cubrimos parte de este material en este tutorial.

Antecedentes

El algoritmo de Shor, desarrollado por Peter Shor en 1994, es un algoritmo cuántico revolucionario para factorizar enteros en tiempo polinomial. Su importancia radica en su capacidad para factorizar enteros grandes exponencialmente más rápido que cualquier algoritmo clásico conocido, amenazando la seguridad de los sistemas criptográficos ampliamente utilizados como RSA, que dependen de la dificultad de factorizar números grandes. Al resolver eficientemente este problema en una computadora cuántica suficientemente potente, el algoritmo de Shor podría revolucionar campos como la criptografía, la ciberseguridad y las matemáticas computacionales, subrayando el poder transformador de la computación cuántica.

Este tutorial se centra en demostrar el algoritmo de Shor factorizando 15 en una computadora cuántica.

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

Finalizamos 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 diseñadas específicamente para factorizar enteros específicos como 15 y 21. Nota: Este tutorial se centra más en la implementación y demostración de los circuitos relacionados con el algoritmo de Shor. Para un recurso educativo en profundidad sobre el material, consulta el curso Fundamentals of quantum algorithms del Dr. John Watrous, y los artículos en la sección de Referencias.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener 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

# Added by doQumentation — required packages for this notebook
!pip install -q numpy pandas qiskit qiskit-ibm-runtime

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: Mapear entradas clásicas a un problema cuántico

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

Problema de estimación de fase

En el problema de estimación de fase, se nos proporciona un estado cuántico $\ket{\psi}$ de $n$ qubits, junto con un circuito cuántico unitario que actúa sobre $n$ qubits. Se nos garantiza 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 cual corresponde $\ket{\psi}$. En otras palabras, el circuito debería producir una aproximación del número $\theta \in [0, 1)$ que satisface $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, deseamos 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 medición sobre $m$ qubits. En el circuito anterior, los $m$ qubits superiores se inicializan en el estado $\ket{0^m}$, y los $n$ qubits inferiores se inicializan en $\ket{\psi}$, el cual se garantiza que es un vector propio de $U$. El primer ingrediente en el circuito de estimación de fase son las operaciones unitarias controladas que se encargan de realizar un retroceso de fase (phase kickback) en su qubit de control correspondiente. Estas unitarias controladas se exponencian de acuerdo con la posición del qubit de control, desde el bit menos significativo hasta el bit más significativo. Dado que $\ket{\psi}$ es un vector propio de $U$, el estado de los $n$ qubits inferiores no se ve afectado por esta operación, pero la información de fase del valor propio se propaga a los $m$ qubits superiores. Resulta que después de la operación de retroceso de fase mediante unitarias controladas, todos los estados posibles de los $m$ qubits superiores son ortonormales entre sí para cada vector propio $\ket{\psi}$ del operador unitario $U$. Por lo tanto, estos estados son perfectamente distinguibles, y podemos rotar la base que forman de vuelta a la base computacional para realizar 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 de dimensión $2^m$. 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 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 Fundamentals of quantum algorithms. Ahora estamos listos para utilizar el circuito de estimación de fase para la búsqueda de orden.

Problema de búsqueda de orden

Para definir el problema de búsqueda de orden, comenzamos con algunos conceptos de teoría de números. Primero, para cualquier entero positivo dado $N$, definimos 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 módulo $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 entero positivo más pequeño $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. Primero, necesitamos definir el operador unitario $U$ que nos permitirá encontrar el orden $r$, y segundo, 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$. Observa que está implícito que estamos tomando el producto módulo $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 propios y valores propios que nos permiten conectar el orden $r$ de $a$ con el problema de estimación de fase. Específicamente, 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 conveniente de vector propio/valor propio 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 así obtener una estimación del orden $r$. Sin embargo, no es fácil hacerlo, y necesitamos considerar una alternativa.

Consideremos qué resultado produciría el circuito si preparamos el estado computacional $\ket{1}$ como estado inicial. Este no es un estado propio de $M_a$, pero es 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 establecemos el estado inicial como $\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 hubiéramos usado $\ket{\psi_k}$ como 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 después de varias ejecuciones independientes, que era nuestro objetivo.

Operadores de exponenciación modular

Hasta ahora, vinculamos el problema de estimación de fase con el 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 restante es encontrar una forma eficiente de definir las 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 solo necesitamos las potencias que son a su vez potencias de 2, podemos hacer esto clásicamente de forma eficiente usando la elevación al cuadrado iterativa.

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

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

Podemos hacer una pausa aquí para discutir un ejemplo específico y construir el circuito de búsqueda de orden para $N=15$. Observa que los posibles valores no triviales de $a \in \mathbb{Z}_N^*$ 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 computacional 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 de la base se reorganizan, por lo que tenemos una matriz de permutación. Podemos construir esta operación en cuatro qubits con compuertas de intercambio (swap). 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)

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)

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

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

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 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)

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

Como podemos ver en la lista de valores de $b$, además de $M_2$ que construimos previamente, también necesitamos construir $M_4$ y $M_1$. Observa 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 de la base computacional 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 puede construirse 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)

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)

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

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

Vimos 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í, dado que $N=15$ requiere solo cuatro qubits, pudimos sintetizar estas operaciones directamente con compuertas SWAP por inspección. En general, este podría no ser un enfoque escalable. En su lugar, podríamos necesitar construir la matriz de permutación explícitamente y usar la clase UnitaryGate de Qiskit y los métodos de transpilación para sintetizar esta matriz de permutación. Sin embargo, esto puede resultar en circuitos significativamente más profundos. A continuación se muestra 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
)

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 conteos con la profundidad del circuito compilado de nuestra implementación manual de la compuerta $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
)

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 más profundo incluso para una sola compuerta $M_2$ en comparación con nuestra implementación manual. Por lo tanto, continuaremos con nuestra implementación anterior de las operaciones $M_b$. Ahora estamos listos para construir el circuito completo de búsqueda de orden utilizando nuestros operadores de exponenciación modular controlada previamente definidos. En el siguiente código, también importamos el circuito QFT de la biblioteca de circuitos de Qiskit, que utiliza compuertas Hadamard en cada qubit, una serie de compuertas controladas U1 (o Z, dependiendo de la fase) y una capa de compuertas de intercambio.

# 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)

Observa que omitimos las operaciones de exponenciación modular controlada de los qubits de control restantes porque $M_1$ es el operador identidad. Ten en cuenta que más adelante en este tutorial, ejecutaremos este circuito en el backend ibm_marrakesh. Para ello, transpilamos el circuito de acuerdo con este backend específico e informamos la profundidad del circuito y los conteos de compuertas.

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
)

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

Paso 3: Ejecutar utilizando primitivas de Qiskit

Primero, discutimos lo que teóricamente obtendríamos si ejecutáramos este circuito en un simulador ideal. A continuación, tenemos un conjunto de resultados de simulación del circuito anterior utilizando 1024 disparos. Como podemos ver, obtenemos una distribución aproximadamente uniforme sobre cuatro cadenas de bits en los qubits de control.

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

plot_histogram(counts)

Al medir 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 encontrar 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 de la siguiente manera.

# 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)

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

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

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Debido a que esto produce fracciones que devuelven el resultado exacto (en este caso, 0.6660000...), puede dar resultados complicados como el anterior. Podemos usar el método .limit_denominator() para obtener la fracción que más se asemeja a nuestro número de punto flotante, con un denominador por debajo de cierto valor:

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

Fraction(2, 3)

Esto es mucho mejor. El orden (r) debe ser menor que N, por lo que estableceremos 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)

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 proporcionaron el resultado correcto: $r=4$, y podemos observar que el algoritmo de Shor para la búsqueda de orden tiene una posibilidad de fallar. Estos resultados incorrectos se deben a que $k = 0$, o porque $k$ y $r$ no son coprimos, y en lugar de $r$, obtenemos un factor de $r$. La solución más sencilla a esto es simplemente repetir el experimento hasta obtener un resultado satisfactorio para $r$. Hasta ahora, implementamos 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 será relacionar el problema de búsqueda de orden con el problema de factorización de enteros. Esta última parte del algoritmo es puramente clásica y puede resolverse en una computadora clásica después de que las mediciones de fase se hayan obtenido de una computadora cuántica. Por lo tanto, diferimos la última parte del algoritmo hasta después de demostrar cómo podemos ejecutar el circuito de búsqueda de orden en hardware real.

Ejecuciones en hardware

Ahora podemos ejecutar el circuito de búsqueda de orden que previamente transpilamos para ibm_marrakesh. Aquí recurrimos al desacoplamiento dinámico (DD) para la supresión de errores, y al gate twirling para propósitos de mitigación de errores. El DD implica aplicar secuencias de pulsos de control cronometrados con precisión a un dispositivo cuántico, promediando efectivamente las interacciones ambientales no deseadas y la decoherencia. El gate twirling, por otro lado, aleatoriza compuertas cuánticas específicas para transformar errores coherentes en errores de Pauli, que se acumulan linealmente en lugar de cuadráticamente. Ambas técnicas se combinan frecuentemente para mejorar la coherencia y la 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))

Como podemos ver, obtuvimos las mismas cadenas de bits con los conteos más altos. Dado que el hardware cuántico tiene ruido, existe cierta fuga hacia 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)

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

Paso 4: Post-procesar y devolver el resultado en el formato clásico deseado

Factorización de enteros

Hasta ahora, discutimos 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 con la factorización de enteros, lo que completa el algoritmo de Shor. Ten en cuenta que esta parte del algoritmo es clásica. Ahora demostramos esto usando nuestro ejemplo de $N = 15$ y $a = 2$. Recuerda que la fase que medimos es $k / r$, donde $a^r \; (\textrm{mod} \; N) = 1$ y $k$ es un entero aleatorio entre $0$ y $r - 1$. 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 continuar y debemos intentar de nuevo con un valor diferente para $a$; de lo contrario, existe una alta probabilidad de que el máximo común divisor de $N$ y ya sea $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 encuentra al menos un factor de $N$. La celda a continuación repite el algoritmo hasta que se encuentra al menos un factor de $N=15$. Utilizaremos los resultados de la ejecución en hardware anterior para estimar 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} ***")

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 ***

Discusión

Trabajo relacionado

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

El trabajo seminal [3] de IBM® demostró el algoritmo de Shor por primera vez, factorizando el número 15 en sus factores primos 3 y 5 utilizando una computadora cuántica de resonancia magnética nuclear (RMN) de siete qubits. Otro experimento [4] factorizó 15 utilizando qubits fotónicos. Al emplear un solo qubit reciclado múltiples veces y codificar el registro de trabajo en estados de mayor dimensión, los investigadores redujeron el número requerido de qubits a un tercio 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 a mitad de circuito con operaciones de realimentación y reciclaje de qubits con operaciones de reinicio. Esta demostración se realizó en una computadora cuántica de trampa de iones.

Un trabajo más reciente [6] se centró en factorizar 15, 21 y 35 en hardware IBM Quantum®. De manera similar a trabajos anteriores, los investigadores utilizaron una versión compilada del algoritmo que empleó una transformada cuántica de Fourier semi-clásica como la propuesta por Kitaev para minimizar el número de qubits físicos y compuertas. 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 involucró el uso de una versión compilada de la rutina de estimación de fase cuántica, y se basó en la demostración previa de [4]. Los autores fueron más allá de este trabajo al usar una configuración de compuertas Toffoli aproximadas con desplazamientos de fase residuales. El algoritmo fue implementado en procesadores cuánticos de IBM utilizando solo cinco qubits, y la presencia de entrelazamiento entre los qubits de control y de registro fue verificada exitosamente.

Escalamiento del algoritmo

Observamos que el cifrado RSA típicamente involucra tamaños de clave del orden de 2048 a 4096 bits. Intentar factorizar un número de 2048 bits con el algoritmo de Shor resultará en un circuito cuántico con millones de qubits, incluyendo la sobrecarga de corrección de errores y una profundidad de circuito del orden de mil millones, lo cual está más allá de los límites del hardware cuántico actual para su ejecución. Por lo tanto, el algoritmo de Shor requerirá ya sea métodos optimizados de construcción de circuitos o una corrección de errores cuántica robusta para ser prácticamente viable en la ruptura de sistemas criptográficos modernos. Te remitimos a [9] para una discusión más detallada sobre la estimación de recursos para el algoritmo de Shor.

Desafío

Felicitaciones por completar el tutorial. Ahora es un buen momento para poner a prueba tu comprensión. ¿Podrías intentar construir el circuito para factorizar 21? Puedes seleccionar un valor de $a$ de tu elección. Necesitarás decidir la precisión en bits del algoritmo para elegir el número de qubits, y necesitarás diseñar los operadores de exponenciación modular $M_a$. Te animamos a que lo intentes por tu cuenta, y luego lee 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. "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer." SIAM review 41.2 (1999): 303-332.
2. IBM Quantum Fundamentals of Quantum Algorithms course by 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. Martin-Lopez, 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. "Realization of a scalable Shor algorithm." Science 351.6277 (2016): 1068-1070.
6. Amico, Mirko, Zain H. Saleem, and Muir Kumph. "Experimental study of Shor's factoring algorithm using the IBM Q Experience." Physical Review A 100.1 (2019): 012305.
7. Skosana, Unathi, and Mark Tame. "Demonstration of Shor's factoring algorithm for N=21 on IBM quantum processors." Scientific reports 11.1 (2021): 16599.
8. Kitaev, A. Yu. "Quantum measurements and the Abelian stabilizer problem." arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig, and Martin Ekerå. "How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits." Quantum 5 (2021): 433.