خوارزمية شور

تقدير الاستخدام: ثلاث ثوانٍ على معالج Eagle r3 (ملاحظة: هذا تقدير فحسب. قد يتفاوت وقت التشغيل الفعلي.)

خوارزمية شور، التي طوّرها Peter Shor عام 1994، هي خوارزمية كمية رائدة لتحليل الأعداد الصحيحة إلى عواملها في زمن متعدد الحدود. تكمن أهميتها في قدرتها على تحليل الأعداد الصحيحة الكبيرة بصورة أسرع بشكل أسي من أي خوارزمية كلاسيكية معروفة، مما يُشكّل تهديداً لأمن أنظمة التشفير المستخدمة على نطاق واسع كنظام RSA، الذي يعتمد على صعوبة تحليل الأعداد الكبيرة. وعبر حل هذه المسألة بكفاءة على حاسوب كمي قوي بما يكفي، يمكن لخوارزمية شور أن تُحدث ثورة في مجالات كالتشفير وأمن المعلومات والرياضيات الحسابية، مما يُبرز القوة التحويلية للحوسبة الكمية.

يركز هذا البرنامج التعليمي على توضيح خوارزمية شور عن طريق تحليل العدد 15 على حاسوب كمي.

أولاً، نُعرّف مسألة إيجاد الرتبة ونبني الدوائر المقابلة لها من بروتوكول تقدير الطور الكمي. بعد ذلك، نُشغّل دوائر إيجاد الرتبة على عتاد حقيقي باستخدام دوائر ذات أقصر عمق يمكن تحويلها. يُكمل القسم الأخير خوارزمية شور بربط مسألة إيجاد الرتبة بالتحليل إلى عوامل.

نختتم البرنامج التعليمي بمناقشة حول التوضيحات الأخرى لخوارزمية شور على عتاد حقيقي، مع التركيز على كل من التطبيقات العامة وتلك المصمَّمة خصيصاً لتحليل أعداد محددة كـ 15 و21. ملاحظة: يركز هذا البرنامج التعليمي أكثر على التطبيق والتوضيح العملي للدوائر المتعلقة بخوارزمية شور. للاطلاع على مرجع تعليمي متعمق حول المادة، يُرجى الرجوع إلى دورة أساسيات الخوارزميات الكمية للدكتور John Watrous، والأوراق البحثية في قسم المراجع.

المتطلبات

قبل البدء بهذا البرنامج التعليمي، تأكد من تثبيت ما يلي:

• Qiskit SDK الإصدار v2.0 أو أحدث، مع دعم التصور البياني
• Qiskit Runtime الإصدار v0.40 أو أحدث (pip install 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$ كيوبت، إلى جانب دائرة كمية أحادية الاتجاه (unitary) تعمل على $n$ كيوبت. ويُفترض أن $\ket{\psi}$ هي متجه ذاتي للمصفوفة الأحادية $U$ التي تصف عمل الدائرة، وهدفنا هو حساب القيمة الذاتية $\lambda = e^{2 \pi i \theta}$ أو تقديرها تقريبياً. بعبارة أخرى، يجب أن تُخرج الدائرة تقريباً للعدد $\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$. المكوّن الأول في دائرة تقدير الطور هو عمليات التحكم الأحادية (controlled-unitary) المسؤولة عن تنفيذ الرفس الطوري (phase kickback) إلى كيوبت التحكم المقابل. هذه العمليات الأحادية المتحكَّم بها تُرفع إلى قوى وفقاً لموضع كيوبت التحكم، بدءاً من البت الأقل أهمية حتى البت الأكثر أهمية. بما أن $\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$ ليس بتكرار دائرة $M_a$ عدد $k$ مرة، بل بدلاً من ذلك بحساب $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}$ بالملاحظة، يمكننا أن نرى أن حالات الأساس تُعاد ترتيبها، لذا لدينا مصفوفة تبادل. يمكننا بناء هذه العملية على أربعة كيوبتات باستخدام بوابات SWAP. في الكود التالي، نبني عمليتَي $M_2$ والـ controlled-$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}$

