Алгоритм Шора

Оцінка використання: Три секунди на процесорі Eagle r3 (ПРИМІТКА: Це лише оцінка. Ваш час виконання може відрізнятися.)

Результати навчання

Після проходження цього посібника користувачі мають розуміти:

• Математичне підґрунтя алгоритму Шора для факторизації цілих чисел
• Як запустити приклад цього алгоритму на реальному обладнанні

Передумови

Ми рекомендуємо, щоб користувачі були знайомі з такими темами перед проходженням цього посібника:

• Основи квантових алгоритмів.
• Оцінка фази та факторизація. Ми розглядаємо частину цього матеріалу в цьому посібнику.

Передумови

Алгоритм Шора, розроблений Пітером Шором у 1994 році, є проривним квантовим алгоритмом для факторизації цілих чисел за поліноміальний час. Його значення полягає у здатності факторизувати великі цілі числа експоненційно швидше, ніж будь-який відомий класичний алгоритм, що загрожує безпеці широко використовуваних криптографічних систем, таких як RSA, які покладаються на складність факторизації великих чисел. Ефективно вирішуючи цю проблему на достатньо потужному квантовому комп'ютері, алгоритм Шора може революціонізувати такі галузі, як криптографія, кібербезпека та обчислювальна математика, підкреслюючи трансформаційну силу квантових обчислень.

Цей підручник зосереджений на демонстрації алгоритму Шора шляхом факторизації числа 15 на квантовому комп'ютері.

Спочатку ми визначаємо проблему знаходження порядку і будуємо відповідні схеми з протоколу квантової оцінки фази. Далі ми запускаємо схеми знаходження порядку на реальному обладnanні, використовуючи схеми найменшої глибини, які ми можемо транспілювати. Остання секція завершує алгоритм Шора, пов'язуючи проблему знаходження порядку з факторизацією цілих чисел.

Ми завершуємо підручник обговоренням інших демонстрацій алгоритму Шора на реальному обладнанні, зосереджуючись як на загальних реалізаціях, так і на тих, що адаптовані для факторизації конкретних цілих чисел, таких як 15 і 21. Примітка: Цей підручник більше зосереджений на реалізації та демонстрації схем, що стосуються алгоритму Шора. Для поглибленого навчального ресурсу з цього матеріалу, будь ласка, звернись до курсу Основи квантових алгоритмів доктора Джона Ватруса та статей у розділі Посилання.

Вимоги

Перш ніж розпочати цей підручник, переконайтеся, що у Вас встановлено наступне:

• Qiskit SDK v2.0 або пізніше, з підтримкою візуалізації
• Qiskit Runtime v0.40 або пізніше (pip install qiskit-ibm-runtime)

Налаштування

# 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

Крок 1: Відобразити класичні вхідні дані на квантову проблему

Алгоритм Шора для факторизації цілих чисел використовує проміжну проблему, відому як проблема знаходження порядку. У цьому розділі ми демонструємо, як розв'язати проблему знаходження порядку за допомогою квантової оцінки фази.

Проблема оцінки фази

У проблемі оцінки фази нам дано квантовий стан $\ket{\psi}$ з $n$ кубітів разом з унітарною квантовою схемою, яка діє на $n$ кубітів. Нам обіцяно, що $\ket{\psi}$ є власним вектором унітарної матриці $U$, яка описує дію схеми, і наша мета полягає в обчисленні або апроксимації власного значення $\lambda = e^{2 \pi i \theta}$, якому відповідає $\ket{\psi}$. Іншими словами, схема повинна вивести апроксимацію числа $\theta \in [0, 1)$, що задовольняє $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ Мета схеми оцінки фази полягає в апроксимації $\theta$ в $m$ бітах. Математично кажучи, ми хотіли б знайти $y$ таке, що $\theta \approx y / 2^m$, де $y \in {0, 1, 2, \dots, 2^{m-1}}$. Наступне зображення показує квантову схему, яка оцінює $y$ в $m$ бітах, виконуючи вимірювання на $m$ кубітах. На наведеній вище схемі верхні $m$ кубітів ініціалізовані в стані $\ket{0^m}$, а нижні $n$ кубітів ініціалізовані в $\ket{\psi}$, який, як обіцяно, є власним вектором $U$. Першим інгредієнтом у схемі оцінки фази є керовані унітарні операції, які відповідають за виконання зворотного впливу фази на відповідний керуючий кубіт. Ці керовані унітарні операції піднесені до степеня відповідно до позиції керуючого кубіта, починаючи від найменш значущого біта до найбільш значущого біта. Оскільки $\ket{\psi}$ є власним вектором $U$, стан нижніх $n$ кубітів не змінюється цією операцією, але інформація про фазу власного значення поширюється на верхні $m$ кубітів. Виявляється, що після операції зворотного впливу фази через керовані унітарні операції всі можливі стани верхніх $m$ кубітів є ортонормальними один до одного для кожного власного вектора $\ket{\psi}$ унітарного оператора $U$. Тому ці стани є повністю розрізнюваними, і ми можемо обернути базис, який вони формують, назад до обчислювального базису, щоб виконати вимірювання. Математичний аналіз показує, що ця матриця обертання відповідає оберненому квантовому перетворенню Фур'є (QFT) у $2^m$-вимірному гільбертовому просторі. Інтуїція полягає в тому, що періодична структура операторів модульного піднесення до степеня кодується в квантовому стані, і QFT перетворює цю періодичність у вимірювані піки в частотній області.

