Algorithme de Shor

Estimation d'utilisation : trois secondes sur un processeur Eagle r3 (REMARQUE : il s'agit uniquement d'une estimation. Votre temps d'exécution peut varier.)

Résultats d'apprentissage

Après avoir suivi ce tutoriel, les utilisateurs devraient comprendre :

• Le contexte mathématique de l'algorithme de Shor pour la factorisation d'entiers
• Comment exécuter un exemple de cet algorithme sur du matériel réel

Prérequis

Nous suggérons que les utilisateurs soient familiers avec les sujets suivants avant de suivre ce tutoriel :

• Fondamentaux des algorithmes quantiques.
• Estimation de phase et factorisation. Nous couvrons une partie de ce contenu dans ce tutoriel.

Contexte

L'algorithme de Shor, développé par Peter Shor en 1994, est un algorithme quantique révolutionnaire pour la factorisation d'entiers en temps polynomial. Son importance réside dans sa capacité à factoriser de grands entiers exponentiellement plus rapidement que tout algorithme classique connu, menaçant ainsi la sécurité des systèmes cryptographiques largement utilisés comme RSA, qui reposent sur la difficulté de factoriser de grands nombres. En résolvant efficacement ce problème sur un ordinateur quantique suffisamment puissant, l'algorithme de Shor pourrait révolutionner des domaines tels que la cryptographie, la cybersécurité et les mathématiques computationnelles, soulignant la puissance transformatrice du calcul quantique.

Ce tutoriel se concentre sur la démonstration de l'algorithme de Shor en factorisant 15 sur un ordinateur quantique.

Tout d'abord, nous définissons le problème de recherche d'ordre et construisons les circuits correspondants à partir du protocole d'estimation de phase quantique. Ensuite, nous exécutons les circuits de recherche d'ordre sur du matériel réel en utilisant les circuits de profondeur minimale que nous pouvons transpiler. La dernière section complète l'algorithme de Shor en reliant le problème de recherche d'ordre à la factorisation d'entiers.

Nous terminons le tutoriel par une discussion sur d'autres démonstrations de l'algorithme de Shor sur du matériel réel, en nous concentrant à la fois sur les implémentations génériques et celles adaptées à la factorisation d'entiers spécifiques tels que 15 et 21. Remarque : ce tutoriel se concentre davantage sur l'implémentation et la démonstration des circuits liés à l'algorithme de Shor. Pour une ressource pédagogique approfondie sur le sujet, consulte le cours Fundamentals of quantum algorithms du Dr. John Watrous, ainsi que les articles dans la section Références.

Prérequis

Avant de commencer ce tutoriel, assure-toi que les éléments suivants sont installés :

• Qiskit SDK v2.0 ou ultérieur, avec le support de visualisation
• Qiskit Runtime v0.40 ou ultérieur (pip install qiskit-ibm-runtime)

Configuration

# 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

Étape 1 : Traduire les entrées classiques en un problème quantique

L'algorithme de Shor pour la factorisation d'entiers utilise un problème intermédiaire connu sous le nom de problème de recherche d'ordre. Dans cette section, nous démontrons comment résoudre le problème de recherche d'ordre en utilisant l'estimation de phase quantique.

Problème d'estimation de phase