لذلك، يمكن بناء هذا التبادل بعملية 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)

ستُفكَّك البوابات التي تعمل على أكثر من كيوبتين إلى بوابات ثنائية الكيوبت.

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 Circuit، التي تستخدم بوابات Hadamard على كل كيوبت، وسلسلة من بوابات controlled-U1 (أو Z، حسب الطور)، وطبقة من بوابات 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)

لاحظ أننا حذفنا عمليات الأس المعياري المتحكَّم بها من كيوبتات التحكم المتبقية لأن $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 primitives

أولاً، نناقش ما يمكن الحصول عليه نظرياً إذا شغّلنا هذه الدائرة على محاكٍ مثالي. فيما يلي مجموعة من نتائج المحاكاة للدائرة أعلاه باستخدام 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$، ويتضح أن خوارزمية Shor لإيجاد الرتبة قد تفشل في بعض الأحيان. هذه النتائج الخاطئة سببها أن $k = 0$، أو أن $k$ و $r$ ليسا أوليَّين نسبياً - وبدلاً من $r$، نحصل على عامل من عوامل $r$. الحل الأبسط لذلك هو تكرار التجربة حتى نحصل على نتيجة مُرضية لـ $r$. حتى الآن، نفّذنا مسألة إيجاد الرتبة لـ $N=15$ مع $a=2$ باستخدام دائرة تقدير الطور على محاكٍ. الخطوة الأخيرة من خوارزمية Shor هي ربط مسألة إيجاد الرتبة بمسألة تحليل الأعداد الصحيحة إلى عواملها الأولية. هذا الجزء الأخير من الخوارزمية كلاسيكي بحت ويمكن حله على حاسوب كلاسيكي بعد الحصول على قياسات الطور من الحاسوب الكمومي. لذلك، نؤجل هذا الجزء الأخير من الخوارزمية حتى بعد أن نُوضّح كيفية تشغيل دائرة إيجاد الرتبة على أجهزة حقيقية.

تشغيلات على الأجهزة الحقيقية

يمكننا الآن تشغيل دائرة إيجاد الرتبة التي ترجمناها سابقاً لـ ibm_marrakesh. نلجأ هنا إلى dynamical decoupling (DD) لقمع الأخطاء، وgate twirling لأغراض تخفيف الأخطاء. تتضمن تقنية DD تطبيق تسلسلات من النبضات التحكمية المضبوطة بدقة على الجهاز الكمومي، مما يُعيّر التفاعلات البيئية غير المرغوبة وظاهرة decoherence بصورة فعّالة. أما gate twirling، فتُعشوئ بوابات كمومية محددة لتحويل الأخطاء المتماسكة إلى أخطاء Pauli، التي تتراكم خطياً لا تربيعياً. كثيراً ما تُجمع هاتان التقنيتان معاً لتعزيز التماسك والدقة في العمليات الحسابية الكمومية.

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

كما نرى، حصلنا على السلاسل البتية ذاتها بأعلى العدادات. نظراً لوجود ضوضاء في الأجهزة الكمومية، يحدث بعض التسرب إلى سلاسل بتية أخرى، يمكننا تصفيتها إحصائياً.

# 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: المعالجة اللاحقة وإعادة النتيجة بالتنسيق الكلاسيكي المطلوب

تحليل الأعداد الصحيحة إلى عواملها الأولية

ناقشنا حتى الآن كيفية تنفيذ مسألة إيجاد الرتبة باستخدام دائرة تقدير الطور. الآن، نربط مسألة إيجاد الرتبة بتحليل الأعداد الصحيحة إلى عواملها الأولية، مما يُكمل خوارزمية Shor. لاحظ أن هذا الجزء من الخوارزمية كلاسيكي. نُوضّح ذلك الآن باستخدام مثالنا $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 ***