Для більш глибокого розуміння того, чому схема QFT використовується в алгоритмі Шора, ми відсилаємо читача до курсу Основи квантових алгоритмів. Тепер ми готові використовувати схему оцінки фази для знаходження порядку.

Проблема знаходження порядку

Щоб визначити проблему знаходження порядку, ми починаємо з деяких концепцій теорії чисел. По-перше, для будь-якого заданого додатного цілого числа $N$ визначимо множину $\mathbb{Z}_N$ як $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Усі арифметичні операції в $\mathbb{Z}_N$ виконуються за модулем $N$. Зокрема, всі елементи $a \in \mathbb{Z}_n$, які є взаємно простими з $N$, є особливими і утворюють $\mathbb{Z}^*_N$ як $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Для елемента $a \in \mathbb{Z}^*_N$ найменше додатне ціле число $r$ таке, що $a^r \equiv 1 \; (\mathrm{mod} \; N)$ визначається як порядок $a$ за модулем $N$. Як ми побачимо пізніше, знаходження порядку $a \in \mathbb{Z}^*_N$ дозволить нам факторизувати $N$. Щоб побудувати схему знаходження порядку зі схеми оцінки фази, нам потрібні дві умови. По-перше, нам потрібно визначити унітарний оператор $U$, який дозволить нам знайти порядок $r$, і по-друге, нам потрібно визначити власний вектор $\ket{\psi}$ оператора $U$, щоб підготувати початковий стан схеми оцінки фази.

Щоб пов'язати проблему знаходження порядку з оцінкою фази, ми розглядаємо операцію, визначену на системі, класичні стани якої відповідають $\mathbb{Z}_N$, де ми множимо на фіксований елемент $a \in \mathbb{Z}^*_N$. Зокрема, ми визначаємо цей оператор множення $M_a$ таким чином, що $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ для кожного $x \in \mathbb{Z}_N$. Зауважте, що неявно ми беремо добуток за модулем $N$ усередині кета з правої сторони рівняння. Математичний аналіз показує, що $M_a$ є унітарним оператором. Крім того, виявляється, що $M_a$ має пари власних векторів і власних значень, які дозволяють нам пов'язати порядок $r$ числа $a$ з проблемою оцінки фази. Зокрема, для будь-якого вибору $j \in \{0, \dots, r-1\}$ ми маємо, що $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ є власним вектором $M_a$, відповідне власне значення якого є $\omega^{j}_{r}$, де $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Спостерігаючи, ми бачимо, що зручною парою власний вектор/власне значення є стан $\ket{\psi_1}$ з $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Отже, якби ми могли знайти власний вектор $\ket{\psi_1}$, ми могли б оцінити фазу $\theta=1/r$ за допомогою нашої квантової схеми і, таким чином, отримати оцінку порядку $r$. Однак зробити це нелегко, і нам потрібно розглянути альтернативу.

Розглянемо, що б дала схема, якби ми підготували обчислювальний стан $\ket{1}$ як початковий стан. Це не є власним станом $M_a$, але це є рівномірною суперпозицією власних станів, які ми щойно описали вище. Іншими словами, виконується наступне співвідношення. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ Наслідком наведеного вище рівняння є те, що якщо ми встановимо початковий стан на $\ket{1}$, ми отримаємо точно такий самий результат вимірювання, як якби ми вибрали $k \in \{ 0, \dots, r-1\}$ рівномірно випадково і використали $\ket{\psi_k}$ як власний вектор у схемі оцінки фази. Іншими словами, вимірювання верхніх $m$ кубітів дає апроксимацію $y / 2^m$ до значення $k / r$, де $k \in \{ 0, \dots, r-1\}$ вибирається рівномірно випадково. Це дозволяє нам вивчити $r$ з високим ступенем впевненості після кількох незалежних запусків, що було нашою метою.

