Algoritmul lui Shor

Estimare de utilizare: Trei secunde pe un procesor Eagle r3 (NOTĂ: Aceasta este doar o estimare. Timpul tău de execuție poate varia.)

Algoritmul lui Shor, dezvoltat de Peter Shor în 1994, este un algoritm cuantic revoluționar pentru factorizarea numerelor întregi în timp polinomial. Semnificația sa constă în capacitatea de a factoriza numere întregi mari exponențial mai rapid decât orice algoritm clasic cunoscut, amenințând securitatea sistemelor criptografice utilizate pe scară largă, cum ar fi RSA, care se bazează pe dificultatea factorizării numerelor mari. Prin rezolvarea eficientă a acestei probleme pe un calculator cuantic suficient de puternic, algoritmul lui Shor ar putea revoluționa domenii precum criptografia, securitatea cibernetică și matematica computațională, subliniind puterea transformatoare a calculului cuantic.

Acest tutorial se concentrează pe demonstrarea algoritmului lui Shor prin factorizarea lui 15 pe un calculator cuantic.

Mai întâi, definim problema găsirii ordinului și construim circuitele corespunzătoare din protocolul de estimare a fazei cuantice. Apoi, rulăm circuitele de găsire a ordinului pe hardware real, folosind circuitele cu adâncimea cea mai mică pe care le putem transpila. Ultima secțiune finalizează algoritmul lui Shor conectând problema găsirii ordinului la factorizarea numerelor întregi.

Încheiem tutorialul cu o discuție despre alte demonstrații ale algoritmului lui Shor pe hardware real, concentrându-ne atât pe implementările generice, cât și pe cele adaptate pentru factorizarea unor numere întregi specifice, cum ar fi 15 și 21. Notă: Acest tutorial se concentrează mai mult pe implementarea și demonstrarea circuitelor privind algoritmul lui Shor. Pentru o resursă educațională aprofundată despre material, consultă cursul Fundamentele algoritmilor cuantici al Dr. John Watrous și articolele din secțiunea Referințe.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele:

• Qiskit SDK v2.0 sau o versiune ulterioară, cu suport pentru vizualizare
• Qiskit Runtime v0.40 sau o versiune ulterioară (pip install qiskit-ibm-runtime)

Configurare

# 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

Pasul 1: Maparea intrărilor clasice la o problemă cuantică

Context

Algoritmul lui Shor pentru factorizarea numerelor întregi utilizează o problemă intermediară cunoscută sub numele de problema găsirii ordinului. În această secțiune, demonstrăm cum să rezolvăm problema găsirii ordinului folosind estimarea cuantică a fazei.

Problema estimării fazei

În problema estimării fazei, ni se dă o stare cuantică $\ket{\psi}$ cu $n$ qubiți, împreună cu un Circuit cuantic unitar care acționează pe $n$ qubiți. Ni se garantează că $\ket{\psi}$ este un vector propriu al matricei unitare $U$ care descrie acțiunea circuitului, iar scopul nostru este să calculăm sau să aproximăm valoarea proprie $\lambda = e^{2 \pi i \theta}$ căreia îi corespunde $\ket{\psi}$. Cu alte cuvinte, circuitul ar trebui să producă o aproximare a numărului $\theta \in [0, 1)$ satisfăcând $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ Scopul circuitului de estimare a fazei este să aproximeze $\theta$ în $m$ biți. Matematic vorbind, am dori să găsim $y$ astfel încât $\theta \approx y / 2^m$, unde $y \in {0, 1, 2, \dots, 2^{m-1}}$. Imaginea de mai jos arată circuitul cuantic care estimează $y$ în $m$ biți prin efectuarea unei măsurători pe $m$ qubiți. În circuitul de mai sus, primii $m$ qubiți sunt inițializați în starea $\ket{0^m}$, iar ultimii $n$ qubiți sunt inițializați în $\ket{\psi}$, despre care se garantează că este un vector propriu al lui $U$. Primul ingredient din circuitul de estimare a fazei îl constituie operațiile unitare controlate, responsabile de efectuarea unui phase kickback la qubit-ul de control corespunzător. Aceste unitare controlate sunt exponențiate în conformitate cu poziția qubit-ului de control, de la bitul cel mai puțin semnificativ la bitul cel mai semnificativ. Deoarece $\ket{\psi}$ este un vector propriu al lui $U$, starea ultimilor $n$ qubiți nu este afectată de această operație, dar informația de fază a valorii proprii se propagă la primii $m$ qubiți. Se dovedește că după operația de phase kickback prin intermediul unitarelor controlate, toate stările posibile ale primilor $m$ qubiți sunt ortonormale între ele pentru fiecare vector propriu $\ket{\psi}$ al unitarului $U$. Prin urmare, aceste stări sunt perfect distinguibile și putem roti baza pe care o formează înapoi la baza computațională pentru a efectua o măsurătoare. O analiză matematică arată că această matrice de rotație corespunde transformatei Fourier cuantice inverse (QFT) în spațiul Hilbert de dimensiune $2^m$. Intuiția din spatele acestui lucru este că structura periodică a operatorilor de exponențiere modulară este codificată în starea cuantică, iar QFT convertește această periodicitate în vârfuri măsurabile în domeniul frecvențelor.