النقاش

أعمال ذات صلة

في هذا القسم، نناقش أعمالاً بارزة أخرى أثبتت خوارزمية Shor على أجهزة حقيقية.

العمل الرائد [3] من IBM® أثبت خوارزمية Shor لأول مرة، إذ حلّل العدد 15 إلى عامليه الأوليين 3 و5 باستخدام حاسوب كمومي بالرنين المغناطيسي النووي (NMR) من سبعة قيوبتات. كما حلّل تجربة أخرى [4] العدد 15 باستخدام قيوبتات ضوئية. من خلال استخدام قيوبت واحد يُعاد تدويره مرات عدة وترميز سجل العمل في حالات عالية الأبعاد، تمكّن الباحثون من تقليص عدد القيوبتات المطلوبة إلى ثلث العدد المطلوب في البروتوكول القياسي، مستخدمين خوارزمية مُجمَّعة بفوتونَين. ورقة بحثية مهمة في إثبات خوارزمية Shor هي [5]، التي تستخدم تقنية تقدير الطور التكرارية لـ Kitaev [8] لتقليص متطلبات القيوبتات في الخوارزمية. استخدم المؤلفون سبعة قيوبتات تحكم وأربعة قيوبتات ذاكرة مؤقتة، إلى جانب تنفيذ مضاعِفات معيارية. غير أن هذا التنفيذ يستلزم قياسات داخل الدائرة مع عمليات تغذية أمامية وإعادة تدوير القيوبتات مع عمليات إعادة ضبط. أُجري هذا الإثبات على حاسوب كمومي بفخ الأيونات.

ركّزت أعمال أحدث [6] على تحليل أعداد 15 و21 و35 على أجهزة IBM Quantum®. وعلى غرار الأعمال السابقة، استخدم الباحثون نسخة مُجمَّعة من الخوارزمية تعتمد تحويل فورييه الكمومي شبه الكلاسيكي كما اقترحه Kitaev، وذلك لتقليص عدد القيوبتات المادية والبوابات المطلوبة. وأجرى عمل أحدث [7] أيضاً إثباتاً تجريبياً مبدئياً لتحليل العدد الصحيح 21. اشتمل هذا الإثبات كذلك على استخدام نسخة مُجمَّعة من روتين تقدير الطور الكمومي، وبنى على الإثبات السابق [4]. تجاوز المؤلفون ذلك العمل باستخدام تهيئة من بوابات Toffoli التقريبية ذات إزاحات طور متبقية. نُفّذت الخوارزمية على معالجات IBM الكمومية باستخدام خمسة قيوبتات فحسب، وتم التحقق بنجاح من وجود التشابك بين قيوبتات التحكم وقيوبتات السجل.

قابلية توسع الخوارزمية

نشير إلى أن تشفير RSA يتضمن عادةً أحجام مفاتيح تتراوح بين 2048 و4096 بتاً. ستُسفر محاولة تحليل عدد مكوّن من 2048 بتاً باستخدام خوارزمية Shor عن دائرة كمومية تحتوي على ملايين القيوبتات، بما في ذلك التكلفة الإضافية لتصحيح الأخطاء، وعمق دائرة يبلغ مليار بوابة، وهو ما يتجاوز حدود الأجهزة الكمومية الحالية. لذلك، ستحتاج خوارزمية Shor إما إلى أساليب محسّنة لبناء الدائرة أو إلى تصحيح أخطاء كمومي موثوق كي تكون ذات جدوى عملية في كسر أنظمة التشفير الحديثة. نُحيلك إلى [9] للاطلاع على نقاش أكثر تفصيلاً حول تقدير الموارد لخوارزمية Shor.

التحدي

تهانينا على إتمام هذا البرنامج التعليمي! الآن وقت مثالي لاختبار مدى فهمك. هل يمكنك محاولة بناء الدائرة لتحليل العدد 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).