Оператори модульного піднесення до степеня

До цього часу ми пов'язали проблему оцінки фази з проблемою знаходження порядку, визначивши $U = M_a$ і $\ket{\psi} = \ket{1}$ у нашій квантовій схемі. Тому останнім інгредієнтом, що залишився, є знаходження ефективного способу визначення модульних експонент $M_a$ як $M_a^k$ для $k = 1, 2, 4, \dots, 2^{m-1}$. Щоб виконати це обчислення, ми виявляємо, що для будь-якого степеня $k$, який ми вибираємо, ми можемо створити схему для $M_a^k$ не шляхом $k$-разової ітерації схеми для $M_a$, а натомість обчислюючи $b = a^k \; \mathrm{mod} \; N$, а потім використовуючи схему для $M_b$. Оскільки нам потрібні лише степені, які самі є степенями 2, ми можемо зробити це класично ефективно, використовуючи ітеративне піднесення до квадрата.

Крок 2: Оптимізувати проблему для виконання на квантовому обладнанні

Конкретний приклад з $N = 15$ і $a=2$

Ми можемо зробити паузу тут, щоб обговорити конкретний приклад і побудувати схему знаходження порядку для $N=15$. Зауважте, що можливі нетривіальні $a \in \mathbb{Z}_N^*$ для $N=15$ є $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Для цього прикладу ми вибираємо $a=2$. Ми побудуємо оператор $M_2$ і оператори модульного піднесення до степеня $M_2^k$. Дія $M_2$ на базисні стани обчислювального базису є наступною. $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}$ Спостерігаючи, ми можемо бачити, що базисні стани перемішуються, тому ми маємо матрицю перестановки. Ми можемо побудувати цю операцію на чотирьох кубітах за допомогою гейтів обміну. Нижче ми конструюємо операції $M_2$ і керованого-$M_2$.

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)

Гейти, що діють на більше ніж два кубіти, будуть далі розкладені на двокубітні гейти.

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

Тепер нам потрібно побудувати оператори модульного піднесення до степеня. Щоб отримати достатню точність в оцінці фази, ми будемо використовувати вісім кубітів для вимірювання оцінки. Тому нам потрібно побудувати $M_b$ з $b = a^{2^k} \; (\mathrm{mod} \; N)$ для кожного $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]

Як ми бачимо зі списку значень $b$, на додаток до $M_2$, який ми раніше побудували, нам також потрібно побудувати $M_4$ і $M_1$. Зауважте, що $M_1$ діє тривіально на базисні стани обчислювального базису, тому це просто тотожний оператор.

$M_4$ діє на базисні стани обчислювального базису наступним чином. $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}$

Тому цю перестановку можна побудувати за допомогою наступної операції обміну.

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)

Гейти, що діють на більше ніж два кубіти, будуть далі розкладені на двокубітні гейти.

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

Ми побачили, що оператори $M_b$ для заданого $b \in \mathbb{Z}^*_N$ є операціями перестановки. Через відносно невеликий розмір проблеми перестановки, яку ми тут маємо, оскільки $N=15$ потребує лише чотири кубіти, ми змогли синтезувати ці операції безпосередньо за допомогою гейтів SWAP шляхом інспекції. Загалом це може бути не масштабованим підходом. Натомість нам може знадобитися побудувати матрицю перестановки явно і використовувати клас UnitaryGate Qiskit і методи транспіляції для синтезу цієї матриці перестановки. Однак це може призвести до значно глибших схем. Приклад наведено нижче.

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

Порівняємо ці значення з глибиною скомпільованої схеми нашої ручної реалізації гейта $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})

Як ми бачимо, підхід з матрицею перестановки призвів до значно глибокої схеми навіть для одного гейта $M_2$ порівняно з нашою ручною реалізацією. Тому ми продовжимо з нашою попередньою реалізацією операцій $M_b$. Тепер ми готові побудувати повну схему знаходження порядку, використовуючи наші раніше визначені керовані оператори модульного піднесення до степеня. У наступному коді ми також імпортуємо схему QFT з бібліотеки схем Qiskit, яка використовує гейти Адамара на кожному кубіті, серію керованих-U1 (або Z, залежно від фази) гейтів і шар гейтів обміну.

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