Dans le problème d'estimation de phase, on nous donne un état quantique $\ket{\psi}$ de $n$ qubits, ainsi qu'un circuit quantique unitaire agissant sur $n$ qubits. On nous garantit que $\ket{\psi}$ est un vecteur propre de la matrice unitaire $U$ qui décrit l'action du circuit, et notre objectif est de calculer ou d'approximer la valeur propre $\lambda = e^{2 \pi i \theta}$ à laquelle $\ket{\psi}$ correspond. En d'autres termes, le circuit doit produire une approximation du nombre $\theta \in [0, 1)$ satisfaisant $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ L'objectif du circuit d'estimation de phase est d'approximer $\theta$ sur $m$ bits. Mathématiquement parlant, nous souhaitons trouver $y$ tel que $\theta \approx y / 2^m$, où $y \in {0, 1, 2, \dots, 2^{m-1}}$. L'image suivante montre le circuit quantique qui estime $y$ sur $m$ bits en effectuant une mesure sur $m$ qubits. Dans le circuit ci-dessus, les $m$ qubits du haut sont initialisés dans l'état $\ket{0^m}$, et les $n$ qubits du bas sont initialisés dans $\ket{\psi}$, qui est garanti être un vecteur propre de $U$. Le premier ingrédient du circuit d'estimation de phase sont les opérations unitaires contrôlées qui sont responsables de l'exécution d'un retour de phase (phase kickback) vers leur qubit de contrôle correspondant. Ces unitaires contrôlés sont exponentiés en fonction de la position du qubit de contrôle, allant du bit le moins significatif au bit le plus significatif. Puisque $\ket{\psi}$ est un vecteur propre de $U$, l'état des $n$ qubits du bas n'est pas affecté par cette opération, mais l'information de phase de la valeur propre se propage vers les $m$ qubits du haut. Il s'avère qu'après l'opération de retour de phase via les unitaires contrôlés, tous les états possibles des $m$ qubits du haut sont orthonormaux les uns par rapport aux autres pour chaque vecteur propre $\ket{\psi}$ de l'unitaire $U$. Par conséquent, ces états sont parfaitement distinguables, et nous pouvons effectuer une rotation de la base qu'ils forment vers la base computationnelle pour effectuer une mesure. Une analyse mathématique montre que cette matrice de rotation correspond à la transformée de Fourier quantique (QFT) inverse dans un espace de Hilbert de dimension $2^m$. L'intuition derrière cela est que la structure périodique des opérateurs d'exponentiation modulaire est encodée dans l'état quantique, et la QFT convertit cette périodicité en pics mesurables dans le domaine fréquentiel.

Pour une compréhension plus approfondie de la raison pour laquelle le circuit QFT est employé dans l'algorithme de Shor, nous renvoyons le lecteur au cours Fundamentals of quantum algorithms. Nous sommes maintenant prêts à utiliser le circuit d'estimation de phase pour la recherche d'ordre.

Problème de recherche d'ordre

Pour définir le problème de recherche d'ordre, nous commençons par quelques concepts de théorie des nombres. Premièrement, pour tout entier positif $N$ donné, on définit l'ensemble $\mathbb{Z}_N$ comme $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Toutes les opérations arithmétiques dans $\mathbb{Z}_N$ sont effectuées modulo $N$. En particulier, tous les éléments $a \in \mathbb{Z}_n$ qui sont premiers avec $N$ sont spéciaux et constituent $\mathbb{Z}^*_N$ tel que $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Pour un élément $a \in \mathbb{Z}^*_N$, le plus petit entier positif $r$ tel que $a^r \equiv 1 \; (\mathrm{mod} \; N)$ est défini comme l'ordre de $a$ modulo $N$. Comme nous le verrons plus tard, trouver l'ordre d'un $a \in \mathbb{Z}^*_N$ nous permettra de factoriser $N$. Pour construire le circuit de recherche d'ordre à partir du circuit d'estimation de phase, nous avons besoin de deux considérations. Premièrement, nous devons définir l'unitaire $U$ qui nous permettra de trouver l'ordre $r$, et deuxièmement, nous devons définir un vecteur propre $\ket{\psi}$ de $U$ pour préparer l'état initial du circuit d'estimation de phase.