Pentru o înțelegere mai aprofundată a motivului pentru care circuitul QFT este utilizat în algoritmul lui Shor, îl îndrumăm pe cititor către cursul Fundamentele algoritmilor cuantici. Suntem acum pregătiți să folosim circuitul de estimare a fazei pentru găsirea ordinului.

Problema găsirii ordinului

Pentru a defini problema găsirii ordinului, începem cu câteva concepte de teoria numerelor. Mai întâi, pentru orice număr întreg pozitiv $N$, definim mulțimea $\mathbb{Z}_N$ ca $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Toate operațiile aritmetice în $\mathbb{Z}_N$ se efectuează modulo $N$. În particular, toate elementele $a \in \mathbb{Z}_n$ care sunt coprime cu $N$ sunt speciale și constituie $\mathbb{Z}^*_N$ ca $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Pentru un element $a \in \mathbb{Z}^*_N$, cel mai mic număr întreg pozitiv $r$ astfel încât $a^r \equiv 1 \; (\mathrm{mod} \; N)$ este definit ca ordinul lui $a$ modulo $N$. Așa cum vom vedea mai târziu, găsirea ordinului unui $a \in \mathbb{Z}^*_N$ ne va permite să factorizăm $N$. Pentru a construi circuitul de găsire a ordinului din circuitul de estimare a fazei, avem nevoie de două considerente. În primul rând, trebuie să definim unitarul $U$ care ne va permite să găsim ordinul $r$, și în al doilea rând, trebuie să definim un vector propriu $\ket{\psi}$ al lui $U$ pentru a pregăti starea inițială a circuitului de estimare a fazei.

Pentru a conecta problema găsirii ordinului la estimarea fazei, considerăm operația definită pe un sistem ale cărui stări clasice corespund lui $\mathbb{Z}_N$, unde înmulțim cu un element fix $a \in \mathbb{Z}^*_N$. În particular, definim acest operator de înmulțire $M_a$ astfel încât $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ pentru fiecare $x \in \mathbb{Z}_N$. Reține că este implicit că luăm produsul modulo $N$ în interiorul ket-ului din dreapta ecuației. O analiză matematică arată că $M_a$ este un operator unitar. Mai mult, se dovedește că $M_a$ are perechi de vectori proprii și valori proprii care ne permit să conectăm ordinul $r$ al lui $a$ la problema estimării fazei. Specific, pentru orice alegere a lui $j \in \{0, \dots, r-1\}$, avem că $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ este un vector propriu al lui $M_a$ a cărui valoare proprie corespunzătoare este $\omega^{j}_{r}$, unde $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Prin observație, vedem că o pereche comodă vector propriu/valoare proprie este starea $\ket{\psi_1}$ cu $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Prin urmare, dacă am putea găsi vectorul propriu $\ket{\psi_1}$, am putea estima faza $\theta=1/r$ cu circuitul nostru cuantic și, prin urmare, am obține o estimare a ordinului $r$. Cu toate acestea, nu este ușor de făcut, și trebuie să luăm în considerare o alternativă.

Să considerăm ce ar rezulta din circuit dacă pregătim starea computațională $\ket{1}$ ca stare inițială. Aceasta nu este o stare proprie a lui $M_a$, dar este superpozița uniformă a stărilor proprii pe care tocmai le-am descris. Cu alte cuvinte, următoarea relație este valabilă. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ Implicația ecuației de mai sus este că, dacă setăm starea inițială la $\ket{1}$, vom obține exact același rezultat al măsurătorii ca și cum am fi ales $k \in \{ 0, \dots, r-1\}$ uniform la întâmplare și am fi folosit $\ket{\psi_k}$ ca vector propriu în circuitul de estimare a fazei. Cu alte cuvinte, o măsurătoare a primilor $m$ qubiți produce o aproximare $y / 2^m$ a valorii $k / r$, unde $k \in \{ 0, \dots, r-1\}$ este ales uniform la întâmplare. Aceasta ne permite să aflăm $r$ cu un grad ridicat de încredere după mai multe rulări independente, ceea ce a fost scopul nostru.

