Algoritmo di Shor

Stima di utilizzo: Tre secondi su un processore Eagle r3 (NOTA: Questa è solo una stima. Il vostro tempo di esecuzione potrebbe variare.)

L'algoritmo di Shor, sviluppato da Peter Shor nel 1994, è un algoritmo quantistico rivoluzionario per la fattorizzazione di interi in tempo polinomiale. La sua importanza risiede nella capacità di fattorizzare grandi interi in modo esponenzialmente più veloce rispetto a qualsiasi algoritmo classico conosciuto, minacciando la sicurezza di sistemi crittografici ampiamente utilizzati come RSA, che si basano sulla difficoltà di fattorizzare numeri grandi. Risolvendo efficacemente questo problema su un computer quantistico sufficientemente potente, l'algoritmo di Shor potrebbe rivoluzionare campi come la crittografia, la sicurezza informatica e la matematica computazionale, sottolineando il potere trasformativo del calcolo quantistico.

Questo tutorial si concentra sulla dimostrazione dell'algoritmo di Shor fattorizzando 15 su un computer quantistico.

Prima definiamo il problema della ricerca dell'ordine e costruiamo i circuiti corrispondenti dal protocollo di stima della fase quantistica. Successivamente, eseguiamo i circuiti di ricerca dell'ordine su hardware reale utilizzando i circuiti di profondità minima che possiamo traspilare. L'ultima sezione completa l'algoritmo di Shor collegando il problema della ricerca dell'ordine alla fattorizzazione di interi.

Concludiamo il tutorial con una discussione su altre dimostrazioni dell'algoritmo di Shor su hardware reale, concentrandoci sia su implementazioni generiche sia su quelle adattate per fattorizzare interi specifici come 15 e 21. Nota: Questo tutorial si concentra maggiormente sull'implementazione e la dimostrazione dei circuiti riguardanti l'algoritmo di Shor. Per una risorsa educativa approfondita sul materiale, si prega di fare riferimento al corso Fundamentals of quantum algorithms del Dr. John Watrous, e agli articoli nella sezione Riferimenti.

Requisiti

Prima di iniziare questo tutorial, assicuratevi di avere installato quanto segue:

• Qiskit SDK v2.0 o successivo, con supporto per la visualizzazione
• Qiskit Runtime v0.40 o successivo (pip install qiskit-ibm-runtime)

Configurazione

# 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

Passo 1: Mappare gli input classici a un problema quantistico

Contesto

L'algoritmo di Shor per la fattorizzazione di interi utilizza un problema intermediario noto come problema della ricerca dell'ordine. In questa sezione, dimostriamo come risolvere il problema della ricerca dell'ordine utilizzando la stima della fase quantistica.

Problema della stima della fase