Pour relier le problème de recherche d'ordre à l'estimation de phase, nous considérons l'opération définie sur un système dont les états classiques correspondent à $\mathbb{Z}_N$, où nous multiplions par un élément fixe $a \in \mathbb{Z}^*_N$. En particulier, nous définissons cet opérateur de multiplication $M_a$ tel que $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ pour chaque $x \in \mathbb{Z}_N$. Note qu'il est implicite que nous prenons le produit modulo $N$ à l'intérieur du ket du côté droit de l'équation. Une analyse mathématique montre que $M_a$ est un opérateur unitaire. De plus, il s'avère que $M_a$ possède des paires vecteur propre/valeur propre qui nous permettent de relier l'ordre $r$ de $a$ au problème d'estimation de phase. Plus précisément, pour tout choix de $j \in \{0, \dots, r-1\}$, nous avons que $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ est un vecteur propre de $M_a$ dont la valeur propre correspondante est $\omega^{j}_{r}$, où $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Par observation, nous voyons qu'une paire vecteur propre/valeur propre pratique est l'état $\ket{\psi_1}$ avec $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Par conséquent, si nous pouvions trouver le vecteur propre $\ket{\psi_1}$, nous pourrions estimer la phase $\theta=1/r$ avec notre circuit quantique et ainsi obtenir une estimation de l'ordre $r$. Cependant, ce n'est pas facile à faire, et nous devons envisager une alternative.

Considérons ce que le circuit donnerait si nous préparions l'état computationnel $\ket{1}$ comme état initial. Ce n'est pas un état propre de $M_a$, mais c'est la superposition uniforme des états propres que nous venons de décrire ci-dessus. En d'autres termes, la relation suivante est vérifiée. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ L'implication de l'équation ci-dessus est que si nous fixons l'état initial à $\ket{1}$, nous obtiendrons précisément le même résultat de mesure que si nous avions choisi $k \in \{ 0, \dots, r-1\}$ uniformément au hasard et utilisé $\ket{\psi_k}$ comme vecteur propre dans le circuit d'estimation de phase. En d'autres termes, une mesure des $m$ qubits du haut donne une approximation $y / 2^m$ de la valeur $k / r$ où $k \in \{ 0, \dots, r-1\}$ est choisi uniformément au hasard. Cela nous permet d'apprendre $r$ avec un haut degré de confiance après plusieurs exécutions indépendantes, ce qui était notre objectif.

Opérateurs d'exponentiation modulaire

Jusqu'à présent, nous avons relié le problème d'estimation de phase au problème de recherche d'ordre en définissant $U = M_a$ et $\ket{\psi} = \ket{1}$ dans notre circuit quantique. Par conséquent, le dernier ingrédient restant est de trouver un moyen efficace de définir les exponentielles modulaires de $M_a$ sous la forme $M_a^k$ pour $k = 1, 2, 4, \dots, 2^{m-1}$. Pour effectuer ce calcul, nous constatons que pour toute puissance $k$ choisie, nous pouvons créer un circuit pour $M_a^k$ non pas en itérant $k$ fois le circuit pour $M_a$, mais plutôt en calculant $b = a^k \; \mathrm{mod} \; N$ puis en utilisant le circuit pour $M_b$. Puisque nous n'avons besoin que des puissances qui sont elles-mêmes des puissances de 2, nous pouvons le faire classiquement de manière efficace en utilisant l'élévation au carré itérative.

Étape 2 : Optimiser le problème pour l'exécution sur du matériel quantique

Exemple spécifique avec $N = 15$ et $a=2$

Nous pouvons faire une pause ici pour discuter d'un exemple spécifique et construire le circuit de recherche d'ordre pour $N=15$. Note que les valeurs non triviales possibles $a \in \mathbb{Z}_N^*$ pour $N=15$ sont $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Pour cet exemple, nous choisissons $a=2$. Nous allons construire l'opérateur $M_2$ et les opérateurs d'exponentiation modulaire $M_2^k$. L'action de $M_2$ sur les états de la base computationnelle est la suivante. $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}$ Par observation, nous pouvons voir que les états de base sont permutés, nous avons donc une matrice de permutation. Nous pouvons construire cette opération sur quatre qubits avec des portes swap. Ci-dessous, nous construisons les opérations $M_2$ et $M_2$ contrôlée.

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)

Les portes agissant sur plus de deux qubits seront davantage décomposées en portes à deux qubits.

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