Operatori de exponențiere modulară

Până acum, am legat problema estimării fazei de problema găsirii ordinului prin definirea $U = M_a$ și $\ket{\psi} = \ket{1}$ în circuitul nostru cuantic. Prin urmare, ultimul ingredient rămas este să găsim o modalitate eficientă de a defini exponențialele modulare ale lui $M_a$ ca $M_a^k$ pentru $k = 1, 2, 4, \dots, 2^{m-1}$. Pentru a efectua acest calcul, găsim că pentru orice putere $k$ pe care o alegem, putem crea un circuit pentru $M_a^k$ nu iterând de $k$ ori circuitul pentru $M_a$, ci în schimb calculând $b = a^k \; \mathrm{mod} \; N$ și folosind apoi circuitul pentru $M_b$. Deoarece avem nevoie doar de puterile care sunt ele însele puteri ale lui 2, putem face acest lucru eficient în mod clasic folosind ridicarea iterativă la pătrat.

Pasul 2: Optimizarea problemei pentru execuția pe hardware cuantic

Exemplu specific cu $N = 15$ și $a=2$

Putem face o pauză aici pentru a discuta un exemplu specific și a construi circuitul de găsire a ordinului pentru $N=15$. Reține că posibilele valori netriviale $a \in \mathbb{Z}_N^*$ pentru $N=15$ sunt $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Pentru acest exemplu, alegem $a=2$. Vom construi operatorul $M_2$ și operatorii de exponențiere modulară $M_2^k$. Acțiunea lui $M_2$ asupra stărilor bazei computaționale este următoarea. $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}$ Prin observație, putem vedea că stările bazei sunt amestecate, deci avem o matrice de permutare. Putem construi această operație pe patru qubiți cu Gate-uri swap. Mai jos, construim operațiile $M_2$ și $M_2$ controlat.

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)

Gate-urile care acționează pe mai mult de doi qubiți vor fi în continuare descompuse în Gate-uri cu doi qubiți.

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

Acum trebuie să construim operatorii de exponențiere modulară. Pentru a obține suficientă precizie în estimarea fazei, vom folosi opt qubiți pentru măsurarea estimării. Prin urmare, trebuie să construim $M_b$ cu $b = a^{2^k} \; (\mathrm{mod} \; N)$ pentru fiecare $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]

După cum putem vedea din lista valorilor $b$, pe lângă $M_2$ pe care l-am construit anterior, trebuie să construim și $M_4$ și $M_1$. Reține că $M_1$ acționează trivial asupra stărilor bazei computaționale, deci este pur și simplu operatorul identitate.

$M_4$ acționează asupra stărilor bazei computaționale după cum urmează. $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}$

Prin urmare, această permutare poate fi construită cu următoarea operație 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)

Gate-urile care acționează pe mai mult de doi qubiți vor fi în continuare descompuse în Gate-uri cu doi qubiți.

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

Am văzut că operatorii $M_b$ pentru un $b \in \mathbb{Z}^*_N$ dat sunt operații de permutare. Datorită dimensiunii relativ mici a problemei de permutare pe care o avem aici, deoarece $N=15$ necesită doar patru qubiți, am putut sintetiza aceste operații direct cu Gate-uri SWAP prin inspecție. În general, aceasta ar putea să nu fie o abordare scalabilă. În schimb, ar putea fi necesar să construim matricea de permutare explicit și să folosim clasa UnitaryGate a Qiskit și metodele de transpilare pentru a sintetiza această matrice de permutare. Cu toate acestea, aceasta poate duce la circuite semnificativ mai adânci. Urmează un exemplu.

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

Să comparăm aceste cifre cu adâncimea circuitului compilat a implementării noastre manuale a Gate-ului $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})

După cum putem vedea, abordarea bazată pe matricea de permutare a rezultat într-un circuit semnificativ mai adânc chiar și pentru un singur Gate $M_2$ față de implementarea noastră manuală. Prin urmare, vom continua cu implementarea anterioară a operațiilor $M_b$. Acum suntem gata să construim circuitul complet de găsire a ordinului folosind operatorii de exponențiere modulară controlați definiți anterior. În codul următor, importăm și circuitul QFT din biblioteca de circuite Qiskit, care utilizează Gate-uri Hadamard pe fiecare qubit, o serie de Gate-uri controlate-U1 (sau Z, în funcție de fază) și un strat de Gate-uri 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)