Nel problema della stima della fase, ci viene dato uno stato quantistico $\ket{\psi}$ di $n$ qubit, insieme a un circuito quantistico unitario che agisce su $n$ qubit. Ci viene promesso che $\ket{\psi}$ è un autovettore della matrice unitaria $U$ che descrive l'azione del circuito, e il nostro obiettivo è calcolare o approssimare l'autovalore $\lambda = e^{2 \pi i \theta}$ a cui $\ket{\psi}$ corrisponde. In altre parole, il circuito dovrebbe restituire un'approssimazione del numero $\theta \in [0, 1)$ che soddisfa $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ L'obiettivo del circuito di stima della fase è approssimare $\theta$ in $m$ bit. Matematicamente parlando, vorremmo trovare $y$ tale che $\theta \approx y / 2^m$, dove $y \in {0, 1, 2, \dots, 2^{m-1}}$. L'immagine seguente mostra il circuito quantistico che stima $y$ in $m$ bit effettuando una misura su $m$ qubit. Nel circuito sopra, i $m$ qubit superiori sono inizializzati nello stato $\ket{0^m}$, e gli $n$ qubit inferiori sono inizializzati in $\ket{\psi}$, che è promesso essere un autovettore di $U$. Il primo ingrediente nel circuito di stima della fase sono le operazioni unitarie controllate che sono responsabili di eseguire un phase kickback al loro qubit di controllo corrispondente. Queste unitarie controllate sono esponenziate in accordo alla posizione del qubit di controllo, variando dal bit meno significativo al bit più significativo. Poiché $\ket{\psi}$ è un autovettore di $U$, lo stato degli $n$ qubit inferiori non è influenzato da questa operazione, ma l'informazione di fase dell'autovalore si propaga ai $m$ qubit superiori. Si scopre che dopo l'operazione di phase kickback tramite unitarie controllate, tutti i possibili stati dei $m$ qubit superiori sono ortonormali tra loro per ogni autovettore $\ket{\psi}$ dell'unitaria $U$. Pertanto, questi stati sono perfettamente distinguibili, e possiamo ruotare la base che formano indietro alla base computazionale per effettuare una misura. Un'analisi matematica mostra che questa matrice di rotazione corrisponde alla trasformata di Fourier quantistica (QFT) inversa nello spazio di Hilbert a $2^m$ dimensioni. L'intuizione dietro questo è che la struttura periodica degli operatori di esponenziazione modulare è codificata nello stato quantistico, e la QFT converte questa periodicità in picchi misurabili nel dominio delle frequenze.

Per una comprensione più approfondita del motivo per cui il circuito QFT è impiegato nell'algoritmo di Shor, rimandiamo il lettore al corso Fundamentals of quantum algorithms. Siamo ora pronti a utilizzare il circuito di stima della fase per la ricerca dell'ordine.

Problema della ricerca dell'ordine

Per definire il problema della ricerca dell'ordine, iniziamo con alcuni concetti di teoria dei numeri. Prima di tutto, per qualsiasi intero positivo dato $N$, definiamo l'insieme $\mathbb{Z}_N$ come $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Tutte le operazioni aritmetiche in $\mathbb{Z}_N$ sono eseguite modulo $N$. In particolare, tutti gli elementi $a \in \mathbb{Z}_n$ che sono coprimi con $N$ sono speciali e costituiscono $\mathbb{Z}^*_N$ come $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Per un elemento $a \in \mathbb{Z}^*_N$, il più piccolo intero positivo $r$ tale che $a^r \equiv 1 \; (\mathrm{mod} \; N)$ è definito come l'ordine di $a$ modulo $N$. Come vedremo più avanti, trovare l'ordine di un $a \in \mathbb{Z}^*_N$ ci permetterà di fattorizzare $N$. Per costruire il circuito di ricerca dell'ordine dal circuito di stima della fase, abbiamo bisogno di due considerazioni. Prima, dobbiamo definire l'unitaria $U$ che ci permetterà di trovare l'ordine $r$, e secondo, dobbiamo definire un autovettore $\ket{\psi}$ di $U$ per preparare lo stato iniziale del circuito di stima della fase.

Per collegare il problema della ricerca dell'ordine alla stima della fase, consideriamo l'operazione definita su un sistema i cui stati classici corrispondono a $\mathbb{Z}_N$, dove moltiplichiamo per un elemento fisso $a \in \mathbb{Z}^*_N$. In particolare, definiamo questo operatore di moltiplicazione $M_a$ tale che $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ per ogni $x \in \mathbb{Z}_N$. Notate che è implicito che stiamo prendendo il prodotto modulo $N$ all'interno del ket sul lato destro dell'equazione. Un'analisi matematica mostra che $M_a$ è un operatore unitario. Inoltre, si scopre che $M_a$ ha coppie autovettore-autovalore che ci permettono di collegare l'ordine $r$ di $a$ al problema della stima della fase. Specificamente, per qualsiasi scelta di $j \in \{0, \dots, r-1\}$, abbiamo che $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ è un autovettore di $M_a$ il cui autovalore corrispondente è $\omega^{j}_{r}$, dove $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Per osservazione, vediamo che una coppia autovettore/autovalore conveniente è lo stato $\ket{\psi_1}$ con $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Pertanto, se potessimo trovare l'autovettore $\ket{\psi_1}$, potremmo stimare la fase $\theta=1/r$ con il nostro circuito quantistico e quindi ottenere una stima dell'ordine $r$. Tuttavia, non è facile farlo, e dobbiamo considerare un'alternativa.