Nous devons maintenant construire les opérateurs d'exponentiation modulaire. Pour obtenir une précision suffisante dans l'estimation de phase, nous utiliserons huit qubits pour la mesure d'estimation. Par conséquent, nous devons construire $M_b$ avec $b = a^{2^k} \; (\mathrm{mod} \; N)$ pour chaque $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]

Comme nous pouvons le voir dans la liste des valeurs de $b$, en plus de $M_2$ que nous avons précédemment construit, nous devons également construire $M_4$ et $M_1$. Note que $M_1$ agit de manière triviale sur les états de la base computationnelle, c'est donc simplement l'opérateur identité.

$M_4$ agit sur les états de la base computationnelle comme suit. $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}$

Par conséquent, cette permutation peut être construite avec l'opération swap suivante.

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)

Les portes agissant sur plus de deux qubits seront davantage décomposées en portes à deux qubits.

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

Nous avons vu que les opérateurs $M_b$ pour un $b \in \mathbb{Z}^*_N$ donné sont des opérations de permutation. En raison de la taille relativement petite du problème de permutation que nous avons ici, puisque $N=15$ ne nécessite que quatre qubits, nous avons pu synthétiser ces opérations directement avec des portes SWAP par inspection. En général, cette approche pourrait ne pas être extensible. Au lieu de cela, nous pourrions avoir besoin de construire explicitement la matrice de permutation et d'utiliser la classe UnitaryGate de Qiskit ainsi que les méthodes de transpilation pour synthétiser cette matrice de permutation. Cependant, cela peut entraîner des circuits significativement plus profonds. Un exemple suit.

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

Comparons ces chiffres avec la profondeur du circuit compilé de notre implémentation manuelle de la porte $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})

Comme nous pouvons le constater, l'approche par matrice de permutation a produit un circuit significativement plus profond même pour une seule porte $M_2$ par rapport à notre implémentation manuelle. Par conséquent, nous continuerons avec notre implémentation précédente des opérations $M_b$. Nous sommes maintenant prêts à construire le circuit complet de recherche d'ordre en utilisant nos opérateurs d'exponentiation modulaire contrôlés précédemment définis. Dans le code suivant, nous importons également le circuit QFT de la bibliothèque de circuits Qiskit, qui utilise des portes de Hadamard sur chaque qubit, une série de portes U1 contrôlées (ou Z, selon la phase) et une couche de portes 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)

Note que nous avons omis les opérations d'exponentiation modulaire contrôlées des qubits de contrôle restants car $M_1$ est l'opérateur identité. Note que plus loin dans ce tutoriel, nous exécuterons ce circuit sur le backend ibm_marrakesh. Pour ce faire, nous transpilons le circuit selon ce backend spécifique et rapportons la profondeur du circuit et le nombre de portes.

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

Étape 3 : Exécuter en utilisant les primitives Qiskit

Tout d'abord, nous discutons de ce que nous obtiendrions théoriquement si nous exécutions ce circuit sur un simulateur idéal. Ci-dessous, nous avons un ensemble de résultats de simulation du circuit ci-dessus utilisant 1024 tirs. Comme nous pouvons le voir, nous obtenons une distribution approximativement uniforme sur quatre chaînes de bits sur les qubits de contrôle.

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

plot_histogram(counts)

En mesurant les qubits de contrôle, nous obtenons une estimation de phase sur huit bits de l'opérateur $M_a$. Nous pouvons convertir cette représentation binaire en décimal pour trouver la phase mesurée. Comme nous pouvons le voir dans l'histogramme ci-dessus, quatre chaînes de bits différentes ont été mesurées, et chacune d'elles correspond à une valeur de phase comme suit.

# 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

Rappelons que toute phase mesurée correspond à $\theta = k / r$ où $k$ est échantillonné uniformément au hasard dans $\{0, 1, \dots, r-1 \}$. Par conséquent, nous pouvons utiliser l'algorithme des fractions continues pour tenter de trouver $k$ et l'ordre $r$. Python dispose de cette fonctionnalité intégrée. Nous pouvons utiliser le module fractions pour transformer un nombre à virgule flottante en objet Fraction, par exemple :

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Comme cela donne des fractions qui retournent le résultat exact (dans ce cas, 0.6660000...), cela peut produire des résultats peu élégants comme celui ci-dessus. Nous pouvons utiliser la méthode .limit_denominator() pour obtenir la fraction qui ressemble le plus à notre nombre à virgule flottante, avec un dénominateur inférieur à une certaine valeur :

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