Зауважте, що ми опустили керовані операції модульного піднесення до степеня з решти керуючих кубітів, оскільки $M_1$ є тотожним оператором. Зауважте, що пізніше в цьому підручнику ми запустимо цю схему на бекенді ibm_marrakesh. Щоб зробити це, ми транспілюємо схему відповідно до цього конкретного бекенда і повідомляємо глибину схеми та кількість гейтів.

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

Крок 3: Виконання за допомогою примітивів Qiskit

Спочатку ми обговоримо, що ми теоретично отримали б, якби запустили цю схему на ідеальному симуляторі. Нижче наведено набір результатів симуляції вищезгаданої схеми з використанням 1024 вимірювань. Як ми бачимо, ми отримуємо приблизно рівномірний розподіл по чотирьох бітових рядках для контрольних кубітів.

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

plot_histogram(counts)

Вимірюючи контрольні кубіти, ми отримуємо восьмибітну оцінку фази оператора $M_a$. Ми можемо перетворити це двійкове представлення в десяткове, щоб знайти виміряну фазу. Як ми бачимо з гістограми вище, було виміряно чотири різні бітові рядки, і кожен з них відповідає значенню фази наступним чином.

# 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

Пригадайте, що будь-яка виміряна фаза відповідає $\theta = k / r$, де $k$ вибирається рівномірно випадково з $\{0, 1, \dots, r-1 \}$. Тому ми можемо використовувати алгоритм ланцюгових дробів, щоб спробувати знайти $k$ і порядок $r$. Python має цю функціональність вбудованою. Ми можемо використовувати модуль fractions, щоб перетворити число з плаваючою комою на об'єкт Fraction, наприклад:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Оскільки це дає дроби, які повертають результат точно (у цьому випадку 0.6660000...), це може давати складні результати, як наведений вище. Ми можемо використовувати метод .limit_denominator(), щоб отримати дріб, який найбільше нагадує наше число з плаваючою комою, із знаменником нижче певного значення:

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

Fraction(2, 3)

Це набагато краще. Порядок (r) повинен бути меншим за N, тому ми встановимо максимальний знаменник як 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

Ми бачимо, що два з виміряних власних значень дали нам правильний результат: $r=4$, і ми бачимо, що алгоритм Шора для знаходження порядку може давати збої. Ці погані результати виникають через те, що $k = 0$, або через те, що $k$ і $r$ не є взаємно простими - і замість $r$ ми отримуємо дільник $r$. Найпростіше рішення цієї проблеми - просто повторювати експеримент, поки ми не отримаємо задовільний результат для $r$. Поки що ми реалізували задачу знаходження порядку для $N=15$ з $a=2$, використовуючи схему оцінки фази на симуляторі. Останнім кроком алгоритму Шора буде зв'язати задачу знаходження порядку з задачею факторизації цілих чисел. Ця остання частина алгоритму є суто класичною і може бути розв'язана на класичному комп'ютері після того, як вимірювання фази були отримані з квантового комп'ютера. Тому ми відкладемо останню частину алгоритму до того моменту, коли продемонструємо, як можна запустити схему знаходження порядку на реальному обладнанні.

Запуски на обладнанні

Тепер ми можемо запустити схему знаходження порядку, яку ми раніше транспілювали для ibm_marrakesh. Тут ми звертаємось до динамічного роз'єднання (DD) для придушення помилок та перекручування вентилів для пом'якшення помилок. DD включає застосування послідовностей точно розрахованих за часом керуючих імпульсів до квантового пристрою, ефективно усереднюючи небажані взаємодії з середовищем та декогеренцію. Перекручування вентилів, з іншого боку, рандомізує специфічні квантові вентилі, щоб перетворити когерентні помилки в помилки Паулі, які накопичуються лінійно, а не квадратично. Обидва методи часто поєднуються для підвищення когерентності та точності квантових обчислень.

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

Як ми бачимо, ми отримали ті самі бітові рядки з найвищими значеннями підрахунків. Оскільки квантове обладнання має шум, існує деяка витоку до інших бітових рядків, яку ми можемо відфільтрувати статистично.

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

Крок 4: Постобробка та повернення результату в бажаному класичному форматі

Факторизація цілих чисел