Consideriamo cosa risulterebbe dal circuito se preparassimo lo stato computazionale $\ket{1}$ come stato iniziale. Questo non è un autostato di $M_a$, ma è la sovrapposizione uniforme degli autostati che abbiamo appena descritto sopra. In altre parole, vale la seguente relazione. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ L'implicazione dell'equazione sopra è che se impostiamo lo stato iniziale a $\ket{1}$, otterremo precisamente lo stesso risultato di misura come se avessimo scelto $k \in \{ 0, \dots, r-1\}$ uniformemente a caso e usato $\ket{\psi_k}$ come autovettore nel circuito di stima della fase. In altre parole, una misura dei $m$ qubit superiori produce un'approssimazione $y / 2^m$ al valore $k / r$ dove $k \in \{ 0, \dots, r-1\}$ è scelto uniformemente a caso. Questo ci permette di apprendere $r$ con un alto grado di confidenza dopo diverse esecuzioni indipendenti, che era il nostro obiettivo.

Operatori di esponenziazione modulare

Finora, abbiamo collegato il problema della stima della fase al problema della ricerca dell'ordine definendo $U = M_a$ e $\ket{\psi} = \ket{1}$ nel nostro circuito quantistico. Pertanto, l'ultimo ingrediente rimanente è trovare un modo efficiente per definire esponenziali modulari di $M_a$ come $M_a^k$ per $k = 1, 2, 4, \dots, 2^{m-1}$. Per eseguire questo calcolo, scopriamo che per qualsiasi potenza $k$ scegliamo, possiamo creare un circuito per $M_a^k$ non iterando $k$ volte il circuito per $M_a$, ma invece calcolando $b = a^k \; \mathrm{mod} \; N$ e poi usando il circuito per $M_b$. Poiché abbiamo bisogno solo delle potenze che sono esse stesse potenze di 2, possiamo farlo in modo classicamente efficiente usando l'elevamento al quadrato iterativo.

Passo 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Esempio specifico con $N = 15$ e $a=2$

Possiamo fare una pausa qui per discutere un esempio specifico e costruire il circuito di ricerca dell'ordine per $N=15$. Notate che i possibili $a \in \mathbb{Z}_N^*$ non banali per $N=15$ sono $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Per questo esempio, scegliamo $a=2$. Costruiremo l'operatore $M_2$ e gli operatori di esponenziazione modulare $M_2^k$. L'azione di $M_2$ sugli stati della base computazionale è la seguente. $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}$ Per osservazione, possiamo vedere che gli stati di base sono mescolati, quindi abbiamo una matrice di permutazione. Possiamo costruire questa operazione su quattro qubit con porte swap. Di seguito, costruiamo l'operazione $M_2$ e la sua versione controllata.

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)

Le porte che agiscono su più di due qubit saranno ulteriormente decomposte in porte a due qubit.

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

Ora dobbiamo costruire gli operatori di esponenziazione modulare. Per ottenere una precisione sufficiente nella stima della fase, utilizzeremo otto qubit per la misura di stima. Pertanto, dobbiamo costruire $M_b$ con $b = a^{2^k} \; (\mathrm{mod} \; N)$ per ogni $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]

Come possiamo vedere dalla lista dei valori di $b$, oltre a $M_2$ che abbiamo costruito precedentemente, dobbiamo anche costruire $M_4$ e $M_1$. Notate che $M_1$ agisce in modo banale sugli stati della base computazionale, quindi è semplicemente l'operatore identità.