Reține că am omis operațiile de exponențiere modulară controlate de la qubiții de control rămași, deoarece $M_1$ este operatorul identitate. Reține că mai târziu în acest tutorial vom rula acest circuit pe Backend-ul ibm_marrakesh. Pentru a face acest lucru, transpilăm circuitul conform acestui Backend specific și raportăm adâncimea circuitului și numărul de Gate-uri.

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

Pasul 3: Execuție folosind primitivele Qiskit

Mai întâi, discutăm ce am obține teoretic dacă am rula acest circuit pe un simulator ideal. Mai jos avem un set de rezultate de simulare ale circuitului de mai sus folosind 1024 de măsurători. După cum se poate observa, obținem o distribuție aproximativ uniformă peste patru șiruri de biți pe qubiții de control.

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

plot_histogram(counts)

Prin măsurarea qubiților de control, obținem o estimare de fază pe opt biți a operatorului $M_a$. Putem converti această reprezentare binară în zecimal pentru a găsi faza măsurată. După cum se poate vedea din histograma de mai sus, au fost măsurate patru șiruri de biți diferite, fiecare corespunzând unei valori de fază după cum urmează.

# 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

Reamintim că orice fază măsurată corespunde la $\theta = k / r$, unde $k$ este eșantionat uniform aleator din $\{0, 1, \dots, r-1 \}$. Prin urmare, putem folosi algoritmul fracțiilor continue pentru a încerca să găsim $k$ și ordinul $r$. Python are această funcționalitate încorporată. Putem folosi modulul fractions pentru a transforma un număr real într-un obiect Fraction, de exemplu:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Deoarece aceasta returnează fracții exacte (în acest caz, 0.6660000...), pot apărea rezultate neplăcute precum cel de mai sus. Putem folosi metoda .limit_denominator() pentru a obține fracția care aproximează cel mai bine numărul nostru real, cu un numitor sub o anumită valoare:

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

Fraction(2, 3)

Mult mai simplu. Ordinul (r) trebuie să fie mai mic decât N, deci vom seta numitorul maxim la 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

Putem vedea că două dintre valorile proprii măsurate ne-au furnizat rezultatul corect: $r=4$, și putem constata că algoritmul lui Shor pentru găsirea ordinului poate eșua. Aceste rezultate proaste apar deoarece $k = 0$, sau deoarece $k$ și $r$ nu sunt coprime — în loc de $r$, obținem un factor al lui $r$. Cea mai simplă soluție este să repetăm pur și simplu experimentul până când obținem un rezultat satisfăcător pentru $r$.

Până acum, am implementat problema găsirii ordinului pentru $N=15$ cu $a=2$ folosind circuitul de estimare a fazei pe un simulator. Ultimul pas al algoritmului lui Shor va fi să legăm problema găsirii ordinului de problema factorizării numerelor întregi. Această ultimă parte a algoritmului este pur clasică și poate fi rezolvată pe un calculator clasic după ce măsurătorile de fază au fost obținute de la un calculator cuantic. Prin urmare, amânăm ultima parte a algoritmului până după ce demonstrăm cum putem rula circuitul de găsire a ordinului pe hardware real.

Rulări pe hardware

Acum putem rula circuitul de găsire a ordinului pe care l-am transpilat anterior pentru ibm_marrakesh. Apelăm la decuplare dinamică (DD) pentru suprimarea erorilor, și la răsucire de porți pentru atenuarea erorilor. Decuplarea dinamică implică aplicarea unor secvențe de impulsuri de control precise în timp pe un dispozitiv cuantic, eliminând eficient interacțiunile nedorite cu mediul și decoerența. Răsucirea de porți, pe de altă parte, randomizează porți cuantice specifice pentru a transforma erorile coerente în erori Pauli, care se acumulează liniar, nu pătratic. Ambele tehnici sunt combinate adesea pentru a îmbunătăți coerența și fidelitatea calculelor cuantice.

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

După cum se poate vedea, am obținut aceleași șiruri de biți cu cele mai mari numărători. Deoarece hardware-ul cuantic are zgomot, există o scurgere către alte șiruri de biți, pe care le putem filtra statistic.

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

Pasul 4: Post-procesare și returnarea rezultatului în formatul clasic dorit

Factorizarea numerelor întregi

Până acum, am discutat cum putem implementa problema găsirii ordinului folosind un circuit de estimare a fazei. Acum, conectăm problema găsirii ordinului la factorizarea numerelor întregi, ceea ce completează algoritmul lui Shor. Reține că această parte a algoritmului este clasică.