Поки що ми обговорили, як ми можемо реалізувати задачу знаходження порядку, використовуючи схему оцінки фази. Тепер ми пов'язуємо задачу знаходження порядку з факторизацією цілих чисел, що завершує алгоритм Шора. Зауважте, що ця частина алгоритму є класичною. Тепер продемонструємо це, використовуючи наш приклад $N = 15$ і $a = 2$. Пригадайте, що виміряна фаза дорівнює $k / r$, де $a^r \; (\textrm{mod} \; N) = 1$ і $k$ - це випадкове ціле число між $0$ і $r - 1$. З цього рівняння ми маємо $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ що означає, що $N$ повинно ділити $a^r-1$. Якщо $r$ також є парним, то ми можемо записати $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Якщо $r$ не є парним, ми не можемо рухатись далі і повинні спробувати знову з іншим значенням для $a$; інакше існує висока ймовірність того, що найбільший спільний дільник $N$ і або $a^{r/2}-1$, або $a^{r/2}+1$ є власним дільником $N$.

Оскільки деякі запуски алгоритму статистично зазнають невдачі, ми повторюватимемо цей алгоритм, поки не буде знайдено принаймні один дільник $N$. Комірка нижче повторює алгоритм, поки не буде знайдено принаймні один дільник $N=15$. Ми використаємо результати запуску на обладнанні вище, щоб вгадати фазу та відповідний дільник на кожній ітерації.

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

Обговорення

Споріднені роботи

У цьому розділі ми обговорюємо інші важливі роботи, які продемонстрували алгоритм Шора на реальному обладнанні.

Основоположна робота [3] від IBM® продемонструвала алгоритм Шора вперше, факторизувавши число 15 на його прості дільники 3 і 5, використовуючи семикубітний квантовий комп'ютер на основі ядерного магнітного резонансу (NMR). Інший експеримент [4] факторизував 15, використовуючи фотонні кубіти. Використовуючи один кубіт, що багаторазово використовувався, та кодуючи робочий регістр у станах вищої розмірності, дослідники зменшили необхідну кількість кубітів до однієї третини від стандартного протоколу, використовуючи компільований двофотонний алгоритм. Значною роботою в демонстрації алгоритму Шора є [5], яка використовує техніку ітеративної оцінки фази Кітаєва [8] для зменшення вимог до кубітів алгоритму. Автори використовували сім контрольних кубітів і чотири кубіти кешу разом із реалізацією модульних множників. Ця реалізація, однак, вимагає вимірювань у середині схеми з операціями прямого зв'язку та повторного використання кубітів з операціями скидання. Ця демонстрація була виконана на іонно-пастковому квантовому комп'ютері.

Більш пізня робота [6] зосередилась на факторизації 15, 21 і 35 на обладнанні IBM Quantum®. Подібно до попередньої роботи, дослідники використовували компільовану версію алгоритму, яка застосовувала напівкласичне квантове перетворення Фур'є, як запропонував Кітаєв, щоб мінімізувати кількість фізичних кубітів і вентилів. Найновіша робота [7] також виконала концептуальну демонстрацію факторизації цілого числа 21. Ця демонстрація також включала використання компільованої версії процедури квантової оцінки фази та була заснована на попередній демонстрації [4]. Автори пішли далі цієї роботи, використовуючи конфігурацію наближених вентилів Тоффолі із залишковими фазовими зсувами. Алгоритм було реалізовано на квантових процесорах IBM, використовуючи лише п'ять кубітів, і заплутаність між контрольними та реєстровими кубітами була успішно підтверджена.

Масштабування алгоритму

Ми зазначаємо, що шифрування RSA зазвичай передбачає розміри ключів порядку від 2048 до 4096 біт. Спроба факторизувати 2048-бітне число за допомогою алгоритму Шора призведе до квантової схеми з мільйонами кубітів, включаючи витрати на виправлення помилок, і глибину схеми порядку мільярда, що виходить за межі можливостей сучасного квантового обладнання для виконання. Тому алгоритм Шора вимагатиме або оптимізованих методів побудови схем, або надійного виправлення квантових помилок, щоб бути практично життєздатним для зламу сучасних криптографічних систем. Ми відсилаємо Вас до [9] для більш детального обговорення оцінки ресурсів для алгоритму Шора.

Завдання

Вітаємо з завершенням підручника! Зараз чудовий час, щоб перевірити своє розуміння. Чи могли б Ви спробувати побудувати схему для факторизації 21? Ви можете вибрати $a$ на свій власний вибір. Вам потрібно буде вирішити щодо бітової точності алгоритму, щоб вибрати кількість кубітів, і вам потрібно буде спроектувати оператори модульного піднесення до степеня $M_a$. Ми рекомендуємо Вам спробувати це самостійно, а потім прочитати про методології, показані на рис. 9 з [6] і рис. 2 з [7].

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

    # Your code here
    pass

Посилання

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.