$M_4$ agisce sugli stati della base computazionale come segue. $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}$

Pertanto, questa permutazione può essere costruita con la seguente operazione swap.

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)

Le porte che agiscono su più di due qubit saranno ulteriormente decomposte in porte a due qubit.

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

Abbiamo visto che gli operatori $M_b$ per un dato $b \in \mathbb{Z}^*_N$ sono operazioni di permutazione. A causa della dimensione relativamente piccola del problema di permutazione che abbiamo qui, poiché $N=15$ richiede solo quattro qubit, siamo stati in grado di sintetizzare queste operazioni direttamente con porte SWAP tramite ispezione. In generale, questo potrebbe non essere un approccio scalabile. Invece, potremmo dover costruire la matrice di permutazione esplicitamente, e usare la classe UnitaryGate di Qiskit e i metodi di traspilazione per sintetizzare questa matrice di permutazione. Tuttavia, questo può risultare in circuiti significativamente più profondi. Segue un esempio.

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

Confrontiamo questi conteggi con la profondità del circuito compilato della nostra implementazione manuale della porta $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})

Come possiamo vedere, l'approccio della matrice di permutazione ha prodotto un circuito significativamente profondo anche per una singola porta $M_2$ rispetto alla nostra implementazione manuale. Pertanto, continueremo con la nostra precedente implementazione delle operazioni $M_b$. Ora siamo pronti a costruire il circuito completo di ricerca dell'ordine utilizzando i nostri operatori di esponenziazione modulare controllati precedentemente definiti. Nel codice seguente, importiamo anche il circuito QFT dalla libreria di circuiti di Qiskit, che utilizza porte Hadamard su ogni qubit, una serie di porte controllate-U1 (o Z, a seconda della fase), e uno strato di porte 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)

Notate che abbiamo omesso le operazioni di esponenziazione modulare controllate dai qubit di controllo rimanenti perché $M_1$ è l'operatore identità. Notate che più avanti in questo tutorial, eseguiremo questo circuito sul backend ibm_marrakesh. Per fare questo, traspiliamo il circuito in accordo a questo backend specifico e riportiamo la profondità del circuito e i conteggi delle porte.

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

Passo 3: Eseguire utilizzando i primitive Qiskit

Innanzitutto, discutiamo cosa otterremmo teoricamente se eseguissimo questo circuito su un simulatore ideale. Di seguito, abbiamo una serie di risultati di simulazione del circuito sopra utilizzando 1024 shot. Come possiamo vedere, otteniamo una distribuzione approssimativamente uniforme su quattro bitstring sui qubit di controllo.

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

plot_histogram(counts)

Misurando i qubit di controllo, otteniamo una stima di fase a otto bit dell'operatore $M_a$. Possiamo convertire questa rappresentazione binaria in decimale per trovare la fase misurata. Come possiamo vedere dall'istogramma sopra, sono state misurate quattro bitstring diverse, e ciascuna di esse corrisponde a un valore di fase come segue.

# 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

Ricordiamo che qualsiasi fase misurata corrisponde a $\theta = k / r$ dove $k$ è campionato uniformemente a caso da $\{0, 1, \dots, r-1 \}$. Pertanto, possiamo utilizzare l'algoritmo delle frazioni continue per tentare di trovare $k$ e l'ordine $r$. Python ha questa funzionalità integrata. Possiamo usare il modulo fractions per trasformare un float in un oggetto Fraction, per esempio:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Poiché questo fornisce frazioni che restituiscono il risultato esattamente (in questo caso, 0.6660000...), può dare risultati complicati come quello sopra. Possiamo usare il metodo .limit_denominator() per ottenere la frazione che assomiglia più da vicino al nostro float, con un denominatore inferiore a un certo valore:

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