Fraction(2, 3)

C'est beaucoup mieux. L'ordre (r) doit être inférieur à N, nous fixerons donc le dénominateur maximum à 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

Nous pouvons voir que deux des valeurs propres mesurées nous ont fourni le résultat correct : $r=4$, et nous pouvons constater que l'algorithme de Shor pour la recherche d'ordre a une chance d'échouer. Ces mauvais résultats sont dus au fait que $k = 0$, ou parce que $k$ et $r$ ne sont pas premiers entre eux -- et au lieu de $r$, nous obtenons un facteur de $r$. La solution la plus simple est de simplement répéter l'expérience jusqu'à obtenir un résultat satisfaisant pour $r$. Jusqu'à présent, nous avons implémenté le problème de recherche d'ordre pour $N=15$ avec $a=2$ en utilisant le circuit d'estimation de phase sur un simulateur. La dernière étape de l'algorithme de Shor sera de relier le problème de recherche d'ordre au problème de factorisation d'entiers. Cette dernière partie de l'algorithme est purement classique et peut être résolue sur un ordinateur classique une fois les mesures de phase obtenues à partir d'un ordinateur quantique. Par conséquent, nous reportons la dernière partie de l'algorithme jusqu'à ce que nous ayons démontré comment exécuter le circuit de recherche d'ordre sur du matériel réel.

Exécutions sur le matériel

Nous pouvons maintenant exécuter le circuit de recherche d'ordre que nous avons précédemment transpilé pour ibm_marrakesh. Ici, nous nous tournons vers le découplage dynamique (DD) pour la suppression d'erreurs, et le twirling de portes pour l'atténuation d'erreurs. Le DD consiste à appliquer des séquences d'impulsions de contrôle précisément temporisées à un dispositif quantique, moyennant efficacement les interactions environnementales indésirables et la décohérence. Le twirling de portes, quant à lui, randomise des portes quantiques spécifiques pour transformer les erreurs cohérentes en erreurs de Pauli, qui s'accumulent linéairement plutôt que quadratiquement. Les deux techniques sont souvent combinées pour améliorer la cohérence et la fidélité des calculs quantiques.

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

Comme nous pouvons le voir, nous avons obtenu les mêmes chaînes de bits avec les comptes les plus élevés. Étant donné que le matériel quantique est bruité, il y a une certaine fuite vers d'autres chaînes de bits, que nous pouvons filtrer statistiquement.

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

Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité

Factorisation d'entiers

Jusqu'à présent, nous avons discuté de la manière dont nous pouvons implémenter le problème de recherche d'ordre en utilisant un circuit d'estimation de phase. Maintenant, nous relions le problème de recherche d'ordre à la factorisation d'entiers, ce qui complète l'algorithme de Shor. Note que cette partie de l'algorithme est classique. Nous démontrons maintenant ceci en utilisant notre exemple de $N = 15$ et $a = 2$. Rappelons que la phase mesurée est $k / r$, où $a^r \; (\textrm{mod} \; N) = 1$ et $k$ est un entier aléatoire entre $0$ et $r - 1$. De cette équation, nous avons $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ ce qui signifie que $N$ doit diviser $a^r-1$. Si $r$ est également pair, alors nous pouvons écrire $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Si $r$ n'est pas pair, nous ne pouvons pas aller plus loin et devons réessayer avec une valeur différente de $a$ ; sinon, il y a une forte probabilité que le plus grand commun diviseur de $N$ et soit $a^{r/2}-1$, soit $a^{r/2}+1$ soit un facteur propre de $N$.