Demonstrăm acum acest lucru folosind exemplul nostru cu $N = 15$ și $a = 2$. Reamintim că faza măsurată este $k / r$, unde $a^r \; (\textrm{mod} \; N) = 1$ și $k$ este un număr întreg aleator între $0$ și $r - 1$. Din această ecuație, avem $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ ceea ce înseamnă că $N$ trebuie să dividă $a^r-1$. Dacă $r$ este și par, atunci putem scrie $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Dacă $r$ nu este par, nu putem merge mai departe și trebuie să încercăm din nou cu o altă valoare pentru $a$; altfel, există o probabilitate mare ca cel mai mare divizor comun al lui $N$ și fie $a^{r/2}-1$, fie $a^{r/2}+1$ să fie un factor propriu al lui $N$.

Deoarece unele rulări ale algoritmului vor eșua statistic, vom repeta acest algoritm până când cel puțin un factor al lui $N$ este găsit.

Celula de mai jos repetă algoritmul până când cel puțin un factor al lui $N=15$ este găsit. Vom folosi rezultatele rulării hardware de mai sus pentru a estima faza și factorul corespunzător în fiecare iterație.

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

Discuție

Lucrări conexe

În această secțiune, discutăm alte lucrări de referință care au demonstrat algoritmul lui Shor pe hardware real.

Lucrarea seminală [3] de la IBM® a demonstrat pentru prima dată algoritmul lui Shor, factorizând numărul 15 în factorii săi primi 3 și 5 folosind un calculator cuantic cu rezonanță magnetică nucleară (RMN) cu șapte qubiți. Un alt experiment [4] a factorizat 15 folosind qubiți fotonici. Prin utilizarea unui singur qubit reciclat de mai multe ori și codificând registrul de lucru în stări de dimensiune mai mare, cercetătorii au redus numărul de qubiți necesari la o treime față de protocolul standard, utilizând un algoritm compilat cu doi fotoni. O lucrare importantă în demonstrarea algoritmului lui Shor este [5], care folosește tehnica de estimare iterativă a fazei a lui Kitaev [8] pentru a reduce necesarul de qubiți al algoritmului. Autorii au folosit șapte qubiți de control și patru qubiți cache, împreună cu implementarea multiplicatorilor modulari. Această implementare necesită însă măsurători de mijloc de circuit cu operații de propagare înainte și reciclarea qubiților cu operații de resetare. Această demonstrație a fost realizată pe un calculator cuantic cu capcană ionică.

Lucrări mai recente [6] s-au concentrat pe factorizarea numerelor 15, 21 și 35 pe hardware-ul IBM Quantum®. Similar lucrărilor anterioare, cercetătorii au folosit o versiune compilată a algoritmului care a utilizat o transformată cuantică Fourier semiclasică, propusă de Kitaev, pentru a minimiza numărul de qubiți fizici și porți. O lucrare foarte recentă [7] a realizat și ea o demonstrație de tip dovadă de concept pentru factorizarea întregului 21. Această demonstrație a implicat de asemenea utilizarea unei versiuni compilate a rutinei de estimare cuantică a fazei și a extins demonstrația anterioară din [4]. Autorii au depășit această lucrare prin utilizarea unei configurații de porți Toffoli aproximative cu defazaje reziduale. Algoritmul a fost implementat pe procesoare cuantice IBM folosind doar cinci qubiți, iar prezența entanglementului între qubiții de control și cei de registru a fost verificată cu succes.

Scalabilitatea algoritmului

Reținem că criptografia RSA implică de obicei dimensiuni de cheie de ordinul 2048 până la 4096 de biți. Încercarea de a factoriza un număr de 2048 de biți cu algoritmul lui Shor va produce un circuit cuantic cu milioane de qubiți, inclusiv costul corecției erorilor, și o adâncime a circuitului de ordinul unui miliard, ceea ce depășește limitele hardware-ului cuantic actual. Prin urmare, algoritmul lui Shor va necesita fie metode optimizate de construire a circuitelor, fie o corecție robustă a erorilor cuantice pentru a fi practic viabil în spargerea sistemelor criptografice moderne. Te îndreptăm către [9] pentru o discuție mai detaliată despre estimarea resurselor necesare algoritmului lui Shor.

Provocare

Felicitări pentru finalizarea tutorialului! Acum este un moment excelent pentru a-ți testa înțelegerea. Poți încerca să construiești circuitul pentru factorizarea lui 21? Poți alege un $a$ la propria ta alegere. Va trebui să decizi asupra preciziei în biți a algoritmului pentru a alege numărul de qubiți, și va trebui să proiectezi operatorii de exponențiere modulară $M_a$. Te încurajăm să încerci singur, și apoi să citești despre metodologiile prezentate în Fig. 9 din [6] și Fig. 2 din [7].

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

    # Your code here
    pass

Referințe

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.