Fraction(2, 3)

Questo è molto più chiaro. L'ordine (r) deve essere inferiore a N, quindi imposteremo il denominatore massimo su 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

Possiamo vedere che due degli autovalori misurati ci hanno fornito il risultato corretto: $r=4$, e possiamo vedere che l'algoritmo di Shor per trovare l'ordine ha una probabilità di fallire. Questi risultati negativi sono dovuti al fatto che $k = 0$, o perché $k$ e $r$ non sono coprimi - e invece di $r$, ci viene dato un fattore di $r$. La soluzione più semplice a questo è semplicemente ripetere l'esperimento finché non otteniamo un risultato soddisfacente per $r$. Finora, abbiamo implementato il problema di ricerca dell'ordine per $N=15$ con $a=2$ utilizzando il circuito di stima di fase su un simulatore. L'ultimo passo dell'algoritmo di Shor sarà quello di mettere in relazione il problema di ricerca dell'ordine con il problema di fattorizzazione di interi. Quest'ultima parte dell'algoritmo è puramente classica e può essere risolta su un computer classico dopo che le misurazioni di fase sono state ottenute da un computer quantistico. Pertanto, rimandiamo l'ultima parte dell'algoritmo fino a dopo aver dimostrato come possiamo eseguire il circuito di ricerca dell'ordine su hardware reale.

Esecuzioni su hardware

Ora possiamo eseguire il circuito di ricerca dell'ordine che abbiamo precedentemente traspilato per ibm_marrakesh. Qui utilizziamo il disaccoppiamento dinamico (DD) per la soppressione degli errori, e al gate twirling per scopi di mitigazione degli errori. Il DD comporta l'applicazione di sequenze di impulsi di controllo temporizzati con precisione a un dispositivo quantistico, mediando efficacemente le interazioni ambientali indesiderate e la decoerenza. Il gate twirling, d'altra parte, randomizza specifici gate quantistici per trasformare gli errori coerenti in errori di Pauli, che si accumulano linearmente piuttosto che quadraticamente. Entrambe le tecniche sono spesso combinate per migliorare la coerenza e la fedeltà dei calcoli quantistici.

# 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

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

Come possiamo vedere, abbiamo ottenuto le stesse bitstring con i conteggi più alti. Poiché l'hardware quantistico ha rumore, c'è una certa dispersione verso altre bitstring, che possiamo filtrare statisticamente.

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

Passo 4: Post-elaborare e restituire il risultato in formato classico desiderato

Fattorizzazione di interi

Finora, abbiamo discusso come possiamo implementare il problema di ricerca dell'ordine utilizzando un circuito di stima di fase. Ora, colleghiamo il problema di ricerca dell'ordine alla fattorizzazione di interi, che completa l'algoritmo di Shor. Si noti che questa parte dell'algoritmo è classica. Ora dimostriamo questo utilizzando il nostro esempio di $N = 15$ e $a = 2$. Ricordiamo che la fase che abbiamo misurato è $k / r$, dove $a^r \; (\textrm{mod} \; N) = 1$ e $k$ è un intero casuale tra $0$ e $r - 1$. Da questa equazione, abbiamo $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ il che significa che $N$ deve dividere $a^r-1$. Se $r$ è anche pari, allora possiamo scrivere $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Se $r$ non è pari, non possiamo andare oltre e dobbiamo riprovare con un valore diverso per $a$; altrimenti, c'è un'alta probabilità che il massimo comun divisore di $N$ e di $a^{r/2}-1$, o $a^{r/2}+1$ sia un fattore proprio di $N$.

Poiché alcune esecuzioni dell'algoritmo falliranno statisticamente, ripeteremo questo algoritmo finché non verrà trovato almeno un fattore di $N$. La cella seguente ripete l'algoritmo finché non viene trovato almeno un fattore di $N=15$. Utilizzeremo i risultati dell'esecuzione hardware sopra per indovinare la fase e il fattore corrispondente in ciascuna iterazione.

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