Puisque certaines exécutions de l'algorithme échoueront statistiquement, nous répéterons cet algorithme jusqu'à ce qu'au moins un facteur de $N$ soit trouvé. La cellule ci-dessous répète l'algorithme jusqu'à ce qu'au moins un facteur de $N=15$ soit trouvé. Nous utiliserons les résultats de l'exécution matérielle ci-dessus pour deviner la phase et le facteur correspondant à chaque itération.

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

Discussion

Travaux connexes

Dans cette section, nous discutons d'autres travaux marquants ayant démontré l'algorithme de Shor sur du matériel réel.

Le travail fondateur [3] d'IBM® a démontré l'algorithme de Shor pour la première fois, en factorisant le nombre 15 en ses facteurs premiers 3 et 5 à l'aide d'un ordinateur quantique à résonance magnétique nucléaire (RMN) de sept qubits. Une autre expérience [4] a factorisé 15 en utilisant des qubits photoniques. En employant un seul qubit recyclé plusieurs fois et en encodant le registre de travail dans des états de dimension supérieure, les chercheurs ont réduit le nombre de qubits requis à un tiers de celui du protocole standard, en utilisant un algorithme compilé à deux photons. Un article significatif dans la démonstration de l'algorithme de Shor est [5], qui utilise la technique d'estimation de phase itérative de Kitaev [8] pour réduire les besoins en qubits de l'algorithme. Les auteurs ont utilisé sept qubits de contrôle et quatre qubits de cache, ainsi que l'implémentation de multiplicateurs modulaires. Cette implémentation nécessite cependant des mesures en milieu de circuit avec des opérations de propagation conditionnelle et le recyclage de qubits avec des opérations de réinitialisation. Cette démonstration a été réalisée sur un ordinateur quantique à piège à ions.

Des travaux plus récents [6] se sont concentrés sur la factorisation de 15, 21 et 35 sur du matériel IBM Quantum®. De manière similaire aux travaux précédents, les chercheurs ont utilisé une version compilée de l'algorithme employant une transformée de Fourier quantique semi-classique telle que proposée par Kitaev pour minimiser le nombre de qubits physiques et de portes. Un travail des plus récents [7] a également réalisé une démonstration de preuve de concept pour la factorisation de l'entier 21. Cette démonstration a également impliqué l'utilisation d'une version compilée de la routine d'estimation de phase quantique, et s'est appuyée sur la démonstration précédente de [4]. Les auteurs sont allés au-delà de ce travail en utilisant une configuration de portes de Toffoli approximatives avec des déphasages résiduels. L'algorithme a été implémenté sur des processeurs quantiques IBM en utilisant seulement cinq qubits, et la présence d'intrication entre les qubits de contrôle et les qubits de registre a été vérifiée avec succès.

Mise à l'échelle de l'algorithme

Nous notons que le chiffrement RSA implique généralement des tailles de clés de l'ordre de 2048 à 4096 bits. Tenter de factoriser un nombre de 2048 bits avec l'algorithme de Shor produira un circuit quantique de millions de qubits, y compris la surcharge de correction d'erreurs, et une profondeur de circuit de l'ordre du milliard, ce qui dépasse les limites du matériel quantique actuel. Par conséquent, l'algorithme de Shor nécessitera soit des méthodes optimisées de construction de circuits, soit une correction d'erreurs quantique robuste pour être pratiquement viable pour le déchiffrement des systèmes cryptographiques modernes. Nous te renvoyons à [9] pour une discussion plus détaillée sur l'estimation des ressources pour l'algorithme de Shor.

Défi

Félicitations pour avoir terminé le tutoriel ! C'est le moment idéal pour tester ta compréhension. Essaie de construire le circuit pour factoriser 21. Tu peux sélectionner un $a$ de ton choix. Tu devras décider de la précision en bits de l'algorithme pour choisir le nombre de qubits, et tu devras concevoir les opérateurs d'exponentiation modulaire $M_a$. Nous t'encourageons à essayer par toi-même, puis à lire les méthodologies présentées dans la Fig. 9 de [6] et la Fig. 2 de [7].

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

    # Your code here
    pass

Références

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.