Discussione

Lavori correlati

In questa sezione, discutiamo altri lavori importanti che hanno dimostrato l'algoritmo di Shor su hardware reale.

Il lavoro seminale [3] di IBM® ha dimostrato l'algoritmo di Shor per la prima volta, fattorizzando il numero 15 nei suoi fattori primi 3 e 5 utilizzando un computer quantistico a risonanza magnetica nucleare (NMR) a sette qubit. Un altro esperimento [4] ha fattorizzato 15 utilizzando qubit fotonici. Impiegando un singolo qubit riciclato più volte e codificando il registro di lavoro in stati a dimensione superiore, i ricercatori hanno ridotto il numero di qubit richiesti a un terzo rispetto al protocollo standard, utilizzando un algoritmo compilato a due fotoni. Un articolo significativo nella dimostrazione dell'algoritmo di Shor è [5], che utilizza la tecnica di stima di fase iterativa di Kitaev [8] per ridurre il requisito di qubit dell'algoritmo. Gli autori hanno utilizzato sette qubit di controllo e quattro qubit cache, insieme all'implementazione di moltiplicatori modulari. Questa implementazione, tuttavia, richiede misurazioni a metà circuito con operazioni di feed-forward e riciclaggio di qubit con operazioni di reset. Questa dimostrazione è stata eseguita su un computer quantistico a trappola ionica.

Lavori più recenti [6] si sono concentrati sulla fattorizzazione di 15, 21 e 35 su hardware IBM Quantum®. Similmente ai lavori precedenti, i ricercatori hanno utilizzato una versione compilata dell'algoritmo che impiegava una trasformata di Fourier quantistica semi-classica come proposto da Kitaev per minimizzare il numero di qubit fisici e gate. Un lavoro più recente [7] ha anche eseguito una dimostrazione proof-of-concept per la fattorizzazione dell'intero 21. Questa dimostrazione ha anche coinvolto l'uso di una versione compilata della routine di stima di fase quantistica, e si è basata sulla dimostrazione precedente di [4]. Gli autori sono andati oltre questo lavoro utilizzando una configurazione di gate Toffoli approssimativi con sfasamenti di fase residui. L'algoritmo è stato implementato su processori quantistici IBM utilizzando solo cinque qubit, e la presenza di entanglement tra i qubit di controllo e di registro è stata verificata con successo.

Scalabilità dell'algoritmo

Notiamo che la crittografia RSA tipicamente coinvolge dimensioni di chiave dell'ordine di 2048-4096 bit. Il tentativo di fattorizzare un numero a 2048 bit con l'algoritmo di Shor risulterà in un circuito quantistico con milioni di qubit, includendo l'overhead di correzione degli errori e una profondità del circuito dell'ordine di un miliardo, che è oltre i limiti dell'hardware quantistico attuale per l'esecuzione. Pertanto, l'algoritmo di Shor richiederà metodi di costruzione del circuito ottimizzati o una robusta correzione degli errori quantistici per essere praticamente praticabile per violare i moderni sistemi crittografici. Vi rimandiamo a [9] per una discussione più dettagliata sulla stima delle risorse per l'algoritmo di Shor.

Sfida

Congratulazioni per aver completato il tutorial! Ora è un ottimo momento per testare la vostra comprensione. Potreste provare a costruire il circuito per fattorizzare 21? Potete selezionare un $a$ di vostra scelta. Dovrete decidere sulla precisione in bit dell'algoritmo per scegliere il numero di qubit, e dovrete progettare gli operatori di esponenziazione modulare $M_a$. Vi incoraggiamo a provarlo voi stessi, e poi leggere le metodologie mostrate in Fig. 9 di [6] e Fig. 2 di [7].

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

    # Your code here
    pass

References

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 Ekera. "How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits." Quantum 5 (2021): 433.